Hydrogeological modeling of the Roussillon coastal aquifer (France): stochastic inversion and analysis of future stresses

Global climate change-induced stresses on coastal water resources include water use restrictions, saline intrusions, and permanently modifying or damaging regional resources. Groundwater in coastal regions is often the only freshwater resource available, so an in-depth understanding of the aquifer, and the aquifer’s response to climate change, is essential for decision-makers. In this study, we focus on the coastal aquifer of Roussillon (southern France) by developing and investigating a steady-state groundwater flow model (MODFLOW 6) and calibrated with PEST++ on a Python interface (FloPy and PyEmu). Model input and boundary conditions are constrained by various scenarios of climate projections by 2080, with model results predicting the aquifer’s response (and associated uncertainty) to these external forcings. Using simple assumptions of intrusion estimates, model results highlight both strong climatic and anthropogenic impacts on the water table. These include aquifer drawdowns reaching several meters locally, and the seawater interface advancing locally several hundred meters inland and rising by several meters. Intrusions of this magnitude risk endangering exploited water wells and their sustainability. Our results demonstrate the critical importance of properly characterizing the geology and its heterogeneity for understanding aquifers at risk because poor predictions may lead to inappropriate decisions, putting critical resources at risk, particularly in coastal environments.


Introduction
It is widely recognized that global warming and its consequences on freshwater resources will be impactful and irreversible (Cisneros et al. 2014;Masson-Delmotte et al. 2021). Increasing trends in temperatures, heat waves and drought in the Mediterranean region already exceed global trends (Cramer et al. 2018), particularly after the 1980s (Lionello and Scarascia 2018), while rising sea levels will impact coastal landscapes and resources (Church et al. 2013).
Coastal aquifers are particularly vulnerable to these issues, with likely decreases in water quantity related to a warmer and drier future, and increasing vulnerability of coastal aquifers to seawater intrusions (Jiao and Post 2019). This inevitably links to a decrease in the availability and quantity of water resources suitable for consumption. Therefore, it is essential to develop a deep understanding of how these forcing factors influence the behavior of coastal aquifers, and to guide the decisions of aquifer managers in achieving sustainable exploitation of water resources in the future.
The study of (future) climatic and anthropogenic impacts on coastal aquifers is not novel (Pinder and Cooper 1970;Oude Essink 1996;Bear et al. 1999;Barlow 2003;Meyer et al. 2018) and is often based on the numerical modeling of density-dependent flow processes (e.g. Carneiro et al. 2010;Comte et al. 2014;Dunlop et al. 2019). Many of these studies consider the case of sedimentary aquifers whose geological structures can be described using numerical parameterization (e.g. Giambastiani et al. 2007;Masterson and Garabedian 2007;Meyer et al. 2018).
Numerical models represent, therefore, an indispensable tool for studying the complexity of coastal systems. Constructing numerical models requires large amounts of data such as topographic maps, geological models, hydrogeological parameters, and boundary conditions, and which are sometimes difficult to obtain. In addition, different types of models (e.g., agronomic and climatic Haj-Amor et al. 2020) are required to couple the influence of other complex processes such as agricultural land use and anthropogenic irrigation strategies groundwater.
One of the key factors for robust simulation of groundwater-seawater interactions is the accurate description of the internal variability of the hydrodynamic properties of the aquifer. Analytical solutions are generally limited to specific cases (Oude Essink et al. 2010), so offer little help in cases where hydrodynamic properties are linked to the spatial organization of lithological facies (Dall'Alba et al. 2020;Micallef et al. 2020), requiring more versatile numerical models (Jiao and Post 2019). Abdoulhalik et al. (2020) carried out several laboratory experiments investigating salt intrusions for the case of a single-layer heterogeneity, and used to calibrate their numerical model SEAWAT (Langevin et al. 2007). They successfully demonstrated the importance of an impermeable layer on the salt wedge and the sensitivity of heterogeneous aquifers to up-coning intrusions. The internal variability of the hydrodynamic properties within aquifers is complex and a major source of uncertainty (Kerrou and Renard 2010; Michael et al. 2016). This uncertainty can be reduced numerically with deterministic approaches like direct parameter estimation software, through calibration and inversion procedures, Doherty et al. 2010;Meyer et al. 2018) or with stochastic approaches (Feyen and Caers 2006;Alcolea et al. 2009).
This work aims to explore the sensitivity of a Mediterranean coastal aquifer, the Roussillon aquifer (France), to future climatic and anthropogenic changes. This aquifer experienced continuous water table declines (up to several meters locally) over the last decades (Chabart 1996;Aunay 2007;Caballero and Ladouche 2015). Global demand for freshwater is expected to increase, or at best stabilize, due to agricultural activity, changes in land use, and ongoing tourism development. The risk of seawater intrusion in the aquifer also increases with rising sea levels (Galassi and Spada 2014) and changing climatic conditions (Terray and Boé 2013). This is particularly acute during high-stress periods in the summer. These saline intrusions could lead to local and temporal limitations of groundwater pumping.
To answer these questions, we constructed a 3D steadystate groundwater model with MODFLOW 6 (Langevin et al. 2017b). The model permits to investigate the impact of these changes on water resources, and identify the dominant parameters using PEST++ . To quantify the uncertainty, we generated 150 different stochastic realizations (derived from the calibrated model) using the first-order second-moment analysis method. We then use this ensemble of models with different scenarios to estimate the possible range of impact of different climatic (decrease of the recharge and sea-level rise) and anthropogenic factors (increased water withdrawals) on the risk of seawater intrusion and resource sustainability by the year 2080. The models are constructed with open-source and free software (MODFLOW and PEST) using a Python interface [FloPy (Bakker et al. 2016) and PyEMU (White et al. 2016)]. This approach facilitates the integration of the groundwater model with other components such as, for example, river discharge or water use in economic models.
In terms of novelties, this study allows to better understand how the Roussillon aquifer may react to new external stresses including the influence of climate change. As compared to previous modeling work in the region (Chabart 1996), the construction of a new model allows integrating new elements such as the importance of recharge from canals and provides a quantification of the uncertainty of the forecasts. This work opens also new perspectives for coupling the groundwater model with a novel economic model to aim for a better management of the water resource.

Description of study area
The study site lies in the Roussillon the plain on the western border of the Mediterranean Sea, in the Pyrénées-Orientales region in southern France (Fig. 1,42.65 N 2.97 E). The surface of the aquifer measures nearly 900 km 2 , a quarter is below 10 m.a.s.l. The relief slopes slightly decrease from west to east, with maximum heights of about 200 m.a.s.l. on the western limit.
The Mediterranean climate includes hot and dry summers and mild winters, with an average minimal temperature in the city of Perpignan of 11.6 • C and an average maximal temperature of 20.5, while the average rainfall is 578 mm (1991-2020. From a hydrological point of view, the plain is crossed by 4 perennial rivers (the Agly, the Têt, the Réart, and the Tech; see Fig. 1), most of which originate from the surrounding massifs. The catchment areas of these rivers range between 700 and 1300 km 2 , except for the Réart river (160 km 2 ). Two other minor and temporary rivers augment the four main rivers, with the Boulès a main tributary of the Têt, and the Canterrane a main ributary of the Réart. It is also important to note the presence of 2 large lagoons, the Canet (7.82 km 2 ) and Salses-Leucate (54 km 2 ) lagoons.
Geology The sedimentary basin of the Roussillon plain is the result of the filling of the Gulf of Lion (south of France, Mediterranean region) after the erosion process that affected the Mediterranean basin during the Messinian crisis Page 3 of 20 201 (Clauzon et al. 2015). This filling occurred following a Gilbert delta structure (Gilbert 1890) in a prograding-aggrading system during the Pliocene (Aunay et al. 2004;Duvail 2008;Clauzon et al. 2015). This sedimentological setting results in a superposition of many prisms composed of clay in the distal parts (bottomsets) and sands in the proximal parts (foresets). However, the area of continental deposits and alluvial fans (topsets) is very heterogeneous with the development of various channels (braids, meanders) and floodplains. A complete description of the geological depositional history can be found in Duvail et al. (2001). Figure 2 shows a conceptual 3D block of the geology and the aquifer.
This sedimentation process was active for several million years and resulted in a lithological spatial heterogeneity that is complex. The Roussillon plain is thus a multilayered basin where three different aquifer units can be identified: the Quaternary (Qt), the Continental Pliocene (CP) and the Sandy Marine Pliocene (SMP). Qt generally follows the riverbeds. It is mainly made up of terraces, alluvial channels and flood plains. These deposits extend in the offshore domain but data are scarce and marine Quaternary is thus poorly documented. The CP is known to be highly heterogeneous with sandy deposits of various sizes alternating with more or less clayey levels over thicknesses of several meters. Some attempts to model the CP heterogeneity with Multiple Points Statistics (MPS) have been presented earlier (Dall'Alba et al. 2020). The SMP is essentially composed of marine sand deposits that extend over distances of hundreds  of meters to kilometers and thicknesses of a few meters to tens of meters and is viewed as mostly homogeneous. At the bottom of the SMP, the Clayey Marine Pliocene (CMP) unit represented in dark/light blue in Fig. 2 is totally impermeable and represents the lower impervious Roussillon plain aquifer's boundary. Hydrogeology The regional Roussillon aquifer can be separated into two major layered water bodies: the Quaternary aquifer (QA) and the Pliocene aquifer (PA), composed of CP and SMP units. If we detail the units of the PA, the top formation, CP, can be described regarding its hydraulic properties as a formation composed mostly of highly permeable, heterogeneously distributed paleochannels in a relatively impermeable to semi-permeable matrix. Seismic campaigns and hydrogeological tests indicate that these paleochannels are likely to be spatially well connected, creating an important available groundwater resource (Dall'Alba et al. 2020). The SMP bottom formation is generally considered to be rather homogeneous and highly permeable (Chabart 1996; Aunay et al. 2004).
The PA appears to be confined over most of the plain due to the fact that its permeable layers are generally embedded into clay rich formations having a low vertical permeability (Chabart 1996; Aunay et al. 2004). Locally, the PA outcrops over the western part of the plain and can thus be unconfined, but for simplicity, we assume they are confined. Table 1 shows the ranges and the arithmetic mean values for the hydraulic conductivity (K) of both aquifers. The values were taken from the public report P22b of the DEM'EAUX project (Dewandel et al. 2022).
From a hydrogeological perspective, the PA appears favorable for protecting from seawater intrusions for three reasons. First, the QA at the top and the transition between the two aquifers act as a natural shield, stopping or mitigating continuous inflows from the sea to the PA. Second, the hydraulic heads in the PA near the coast are higher than the topography, which results in groundwater flow toward the sea. Finally, the hydraulic gradient, between the two aquifers, generally induces flows from PA to QA on the coast, reducing the risk of indirect SI from QA to PA.
However, these protections are reduced and could potentially be problematic in a near future, because the hydraulic head in some area of the Pliocene aquifer underwent a nearly constant decrease over several decades, primarily due to pumping withdrawal (Chabart 1996;Aunay et al. 2006;HYDRIAD 2014b;Caballero and Ladouche 2015). The increase of the intake, mainly during summer, weakens the natural vertical flow and could locally reverse it, which could lead to local SI, either laterally or vertically from QA to PA.
Laterally, most of the limits of the study zone ( Fig. 1) are considered as no flow boundaries. However, some zones are known to contribute to the system. The karstic Corbières massif (northwest) is estimated to supply nearly 200 l/s to QA and PA (Ladouche and Dörflinger 2004). Other massifs, such as Aspres, located southwest of the area, are also suspected to be a supply zone for the aquifer (Chabart 1996), but the associated flows are unknown.
The water budget is not known with certainty. Important questions remain regarding the interactions between the aquifer and the rivers (HYDRIAD 2014b), with further uncertainty in withdrawals for irrigation and agricultural withdrawal. The most recent estimates, based on land use and agriculture type mapping (BRLI 2018b; Terrasson et al. 2014), show that the total amount of water entering (or leaving) the aquifer is greater than 117 m 3 /year (> 3.7 m 3 /s), providing an order of magnitude of the aquifer water budget.

Methods
In the following section, we provide a detailed description of the methods and tools used to create and identify parameters of the hydrogeological model of the Roussillon aquifer. We discuss the construction of the models, boundary conditions (BCs), and inversion strategies, as well as provide an overview of the different climatic scenarios that were investigated.

Model construction
The MODFLOW groundwater modeling software is widely used and recognized in groundwater studies (McDonald et al. 2003;Harbaugh 2005;Langevin et al. 2017a). In this paper, we use the MODFLOW 6 version and the opensource Python package FloPy (Bakker et al. 2016(Bakker et al. , 2020) to build the model, run it, and analyze the results. Figure 3 presents the extension of the model. The offshore limit is fixed at the limits of the known geological model (around 25 km from the seashore). The horizontal resolution of the model is set to 300 × 300 m inland and 600 × 300 m at sea. The difference in resolution allows refining the mesh around active groundwater flow dynamics, or better data availability. This results in a model composed of 182 rows ×183 columns and a total of 70,755 active cells. The upper layer of the model is set using a DEM of the area with an original resolution of 5 × 5 m, resampled to the resolution of the model. It is built using the DIS MODFLOW 6 package (Fig. 3). The model is divided vertically into five layers:

Discretization and parameters
• Qm: The Marine Quaternary (offshore only) is a thin 2 m interface layer laying on top of the offshore part of the Quaternary.
• Qt: Quaternary aquifer, its bottom corresponds to the base of the Quaternary and its top to the DEM on land and the base of Qm offshore. • Qt-CP interface: A semi-permeable layer that separates the Qt and CP aquifers (Vandenbeusch and Marchal 1980;Chabart 1996). • CP: The main Pliocene aquifer. Its bottom is defined by the CP surface from the geological model. • SMP: The marine sands from the Pliocene. Its bottom is defined by the SMP surface from the geological model.
The minimum cell thickness is fixed to 2 m to avoid numerical artifacts. Figure 3 shows a map of the 3D model mesh and the different formations. Conductivity (K) of these five layers are initially assumed to be homogeneous and isotropic with a fixed value determined by prior information derived from the results of existing pumping tests ( Table 2). The impermeable layer between QA and PA was defined with a thickness of 2 m and a conductivity of 0.25×10 −8 m/s (following Chabart 1996). These homogeneous layers were then manually and locally calibrated (using a least-square optimization) and adjusted before performing an automatic parameter identification with PEST ++.

Boundary conditions
A constant head boundary condition (BC) is prescribed in the area of the model covered by the sea using the constant head package (CHD). Since no data are available offshore and since we are mainly interested in the risk of intrusions near the coast, we kept the offshore part of the model very simple. A unique constant head value of 5 × 0.025 = 0.0125 ms (equivalent freshwater head) is prescribed for the entire offshore area. This choice allows us to minimize unrealistic groundwater flow circulations in the offshore part of the model due to the complex bathymetry.
Environmental Earth Sciences (2023)  The rivers (Agly, Têt, Boulès, Réart, Canterrane, and Tech) are included using a Cauchy BC using the RIV package. The river stage is defined according to the DEM and the river bottom is set to 1.5 m below the river stage (based on existing data and expert knowledge). The initial conductance parameter for all river beds is fixed to 0.005 m 2 /s based on the conductance equation combining the geometry (width, thickness and length) and the hydraulic conductivity of the riverbed (see Langevin et al. 2017b).
Lagoons are implemented using the MODFLOW 6 Drain package (DRN), as we assume that these water bodies act only as outlets for the aquifer. This is based on a technical report of the DEM'EAUX project that indicates that no flow from the lagoons to the aquifer can be observed (Caballero and DBLadouche 2022). A fixed altitude is set at 0.1 m with a conductance value of 1 m 2 /s. Along the lateral boundaries, no flux BC is fixed except in two places (green lines, Fig. 1).
• A Neumann BC is prescribed near the Corbières massif (north of the Roussillon plain) with a flux of 200 l/s homogeneously distributed over the cells in the PA (Ladouche and Dörflinger 2004). • Another Neumann BC is prescribed in the limestone massif of Aspres, but the magnitude of the flux is unknown. There are suspicions about this flux (Chabart 1996) but no detailed field investigation has been done in this area. In our study, an initial estimate of 300 l/s is used to constrain the model to respect head observations. The flux is distributed homogeneously on the cells of both aquifers (QA and PA).
Recharge The recharge is made up of two parts; one accounts for precipitations and the other accounts for the water that comes from irrigation of agricultural parcels and irrigation channels.
The recharge from precipitation is estimated using simple soil water budget approaches to calculate effective precipitation (Lanini et al. 2019; Healy 2010) combined with an infiltration coefficient related to the index of development and persistence of river networks (IDPR Mardhel et al. 2021) and a cartographic index that provides a qualitative estimate of infiltration versus runoff (Martinsen et al. 2022). As the initial data (SAFRAN meteorological data reanalysis, Vidal et al. 2010) are coarse (grid of 8 km by 8 km), we applied a moving average to obtain the final recharge to avoid spatial artifacts.
In the plain, irrigation can be decomposed into two components: one from canals and one from aquifer withdrawals by pumping. We assumed that aquifer irrigation has little losses because farmers use drip irrigation most often in this case (HYDRIAD 2013). On the contrary, the seepage from agricultural channels has been shown to be very important (2-5 m 3 /s, BRLI 2018b). It locally exceeds several times the "natural" recharge (Adombi 2019).
The recharge from the channels is estimated by zones and is based on the water needs (depending on the type of crop and the local climatic conditions), the irrigation rate of the channels (irrigation water proportion of channel origin) and a multiplicative coefficient that represents the efficiency of the channels estimated by previous technical studies (BRLI 2018a, b). It represents the ratio between the recharge from the canals and the water needs. It is set to two for the simulation, based on the estimations made by BRLI (2018b). Practically, it indicates that if we consider an agricultural area, solely irrigated by channels, that needs 10.000 m 3 of water per year, a total inflow of 20.000 m 3 per year will be established.
Withdrawals The main withdrawals in the plain can be decomposed into two parts, one for agricultural purposes and the other for drinking water supply.
The agricultural withdrawals are not known accurately and, therefore, must be estimated. We use a method similar to that for irrigation recharge, but assuming a low level of leaks ( ∼ 20%). The water needs and the well irrigation rate (irrigation water proportion of groundwater origin) were defined by zones (BRLI 2018a), combined with leakage coefficients, leading to rough estimates of withdrawals. Finally, these withdrawal values are redistributed between the different aquifers (QA and PA) using estimates made by HYDRIAD (2013) that indicate, by areas, the proportion of water withdrawn from each aquifer. All these BCs are implemented with the RCH package.
Finally, we use groundwater withdrawal data regarding drinking water supply and industrial use collected by the Rhône-Méditerrannée-Corse Water Agency (AERMC) and the Roussillon Plain Association for the Protection and Management of Groundwater (SMNPR). This database contains all annual registered withdrawals from 1987 to 2018. For the reference year for the model (2013), it identified 152 producing wells distributed over the study area. These withdrawals were implemented with the WEL package.

Parameter identification
This section describes the inversion process and how the model was tested. The hydrogeological model is inverted using the PEST++ software ) and the PyEMU interface (White et al. 2016).

Parameters
Hydraulic conductivity values (K) of the QA and PA layers have been parametrized using the Pilot Point method (PiPM) (Marsily et al. 1984;Certes and de Marsily 1991;RamaRao Page 7 of 20 201 et al. 1995;Fienen et al. 2009;Doherty et al. 2010). The 280 (118 in Qt and 162 in CP) Pilot Points (pp) are disposed uniformly in these two layers with a fixed spacing distance of 2400 m. The upper and lower limits, as well as the initial values, come from the literature (Chabart 1996, Aunay 2007, Dewandel et al. 2022. The other layers composing the aquifer (Marine Quaternary, Qt-CP interface, and SMP) have been set with fixed K values since they are considered much more homogeneous.
Regarding the rivers, their values of conductance coefficients are poorly known. However, field observations suggest that the exchanges with the aquifer are significant (HYDRIAD 2014b;BRLI 2018a). In view of the limited information available on conductance values and for simplicity, homogeneous values were used for each river and estimated during the inversion process.
Since the recharge from the agricultural parcels carries a large uncertainty, it was decided to estimate them as well during the inversion using a direct multiplicative coefficient of the initial recharge. This coefficient is allowed to vary between 0.8 and 2 for each irrigation zone.
Around 300 parameters have to be identified. They are listed in Table 2, with their associated constraints, used to control the range of acceptable values. A regularization term is added for each parameter following the preferred value approach. This helps the parameters to not deviate too much from a specified value (here the initial value) (Doherty et al. 2010).

Observations
A total of 123 observation points of the groundwater level are available (75 in QA and 48 in PA). Preliminary tests showed difficulties to reach the inversion with these scarce data sets only. To address this issue, observation points were added based on available piezometric maps (2013 campaign for the QA levels and 1990 campaign for the PA). The 1990 map is used for the Pliocene levels because the most recent one (2013) has very few data and shows artifacts that are questionable. However, since some of the observations in PA are separated by 23 years, it is likely that they represent a different hydrological setting. Aware of this issue, we have decided to give a smaller importance to the observations coming from the piezometric maps (even for the ones from the QA piezometric map).
Finally, to assess the performance of the calibrated models, we used two different criteria: the Root Mean Squared Error (RMSE) and the Mean Absolute Error (MAE). We only considered the real 123 observation points with equal weights. The observations from the piezometric maps were left apart to ensure that the calibrated model has been improved on the real data, even using another source of observations. Moreover, it gives a better representation of the errors than if all the observations were included.

Inversion and uncertainty estimation
The identification of parameters was carried out using PEST++ software ). PEST and PEST++ minimize an objective function defined as the combination of two components (Doherty 2015): where Φ m (Eq. 2) is the measurement misfit component, Φ r is the regularization component, and is the regularization factor. PEST++ automatically adjusts during inversion by aiming at a certain value of the user-defined objective function (PHIMLIM control variable).
Φ m measures the discrepancies between the actual measurement data and the model output and is expressed as where for a specific observation point i, w i is the weight associated with the observation, h o i represents the head measurement and h c i the simulated head. n is the number of observation points. In our case, the weights are set to 1 for the actual measurements and to 0.5 for the head estimations derived from the piezometric maps.
Φ r measures the differences between the estimated parameters and their preferred values or conditions (Doherty 2015). Generally, it is computed as a weighted least squares of these differences. In our case, the weights are identical within parameters groups and automatically defined by PyEMU to give a relatively similar importance to each parameter group. (for example, all pilot points have the same weight, but it is different from the weight of the irrigation zone).
The objective function is minimized in a deterministic manner using the Levenberg-Marquardt algorithm to obtain the parameter set that best fits the data and regularization criteria. The termination criteria and the PEST control variables were defined using default and/or recommended values for typical parameter identification problems (White et al. 2016;Doherty 2018).
In addition, PEST++ incorporates tools to quantify the parametric uncertainty. We used the first-order secondmoment analysis method (FOSM) available in PEST++. It assumes that the relation between parameters and outputs of the model can be approximated by a linear expression that is locally valid around the optimal value obtained by minimizing the objective function. It also assumes that the noise of the measurements and the uncertainty can be represented reasonably well with a multi-Gaussian distribution. Using Environmental Earth Sciences (2023)  these assumptions, it is possible to estimate a posterior probability distribution for the parameters and generate samples from that distribution. In practice once calibrated, we generated 150 parameter sets and used them to build 150 models for uncertainty estimation. To represent the final uncertainty on the different predictions, we post-processed the simulations and computed the standard deviation or percentile intervals of the predictions.

Climatic sensitivity
After estimating the hydraulic parameters of the models, the next step of the study is to study the sensitivity of the water resource under different climatic scenarios. These scenarios are then applied on the set of models obtained from the uncertainty analysis. We define three plausible scenarios for the year 2080: • Scenario 1: decrease of the recharge; • Scenario 2: decrease of the recharge + sea-level rise of 35 cm; • Scenario 3: decrease of the recharge + sea-level rise of 35 cm + water extraction increase of 20%.
Reduction of recharge from precipitation is based on the expected evolution of effective rainfall by 2080 for the Roussillon region obtained using an ensemble of simulations based on the projections of five global climate models and two downscaling methods for the RCP8.5 emission scenario (Lanini et al. 2019). These projections highlight a significant decrease in annual recharge due to both a decrease in precipitation and an increase in temperature (Terray and Boé 2013;Caballero and Ladouche 2015). The global reduction in the plain is estimated to be around 30%. Locally, no area appears to be affected differently, the reduction is quite homogeneous. We can point out that the southern to southeastern part of the area is likely to have very little recharge by 2080 (around 40 mm/year), with probable implications on groundwater levels and interactions with rivers. The reduction of the recharge is implemented in the model by reducing the effective rainfall which, after the application of the infiltration coefficient, corresponds to the natural recharge in the models. The rise in sea level is also taken into account in the analysis and estimated based on the calculations of Galassi and Spada (2014). According to their estimates, an increase of about 30-40 cm is expected for 2080. Therefore, we took a value of 35 cm.
Regarding the withdrawals for drinking water and industrial usages, their evolution is hard to predict, as not only the population and demand will change, but also restriction policies on water use are often subject to modifications. To test, however, its potential influence on the system, an increase in withdrawals was included in Scenario 3. It corresponds to a situation in which the population will continue to grow and, therefore, the water requirements will increase significantly. Increased temperature will lead to increased evapotranspiration and will require increased irrigation to compensate for this deficit. All of this will necessarily have an impact on withdrawals, which will have to increase to meet these needs. In this scenario, the 20% increase in withdrawals is arbitrary, but it is considered a plausible assumption, as it stays within the range of withdrawal variations observed over the past decades (which have increased by 50% between 198250% between and 200450% between , Terrasson et al. 2014Aunay 2007), even if this upward trend is slowing down since the 2000s and is on the same order of magnitude as scenarios built on the territory for 2030.
The results of the climate sensitivity approach are interpreted in terms of head reduction and potential seawater intrusion. The position of the seawater interface is estimated using the Ghyben-Herzberg equation (Ghyben and Drabbe 1889;Herzberg 1901) from the heads computed with the model in the QA. The altitude of the interface in each cell is then calculated as follows: where h s [m] is the altitude of the interface, h Q a [m] are the heads in the QA and BC sea [m] is the BC applied in the sea. This condition increases with increasing sea level and must be taken into account for forecasts. The use of this simple relation greatly reduces the complexity of the models and facilitates uncertainty analysis, but it can only be considered as a first approximation. However, it will permit to couple this model in the near future with an economic model and carry out a coupled hydro-economic optimization that would be much more difficult to pursue using a fully densitydependent transport model due to a lack of computational resources and consequently run-time issues.

Parameter identification
The parameter identification took nearly 3000 model runs and 10 PEST++ iterations to obtain a final solution that respects (at best) both field measurements and a priori knowledge.
Parameters Figure 4 shows the estimated K for both aquifers after the inversion. Strong variations can be observed such as a low conductivity zone in the southern part of the QA or conversely, a highly transmissive zone in the north. For the PA, we cannot determine distinct zones but the estimation show a rather large variability with values ranging from -6 to − 3.5 log 10 [ m s ] . Generally, K is higher in the PA than in the QA, which is consistent with measurements and knowledge about these aquifers (Chabart 1996;Duvail et al. 2001;Aunay 2007;HYDRIAD 2014b). The uncertainty is also quite high, especially near the sea, where the hydraulic conductivity has less influence on the heads. It should be noted that the influence of the pilot points on the estimation are visible on the K maps for both the PA and the QA layer, which could indicate that the number of pp was not sufficient to capture the rather complex geology of this aquifer.
The K values for the other layers are shown in Fig. 5G-I. Their estimations remained close to their initial values, except for the SMP where the median is around one order of magnitude smaller, indicating than the SMP could be less permeable than initially assumed. The results of the inversion show that the available data do not constrain well the conductance values for the rivers (Fig. 5A-F). The remaining uncertainty is large and the estimated values are often reaching the limits (Canterrane, Tech and Boulès).
Concerning the agricultural recharge (Fig. 6), it has been increased in the Têt and Tech upper valleys where the multiplication coefficient reach the fixed limit of 2. The fact that the optimization process always selects the maximum limit for these parameters may imply that the conceptual model and the BCs do not represent properly the reality. Nevertheless, these coefficients allow the recharge to reach 800 mm/ year in some areas (Fig. 6). These values were already suggested by Adombi (2019), which points out the importance of the recharge from canals to the system. However, the other zones have a value close to one, probably due to the fact that they have less impact on the total recharge.

Head distribution
The computed heads from the inverted model show a good agreement with the heads from the piezometric map ( Fig. 7A) for both aquifers. This is confirmed by the scatter plot (Fig. 7C). The best model have a MAE of 1.5 m and a RMSE of 2.3 m. The spatial distribution of these errors (Fig. 7D) show no trend and are thus well dispersed with no area badly simulated. However, we can point out that the simulated heads in the PA are often slightly too high, a fact confirmed with the mean of the error residuals that is around + 0.5 m (ideally it should be 0). This is not observed for QA.
To complete these results, the related uncertainty was derived from the 150 simulated models (Fig. 7B). Assuming that the heads follow a Gaussian distribution, we show the 95% confidence interval (i.e. two times the standard deviation ). This error can reach up to 10 m, especially near the inland boundaries, but was capped to 5 m for visualization. Generally, it ranges from 1 to 2 m and is relatively constant over the all domain.
The simulation in the offshore part of the model is not shown due to the fact that there's little to none effect of this model area on heads distribution.

Water budget
The computed water budget for the current state (Fig. 8A) indicates that the total flux through the system is relatively high (5.3 m 3 /s ≈ 165 Mio. m 3 /year). This estimation is greater than what was previously estimated (around 80-120 Mio. m 3 /year according to Chabart (1996) andHYDRIAD (2014b). This difference is mainly due to the considered area that is larger in this study as compared to Chabart (1996) but it is also due to a difference in the modeling approach that allows us to better capture the flux from the rivers, and others hydrological entities (lagoons, sea). Indeed these features have complex interactions with the aquifer and can act as a main outlet (e.g. the Têt) or entry as the Boulès (BOL) or both as the Canterrane (CANT). HYDRIAD (2014b) limited the infiltration where the PA is outcropping, assuming the formation to be impermeable on the top. They also integrate the impermeability of the urban area which in final lead to a recharge that was 20-30% smaller than in our study.
Finally, the uncertainty for the different water balance components was estimated using the 5th and 95th percentile over the 150 simulated models. The uncertainty is quite large (10%) due to multiple sources of uncertainty (which are the estimated parameters during the inversion process i.e. river conductance, hydraulic conductivity and the coefficient of recharge for the irrigation zones), without forgetting that some have not been taken into account such as the recharge.

Climatic forecasts
In this section, we present the climatic forecasts for 2080 in order to identify critical areas (regarding to the projections) and estimate the sensibility of the aquifer to climatic and anthropogenic response. Today scenario refers to the inverted model described in the previous section.

Head reduction and budgets
In this region, the two major consequences of climate change will be a significant reduction of the precipitations and a rise of the temperatures. The recharge will decrease and consequently the total water budget of the aquifer will be affected (Fig. 8). The model shows that the reduction of the rain and recharge ( −30% ) significantly impacts the exchange with the different rivers. For today, rivers release 0.8 ± 0.4 m 3 /s and capture 1.4 ± 0.3 m 3 /s of water (Fig. 8A), while these numbers would become 1.1 ± 0.4 m 3 /s and 0.9 ± 0.2 m 3 /s, respectively, for Scenario 3 (Fig. 8D). These values were obtained by summing all the inflows and outflows of the rivers for a specific scenario.
Under the selected scenarios, the rivers would contribute more to the aquifer and partly compensate the reduction of recharge from natural recharge. The rivers would also receive less water from the aquifer due to the recharge decrease. The simulated uncertainties about these results remain, however, large.
The model also shows (Fig. 9) that both aquifers would be affected in a similar manner, undergoing a strong decrease in groundwater levels. The heads reduction could reach up to 10 m with a strong anthropogenic influence. Figure 9 does not show scenario 2 as it influences the heads only near the coast by increasing them to around 0.5 m. The anthropogenic influence is clearly visible in the results of Scenario 3, which would lead to a water table drop of up to 5-10 ms depending on the location. The uncertainties are in order as those of the models for today (Fig. 7B). An interesting fact to note is that the heads near the shoreline are affected differently in both aquifers. In PA, we see a reduction of the heads up to 1 m, even in the offshore area. However, in QA, the heads rise due to the sea-level rise. Such changes can durably affect the hydrodynamics of the system and the exchanges between these aquifers. It is also possible to notice the effect of the rivers that minimize the climatic effects, as they cannot dry (the rivers are implemented with a head-dependent condition) and thus continuously contribute to the aquifer, leading to a small decrease of the heads.

Saltwater/freshwater interface
The modeled seawater/freshwater interface (Fig. 10) for Scenario 3 and today is shown as the median position along the 150 models. The rise of the sea level combined with the reduced recharge and increased extraction would lead to a progression of the interface inland. This movement would be limited or amplified depending on the steepness of the QA base, as shown in Fig. 10. The figure also shows that the uncertainty is large because the saltwater intrusion process (modeled here by the Ghyben-Herzberg relation) greatly amplifies the changes and uncertainties that affect the freshwater heads. Globally, the interface would move forward inland of several hundred meters, but depending on the steepness of the QA base, it could reach up to 1 km. Generally, the lower the slope, the further the salt wedge can penetrate inland. The interface could also rise vertically up to 10-15 m, according to the medians. Generally, the 10 percentile position for scenario 3 would be equal (or close) to the median position for today.
We summarized those results by producing probability maps for each scenario allowing us to assess the spatial evolution of the risk of seawater intrusion. Figure 11 shows the probability of observing that the interface moves north to a specific location in the different scenarios. Practically, it indicates where the interface is more likely to be present in the future and with which probability compared to today. These maps have been produced by comparing the ensemble of models for each scenario with the reference ("present"). Practically, for each model of each scenario, a map of the absence / presence of saline water in QA was established. Then, for each scenario, these maps were combined to obtain the probability of the presence of saline water in the QA at a specific location (x, y). Hence, by subtracting the probability map of the present day against the ones of the scenarios, one can obtain the probability of observing an advance of the interface (Fig. 11).
The maps show a significant advance of the interface, the increase of probability of occurrence can reach more than 50% in some areas. The two circular regions affected in the southeast (approximately 700,000/6,170,000) where the probability becomes close to 1, are probably artifacts due to pumping wells. The permeability in these areas is rather low for the QA (10 −5 m/s, Fig. 4A, B) which can explain this exaggerated drawdown.

Flow regime variations
Due to the modifications in heads in the aquifers caused by the climatic (and anthropogenic) stresses, it is interesting to investigate the evolution of the exchanges between QA and PA. Figure 12 compares the vertical fluxes simulated for today's conditions and for Scenario 3, since the two other scenarios present similar variations but with lower intensity. Under current conditions, the aquifer fluxes are mostly descending, except near the Têt river and the coast where there is a strong ascending flux and it naturally protects the PA from seawater intrusions. Under Scenario 3, the ascending flux tends to decrease, and, on the contrary,  Figure 13 shows the probability of observing a descending flux today. The probability increases successively for each scenario, meaning that each factor (climatic and anthropogenic) has significant impact on the interface. The most problematic areas are those close to the sea and where the flow becomes descending (yellow-red in Fig 13) and enhances the risk of vertical seawater intrusion (VSI). We can identify, for example, one zone in the north which corresponds Fig. 11 Probability of observing an advance of the interface for all scenarios compared to today. These plots show how the interface may advance inland in QA and with which probability. They are obtained by subtracting the probability maps of observing the interface for today and a specific scenario. The gray line represents the median position of the interface for today (i.e. having a 50% probability of observing the interface). The black line represents the median for Scenario 3

Discussion
Parameter identification The first thing to be reminded is the objective of this study and the modeling, namely to analyze and estimate the sensitivity of a coastal regional aquifer to climate changes. For this purpose, we had to proceed with parameter identification, even if this process is often difficult (Anderson et al. 2015;Doherty 2015;Hill and Tiedeman 2006) and requires several compromises to be made, notably in the conceptual model, the choice of parameters, the computational costs, etc. In this case, we mainly oriented the calibration to the hydraulic conductivity distribution and to a lesser extent to the parameters controlling the boundary conditions. The parameter identification is judged satisfactory in terms of heads reproduction with an RMSE ranging from 2.3 to 2.9 m in the PA. The simulated heads in the southeast region (Fig. 7A) are larger than the values reported in the piezometric map. This area is highly uncertain, since a lateral inflow is suspected (HYDRIAD 2014a), but the data are not sufficient to accurately estimate it. It should also be noted that the piezometric map (for the PA) in this area is mainly determined from an interpolation which may not properly represent the real water table, due to an insufficient number of observation points. Nevertheless, our simulations are consistent with the piezometric map of the QA, which also shows a high head at this location. An external flow is thus not necessary, but the calibration needs to lower the K values to increase the heads. The simulated heads in the PA are often slightly too high (mean residual error = + 0.5 m). This issue can probably be attributed to the use of an outdated piezometric map (1990) as additional observation points for the inversion. As the heads drop since decades in the aquifer, they were slightly higher in 1990 than today, especially in the PA. Despite the fact that we used a weight factor of 0.5 for each point, it is probable that it slightly biased the inversion.
In addition, some parameters remain highly uncertain, such as the conductance of the rivers or the irrigation recharge zones. For rivers, this can be easily explained with our choice to impose only one conductance value per river. It is probably not sufficient to capture the complex hydrology of these streams and their interactions with the aquifer. The river bed K is known to have a wide range of values (Calver 2001). Further studies and in situ measurements would be necessary to better constrain these exchanges. However, for a regional purpose, assuming homogeneous values was considered sufficient.
Model limitations and improvements To study the impact of climate change on the Roussillon aquifer in a reasonable amount of time, a few simplifications in the model design were needed, such as the use of a steady-state model. An important simplification is that even with zero-recharge conditions, the rivers cannot become dry since we impose a fixed level as BC. The MODFLOW RIV package that we used accounts for the fact that the flow can be limited if the groundwater levels go below the bottom of the river (instead of continuously increasing), but it does not offer the possibility of drying the river and stopping the recharge completely. Therefore, a possible improvement could be to use the more advanced package SFR (Prudic et al. 2004;Leaf et al. 2021) that allows the implementation of the rivers in a more realistic manner. However, the implementation of this package would require additional data and parameters, not necessarily available, that would need to be given to the model. The way we implemented the rivers tends to mitigate the impact of the reduction of the recharge in the climate scenarios and is, therefore, conservative. We are thus confident that we do not overestimate the impact. When the models indicate a significant impact, it is likely that, in reality, the impacts will be even greater for a given climate and extraction scenario. Finally, the model does not simulate the transient and density transport processes. Therefore, it cannot fully capture the complex circulation of salted water in the different layers of this coastal aquifer. However, the simulation of these processes would require more parameters and their identification would involve a large computational cost, Fig. 13 Probability of observing an inversion of the vertical flux between QA and PA according to Scenario 3 (with respect to today). The yellow-red area is for regions where the flux is estimated to be ascending today but will have a certain probability to become descending. The blue regions represent the exact inverse. The solid gray line displays the coastline which would make the uncertainty analysis much more difficult or even intractable.
Sensitivity to climate change and sustainability In terms of head reduction, the aquifer seems to be highly sensitive to future climate scenarios. This should be retained as the observed decrease in both aquifers (Fig. 9) is quite high where some regions suffer a water table drop of up to 10 m, which will have a significant impact on the water use for these areas. Such head reductions are plausible because they are of the same order of magnitude as the largest changes observed on the piezometric time series (10 m) in the past. Although they may seem significant, the projected future decreases can be considered conservative. Indeed, even Scenario 1, with the smaller impact of climate changes, shows an important decrease of the heads in a large part of the aquifer. Furthermore, the model assumes that the water levels in the rivers are not affected by climate change. This assumption is conservative and has a strong impact on the results; moreover, if we consider that we are using a RCP 8.5 scenario. The rivers are likely to become dry under the climatic conditions of this kind of scenarios, sporadically or permanently. But, as rivers cannot become dry (in the model), they continuously contribute to recharge the aquifer despite the fact that the water table could be several meters below the river. Normally, we would expect a decrease of the river stages and thus a reduction of the recharge from the rivers.
As already mentioned, projected head changes will affect the hydrodynamics of the system and the vertical exchanges between the aquifers (Fig. 12). Compared to the decrease in heads, these changes seem less significant, but they have important consequences. Indeed, it is important to remember that the projections are made in steady state and that, therefore, the simulated flows are permanent. If the model predicts a downward flow in a steady-state regime, this implies that this descending flow is likely to be "permanently" present. This is particularly important for flows near the coast where a natural upward flow (blue, Fig. 12) is present and has a protective effect on PA, preventing vertical seawater intrusion (VSI) from QA. In Scenario 3, this ascending flow would significantly decrease and would even have some probabilities to reverse (Fig. 13), in particular in the north (700,000/6,180,000-6,190,000) and in the south (700,000/6,165,000). Inversions of the flux in the northern part of the aquifer are already reported by some studies (HYDRIAD 2014a), but these are generally restricted to the summer period, under high tourism activity. In these cases, the modeled semi-permeable layer between QA and PA, prevents any intrusions to occur during the period. However, if head changes occur permanently, the quality of PA freshwater could be significantly affected by downward seawater flows.
Furthermore, each projected climate scenario would result in an inland progression of the salt/freshwater interface estimated with the Ghyben-Herzberg relation (Figs. 10 and 11). Near the coast, the reduction of the recharge will mostly reduce the heads in the PA; the sealevel rise will translate upwardly the whole system. Finally, the pumping increase (Scenario 3) will also contribute to reducing the heads in the QA (and the increase of the interface). Thus, even with a conservative scenario (like the 1st or the 2nd), the risk of intrusion is significantly increased, and combined with the risk of flow inversion, VSI from QA to PA can occur more easily. As the PA is mostly used for drinking water, especially near the coast, those scenarios suggest that it will probably be more difficult to deal with future demands related to drinking water, agriculture, or tourism.
Finally, this study highlights that the Roussillon aquifer is strongly controlled by human activities, as reflected in the differences between scenarios 3 and 1. It shows an increase in the risk (and various impacts) that seems to be of the same order of magnitude as the one due to recharge reduction and/or sea-level rise. However, locally, the climatic impacts affect the aquifer differently. In fact, the rise in sea level greatly affects the aquifer near the coast, while the withdrawals have a rather homogeneous impact over the plain because we assumed a uniform increase of 20%. It is difficult to clearly distinguish the effect of each factor as they are coupled. More scenarios could be investigated (only increase in pumping, only rise in sea level, etc.) to distinguish their respective impacts. However, scenario 3 still indicates strong human control over the water resource under changing climatic conditions. The anthropogenic footprint can also be seen in the total budget (Fig. 8) and is consistent with previous studies (Caballero and Ladouche 2015). The withdrawals represent nearly 40% of the total flows and irrigation from irrigation channels has a contribution that can locally exceed the natural recharge by 2-3 times. This human control exists for centuries (Caucanas 1988) and is an integral part of the hydrodynamic system today. In order to maintain sustainable water resources and local economic activities, a strong monitoring of water levels and salinity will be necessary, as well as a better redistribution of the water resource on the plain (e.g. by bringing more water from the land to the coast). Investigating how to optimize this redistribution and the entire water supply system can be done by coupling groundwater models (such as the one presented in this document) with hydro-economic models (Girard et al. 2015).

Conclusion
In this study, the response of the Roussillon aquifer to climate change was investigated by building and calibrating a new numerical groundwater model and forcing it with simplified future stress scenarios. Three climatic and anthropogenic scenarios were defined on the basis of the available future precipitation projections and the observed water demand trends. The proposed workflow has proved to be efficient to quantify the sensitivity of the system at a regional scale and the corresponding uncertainty. The main novel results are the following.
• The Roussillon aquifer is likely to suffer a median head reduction ranging between 0.5 m and several meters (up to 5-10 m) by 2080, particularly in inland areas. • The withdrawals are an important driving force of the hydrodynamic system. A moderate increase of them could significantly affect the water resource and increase the risk of seawater intrusion (by 10-20%) through downward fluxes from the superficial Quaternary aquifer to the deep confined Pliocene aquifer. • Under projected climatic conditions, the ensemble of models projects a weakening, or even locally an inversion, of the ascending flux from the Pliocene aquifer to the Quaternary aquifer along the coast. This led to a higher risk of downward vertical seawater intrusion in the Pliocene aquifer, the main source of drinking water in the region.
Due to its regional importance for irrigation and drinking water production and its sensitivity to climate change, this work shows that the Roussillon aquifer should be monitored and studied with attention in the future. Several aspects of its behaviors must still be investigated to refine the forecasts presented in this paper. Finally, despite the extensive anthropic use of this aquifer for centuries, it is more than likely that humans will have to deal even more with a very important actor in the future, climate.