The Impact of Predicted Climate Change on Groundwater Resources in a Mediterranean Archipelago: A Modelling Study of the Maltese Islands

: The effects of changes in climate predicted for 2100—reduction in recharge, increase in water demand and sea-level rise—on groundwater volume and saltwater intrusion have been quantiﬁed in the Maltese Islands, an archipelago located at the center of the Mediterranean Sea. A three-dimensional density dependent and heterogeneous model, working in transient conditions, was developed based on morphological and geological information. The hydraulic conductivity and porosity of the lithological formations were derived from previous tests and studies conducted on the islands. The complex fault system intersecting the area has also been included in the model. The results show that among the three considered factors affecting groundwater resources, the most signiﬁcant is the increase in water demand, which is closely followed by the decrease in groundwater recharge. Sea-level rise plays a marginal role. The 80-year simulation period showed that these combined impacts would cause a loss of more than 16% of groundwater volume. The obtained present-day were used in the calibrated as initial for the of The simulation results show that three is the in the in groundwater Sea-level rise plays a marginal role. The 80-year simulation that these a more than Based on the modelling results, we suggest


Introduction
The predicted effects of climate change will worsen current environmental problems globally such as land use issues, pollution, water scarcity and biodiversity decline [1]. The most vulnerable societies are characterized by insufficient observation systems and monitoring equipment, and they also lack impact-based models. This is particularly the case of the Mediterranean region. Here, average annual temperatures are now approximately 1.5 • C higher than during the period of 1880-1899, which is well above current global warming trends [1]. It is expected that the future average annual temperature increase will be 2.2 • C in 2040, possibly exceeding 3.8 • C in 2100 in some Mediterranean regions [2,3]. In addition, climate models suggest scenarios with reduced rainfall in the coming decades [4]. The combination of warming and reduced rainfall will result in drier conditions. The frequency and intensity of droughts have already increased significantly in the Mediterranean region [5]. Scenarios with temperature increases of 2-4 • C up to the 2080s for Southern Europe would imply more widespread decreases in precipitation of up to 30% in comparison to present levels [6]. As a result of all of these factors, water availability in the Mediterranean Basin is likely to undergo the largest decrease in the world by 2-15% for 2 • C of warming [7,8]. Water scarcity is expected to be exacerbated by two additional factors. First, populations are expected to grow significantly, especially in countries that are already facing water stress. Irrigation, for example, represents 50-90% of the total Mediterranean water demand [9]. This demand, which is projected to increase by 4-18% by the end of the century due to climate change alone, may increase by 22-74% if population growth is taken into consideration [10]. The second factor is sea-level rise due to global warming, which in the last two decades has increased by about 3 mm per year [11]. Depending on the method used, the global mean sea-level rise in 2100 will range from 52 to 190 cm [12,13].
Improved understanding on how climate change affects the hydrological system, including groundwater, is necessary for developing suitable adaptation strategies [14]. Numerical models have played a key role in assessing how climate change is expected to influence groundwater resources. Groundwater over-extraction is a widely known cause of saltwater intrusion, which has been extensively studied under different pumping regimes and conditions (e.g., [15]). Aquifer salinization can result in the quality loss of precious water resources to the point that their use for drinking or irrigation purpose is no longer possible. In the first case, salt contamination exceeds the standards set to preserve people's health, in the second case, soils fertility can be strongly compromised. The impact of climate change on saltwater intrusion has been addressed by a number of studies in the past decade. Site specific experimental studies (e.g., [14,16,17]) and synthetical simulated cases (e.g., [18][19][20]) have been performed for different confined and unconfined aquifers using more or less complex models, such as density-dependent flow and solute transport models, numerical models based on the Ghyben-Herzberg assumption and even analytical solutions. The results are extremely variable since they are strongly influenced by sitespecific settings. Hydrogeological properties, climate forcing and boundary conditions contribute to a wide range of possible scenarios. Werner et al. [18] together with Morgan and Werner [21], for example, found that aquifers in which groundwater discharge flux to the sea is set to be constant are less sensitive to sea-level rise with respect to constanthead governed systems. The impact of other factors such as changes in precipitation and temperature, which influence evapotranspiration and groundwater recharge, has also been assessed, but it is not obvious, within a specific context, whether they will have a cumulative effect or will offset each other. Green and MacQuarrie [22] investigated the effect of projected sea-level rise and climate change on recharge and of groundwater extraction rates on seawater intrusion. The setting in this case was an unconfined sandstone aquifer located on the coast of the Gulf of St. Lawrence (Canada). Although the relative importance of the three factors changes according to the specific location, they found that sea-level rise had the least impact on seawater intrusion. The effect of declining recharge was most significant at shallow to intermediate depths along the transition zone, while the impact of increased pumping rates was limited to the area close to the well and at the same depth of extraction.
Providing access to a potable water supply is a great challenge for small island communities where groundwater resources are constantly being threatened by saltwater intrusion. However, there is a paucity of modelling studies of the effects of changing climate on island water resources (e.g., [16,21]), particularly in the Mediterranean Sea. The Maltese archipelago has one of the highest population densities in the world and is one of the poorest regions globally in terms of water resources per inhabitant [23]. Potable water supply is mainly dependent on groundwater abstraction from limestone aquifers, and seawater upconing into the pumping wells is a serious problem [24]. In this study we develop a three-dimensional density dependent and heterogeneous model, working in transient conditions, in order to quantify the potential impact of predicted climate change (sea-level rise, groundwater recharge decrease and water demand increase) up to 2100 on the freshwater resources of the Maltese archipelago.

Geology
The Maltese archipelago comprises three main islands located in the Strait of Sicily, around 95 km south of Sicily. The Maltese Islands comprise one of the few emergent areas of the Pelagian carbonate platform. Outcropping across the Maltese Islands includes five shallow water, Oligo-Miocene pre-rift to syn-rift sedimentary successions [25]. These in-clude the following: (i) the Lower Coralline Limestone (LCL)-a 1000 m thick succession of algal foraminiferal limestone; (ii) Globigerina Limestone (GL)-fine-grained biomicrites up to 207 m thick, which are further subdivided into Upper (UGL), Middle (MGL) and Lower (LGL) Globigerina Limestone; (iii) Blue Clay (BC)-up to 75 m of slightly consolidated hemipelagic clays and marlstones; (iv) Greensand (Gs)-poorly cemented bioclastic and glauconitic limestones that are up to 11 m thick; (v) Upper Coralline Limestone (UCL)-an up to 100 m thick shallow water reef complex ( Figure 1). The entire sequence is disrupted by two normal faults systems ( Figure 1) [26]. The first is a system of syn-depositional ENE-WSW trending faults, which were active intermittently between the Early Miocene and mid-Pliocene. The second system consists of NW-SE trending faults that were formed or re-activated in the late Pliocene to Quaternary [27][28][29][30][31] (Figure 1).

Hydrogeology
The carbonate lithologies described in the previous section host three types of meteorically recharged groundwater bodies. The largest is the mean sea-level groundwater body hosted in the LCL and GL formations, which appears as a freshwater Ghyben-Herzberg lens floating above seawater [33]. This groundwater body provides the primary source of potable water and has a residence time of 15-40 years [34]. The aquifer hosting this groundwater body is in contact with the sea along the entire coastline and extends well below sea-level. The other two groundwater bodies have smaller extensions, and they are restricted to specific zones of the Malta island hosted within UCL and are either perched over the impermeable BC formation above sea-level or in valleys below sea-level. Hydraulic conductivity and porosity of these formations are known from previous tests and studies conducted on the islands (Table 1), whereas the specific storage was derived from the literature [35]. There are no significant surface water bodies across the Maltese Islands; thus, recharge of these three groundwater bodies is largely dependent on a mean annual precipitation of 553.1 mm recorded for the period of 1951-2010 [44]. The Maltese Islands have a semi-arid Mediterranean climate characterized by hot, dry summers and mild, wet winters. The most significant rain events usually occur between October and February, whereas the main tourist season coincides with the driest months [34]. This combination of events causes periods of unsustainable exploitation with aquifers that are not fed by any superficial recharge but undergo significant withdrawals for the sudden increase in population. The bottom of the main (mean sea-level) groundwater body, which is the subject of this study, is non-constant in time and defined by a mixing zone between fresh and saltwater. The extent and thickness of the mixing zone is controlled by the fluctuation of fresh groundwater level caused by abstractions, recharge or even atmospheric pressure variation [33]. Groundwater flow regimes are also influenced by the dense system of faults intersecting the archipelago since they represent an element of heterogeneity and discontinuity in the geological formations (Figure 1), whereas land use, abstraction for public supply and irrigation and recharge from distribution pipe leakage influence the groundwater budget ( Table 2). These data were obtained by putting together the results of the groundwater balances performed in the documents mentioned in the title of Table 2, part of which are peerreviewed articles [34] or official records of Maltese governance [45,46].
The main aquifer hydraulic head data were collected in 1944 and 1990 by the Bureau de Recherche Geologique et Miniere [47] which used two different sets of monitoring wells, which consist of a total of 52 observation points. In 2014, hydraulic head data were recorded by the Malta Resources Authority (MRA) who monitored a third set of 17 observation points. In 1944, measured groundwater levels ranged between a maximum of around 3.5 m and a minimum slightly below 1 m above sea level near the coast. Between 1990 and 2014, a general decrease in aquifer levels can be observed with time (the maximum level registered a drawdown of 1 m).

Modeling
The universal groundwater model MODFLOW-2005 [48] and the Seawater Intrusion package SWI2 [49] were used to simulate variable density flow in the mean sea-level aquifer. The SWI2 Package is designed to simulate regional seawater intrusion in coastal aquifer systems. Variable-density flow is represented by discrete zones of uniform density. The surfaces bounding the zones are density isosurfaces. Water density does not vary in the freshwater and saltwater zones. The governing equations for vertically integrated variabledensity groundwater flow adopted by the SWI2 package have been extensively described by Bakker et al. [49]; thus, only a brief description on how the package works will be provided here. The numerical approach disregards diffusion and dispersion phenomena, so the resulting differential equations are equivalent in form to those implemented by the MODFLOW code for uniform-density flow. The density effects are incorporated into MODFLOW by the inclusion of pseudo-source terms relative to the groundwater flow equation, and there is no need to solve a separate advective-dispersive transport equation. Vertical and horizontal movements of defined density surfaces are calculated separately by using a combination of fluxes. The latter is estimated by solving the groundwater flow equation and by using a simple tip and toe tracking algorithm. The main advantage provided by this approach is that fewer model cells are required, since aquifer domains can be represented by a single layer. All these factors result in substantial model run-time savings, which is an advantage when dealing with a regional aquifer.  [32] with legacy [50] and recently acquired borehole logs (from a 2017 geotechnical study to assess the feasibility tunnel construction between Malta and Gozo). The entire domain was horizontally discretized by a mesh of 100 m × 100 m square elements rotated by 53 • . The vertical extent of the model is covered by four layers with different thicknesses, which are designed to replicate the heterogeneity created by the geological formations across the Maltese Islands. The model layers are set to be "convertible," which means that the model checks the hydraulic head distribution to determine if it is above the cells of a specific layer. In this case, the layer is considered by the model as confined or under pressure, otherwise the model treats it as phreatic. The bottom layer has a thickness of 500 m and represents the predominant geological formation of LCL, whereas the remaining three overlying layers have a total thickness of only 5 m in which five more formations/members are condensed together with the previous one (LGL, MGL, UGL, BC and UCL). A much denser vertical discretization would have been necessary to couple groundwater flow and dispersive solute transport codes. The horizontal hydraulic conductivity and porosity, which was set to be equal to the specific yield, together with the storage coefficient of the geological formations are derived from Table 1 (Figure 2). Vertical hydraulic conductivity has been set to be a tenth of the horizontal one, as has been applied in many modelling applications, in order to account for anisotropy (e.g., [16]).

Model Domain, Properties and Boundary Conditions
A further element of heterogeneity that has been included into the model comprises faults that intersect the area (Figure 3). The faults are considered vertical because they are normal faults dipping at 65-70 • , which is difficult to model in consideration of the grid cell dimensions. Faults play a key role in the model calibration process. The representation of their conformation and dimension is directly linked to that of the mesh cells. In order to ensure a more realistic performance of their interaction with fresh and saltwater flow, their hydrodynamic properties have been specifically estimated. This aspect is discussed in more detail in the following sections.  The hydrological and hydro-geological forcings derived from Table 2 were assigned to the transient model according to the relevant time period. Natural aquifer recharge from rainfall has been homogeneously distributed on the island's domain, and urban leakage recharge has been added to this component. The annual volume of groundwater withdrawn for irrigation has been deducted from the recharge volume, since the exact numbers and positions of the involved wells are not always precisely known. The distribution of public supply facilities is known, however, so that all the wells and the related flow rates have been included and modelled. The pumping points are placed at a depth ranging between 0 and −10 m with respect the mean sea level. The observation wells monitored by the BRGM and the MRA (see Section 2.2) were also precisely placed within the model and taken into consideration during the calibration procedure. A constant head boundary condition of 0 m pertaining to the sea has been assigned to the cells. An initial groundwater hydraulic head, corresponding to the top elevation of the model, has been set at the inland cells. An associated initial depth for the freshwater/saltwater interface ("ζ surface" from now on) has been defined by using the Ghyben-Herzberg analytical solution.

Model Calibration
A long-term transient simulation (ca. 1000 years) was run to reach present-day conditions by utilizing a complex calibration procedure. The simulation time has been divided into five different chained intervals belonging to the same model provided with the aforementioned initial conditions. Several centuries were run in this first temporal window to reach steady state conditions in terms of both groundwater heads and ζ surface position. The recharge value recorded in 1944 has been used for this purpose. The model status reached after this first long-time interval, in turn, provides the initial condition for the subsequent linked timestep, and this process recurs sequentially and automatically for the other time intervals. They cover the periods 1944-1977, 1977-1999, 1999-2004 and 2004-2016 and involve all the respective components presented in Table 2 in the calculations. A time step of 365.25 days (one year) has been employed for the entire simulation.
The model thus set has been run into an iterative optimization procedure for the estimation of the hydraulic conductivity of the faults. The automatic calibration was driven by UCODE-2005 [51]. This parameter has been managed as follows: horizontal and vertical hydraulic conductivity values are assigned to faults equivalent to that of the geological formations in which they occur. Two different constants, set for the entire fault system and possessing an initial unitary value, were respectively multiplied by the aforementioned horizontal and vertical conductivities. The hydraulic head values, which were estimated from the model observation points and representing the monitoring wells, were compared with the measured data. The procedure changed the values of the two constants until the deviation between calculated and observed values reaches a minimum value [52]. The hydraulic conductivities of the faults are ultimately a multiple or a fraction of the conductivity of the lithology hosting them. This procedure was carried out using hydraulic heads measured in 1944, 1990 and 2014 and by enforcing the model to reproduce the observed groundwater level distributions precisely in those years and as closely as possible.

Climate Change Simulation Strategy
Once the calibration process was successfully concluded, the simulated present-day conditions have been used as initial conditions for the climate change scenarios. Different scenarios, with increasing severity conditions, were considered within the variation ranges estimated for sea-level rise, decrease in precipitation and increased water demand expected by 2100 and described in the Introduction. Specifically, a sea-level rise of 0.5, 1.2 and 1.9 m has been considered, together with a decrease in groundwater recharge of 10, 20 and 30% and an increase in water demand of 22% (4% for irrigation abstraction + 18% for public supply) and 74% (18% for irrigation abstraction + 56% for public supply). An overall set of 35 transient simulations, covering a period of 80 years (2020 to 2100), has been carried out by combining these three effects (Table 3). Climate change impacts have then been quantified in terms of ζ surface displacements and related freshwater volumes losses.

Model Calibration
The first modeling stage, represented by the inversion of the hydraulic head data, was performed in order to estimate two different multiplicative constants, respectively, associated with horizontal (K H ) and vertical (K V ) hydraulic conductivities of the lithologies affected by faults. Figure 4 shows the outcomes of the model at the end of the iterative optimization procedure. Values of 0.091 × K H and 50.12 × K V were obtained. The goodness of the fit was determined at the end of the inversion procedure via the Nash-Sutcliffe index (NSI) for all observed scenarios. This index establishes the relative magnitude of residual variance (noise) compared to measured data variance (information) [53]. An NSI value of 0.914 was calculated for the 1944 dataset: 0.829 for the 1990 observations and 0.871 for the 2014 data. It can be observed that the location and the number of the monitoring wells in the maps changes in the three measurements campaigns. This aspect was specifically addressed in the modeling phase and considered during the calibration process. Figure 5 shows the map of the ζ surface at the end of the optimization procedure depicting the spatial distribution of the bottom of the aquifer in 2014.

Climate Change Impacts
The second simulation phase started from the calibrated model (present day conditions) to predict the effects of different climate change scenarios (sea-level rise, recharge decrease and water demand increase) on the groundwater resources of the Maltese Islands. Initially, the effects of each forcing, undergoing the range of variations expected by 2100, were investigated separately. The results shown in Figure 6 are the average changes in the position of both the water table of the aquifer and the freshwater/saltwater interface (ζ surface) calculated for 2100, relative to the levels in the same year, without any climate change. It can be observed that sea-level rise, adopted in the simulations, generates a lifting up of the ζ surface of a few centimeters, while the water table gains an average hydraulic head increase of almost 2 m. On the other hand, a completely different behavior of the system is obtained when either the effects of the forecasted reduction in water recharge, caused by the increasingly scarce rainfall, or when the increase in groundwater exploitation, due to a growing water demand, are considered. In these two cases, the water table is subject to a general drawdown with respect to the reference condition, and this is highlighted by the negative hydraulic head variations ( Figure 6). The ζ surface, in return, moves up several meters causing significant seawater intrusion. The two different perturbation rates (recharge decrease and water demand increase), caused by the climate change, appear to have a similar impact with slightly more pronounced effects in the second case. Figure 7, obtained by means of ordinary Kriging performed on the results of the simulated scenarios (from 18 to 35) described in Table 3, shows the extent of groundwater body reduction in terms of percentage of groundwater volume loss evaluated by the model in 2100. This volume loss was determined, once more, by taking as a reference the position of the saltwater-freshwater interface (ζ surface) calculated in 2100 with present day conditions and by comparing the latter with that produced by different scenarios generated, this time, by the combination of climate change effects. The effects produced by the combination of the three forcings reflect the results shown in Figure 6. If we consider the impact generated by the sea-level rise, we note that the contour lines, which describe freshwater volume loss (percentage), are almost parallel to the horizontal axis. The slight slope of these contour lines associated with the color change point out, without significant variations, that the higher the sea-level rise, the lower the freshwater loss. This is consistent with the small positive variation provoked by sea-level rise on the ζ surface average position and the average increase in the water table hydraulic heads in Figure 6. The decrease in groundwater recharge has, instead, an important effect since it can result in a loss of freshwater volume of up to 5-6% for all the considered scenarios. Combining this with the component represented by the increase in water demand (22% and 74%) results in a rise of the percentage of lost freshwater volume by an additional~4 and 10% respectively, if compared to the loss recorded without this occurrence. The combined impacts, in the worst scenario, would cause a loss of more than 16% of the volume of the groundwater body. Figure 8 shows vertical cross sections highlighting groundwater body contraction caused by the worst climate-driven changes, derived from the plot in Figure 6 (30% decrease in groundwater recharge, 74% increase in water demand and a sea-level rise of 0.5 m), with respect to the reference condition (no climate change effects) in 2100. Red areas represent groundwater zones included between the two ζ surfaces determined by the conditions described above and provide a picture of the saline intrusion caused by the aforementioned factors. The most severe saltwater intrusion is observed in the points where the section lines cross the extraction wells (actual public supply facilities), which are represented on the map by small squares with a yellow perimeter (Figure 8). It can be observed that, in some places, the saltwater interface, determined by the worst conditions, is significantly close to the water table, clearly affecting the wells located there and, thus, making them unusable. This is clearly shown in the section line A-A (Figure 8). The southeast side of the island of Malta is characterized by a high pumping well density, and this is reflected, on the right side of the vertical representation of the aforementioned section, with a significant rise in the saltwater interface. The central part of the island instead has a lower density of wells and, in comparison, shows a smaller extent of saltwater intrusion.

Influence of Faults on Groundwater Modelling Results
The calibration procedure of the model, which involved the joint evaluation of hydraulic heads in three different periods each of which characterized by distinct sets of observation points, has delineated the role of the complex fault system on groundwater circulation. The geological reconstruction of the archipelago exclusive of faults did not allow the model to reproduce the piezometric trend observed in the three measurement campaigns. The introduction of faults and the estimation of their hydraulic properties, however, allowed the model to generate results that are closer to the monitored behavior. The ratio between the horizontal and vertical hydraulic conductivity of the faults turned out to be about 1:500, whereas the anisotropy in the rest of the domain was set to 1:10. This is the first important outcome of the study, and the adopted approach can also be transferred to other islands where faults are widespread. Lotti et al. [24] developed a 2-D model of the Mean Sea Level Aquifer focusing only on the island of Malta and including the main faults only. Their simulations were performed with the same code adopted in this study, and the calibration process was driven by using data collected in the same time window. The faults in the island were accounted for by applying HFB-Horizontal Flow Barrier package of MODFLOW-while, in our case, faults are represented by means of sudden variation in the porous media hydraulic properties. Fault isotropic hydraulic conductivities estimated in Lotti et al. [24] range, for most of the cases, between an interval comprised by 8 × 10 −8 and 1 × 10 −2 m/d, which is not far from that resulting from this study and presenting horizontal conductivity values ranging from 2 × 10 −6 and 1.7 × 10 −2 m/d. Moreover, the maps of the ζ surface, resulting at the end of calibration processes (present day conditions, time at which their study stops), are comparable. A maximum depth of −160 m asl can be, in fact, observed in both results, and the point-to-point discrepancies encountered in the maps can be attributed to different configurations of the models (2D vs. 3D) and to the higher complexity of the fault system included in this study. The greater complexity of the model developed in this study is justified by its use as a forecasting tool for assessing the impacts of climate change. The more responsive to reality the adopted predictive instrument is, the more accurate the results will be.

Impact of Predicted Climate Change on Groundwater Resources
The volume of freshwater stored in the main aquifer of the Maltese Islands underwent a visible reduction from 1944 to 2014, as suggested by the progressive drawdown in hydraulic head (Figure 4). For example, the area in the south of Gozo presents, in the three different periods for which the model was calibrated, groundwater levels becoming increasingly close to the average sea level. This is due to the fact that this area hosts a greater concentration of extraction wells of the public supply facilities (Figure 8). This should already constitute a wake-up call for the severe saltwater intrusion in that area for which corrective actions should be undertaken.
The impact of climate change on groundwater resources has been evaluated starting from the situation in 2014 projected to present day conditions via forward modelling up to 2020. The increase in water demand and the decrease in the aquifer recharge play a strong role in predicted freshwater volume loss. This aligns with the conclusions of Kløve et al. [54], who concluded that the effects of climate variability on groundwater are likely more pronounced in smaller aquifers. In areas where groundwater systems are likely to be affected by climate change, sustainable groundwater management is challenging, and modeling/forecasting tools, such as the one developed in this study, could be of help to support decisions to improve the long-term sustainability of groundwater extraction [55].
Sea-level rise does not result in significant volume reduction in the groundwater body (Figures 6 and 7); on the contrary, it results in an increase in the average level of the aquifer. The saltwater-freshwater interface is not very sensitive to sea-level rise as its position does not change much after a modeling time of 80 years. As a result, the final aquifer volumetric balance slightly increases with an increase in the sea-level. Although this result suggests a positive effect of sea-level rise, one needs to consider that the investigated simulation time does not allow the interface to reach steady state conditions. If we were to model the same scenarios with much longer simulation times, the ζ surface position would likely reach higher elevations. Furthermore, the Maltese Islands are characterized by gently to steeply sloping coasts, and this limits the extent by which the aforementioned phenomenon reduces the onshore terrain. The hinterland portion, gained by the sea with the highest considered level rise, is well below the cell extension of the model mesh (a coast belt with a thickness of about 2-3 m compared to 100 m mesh elements). This aspect was, therefore, neglected; thus, the extension of the modeling domain was kept constant during the simulation of the various scenarios. In other circumstances, a significant loss of land to sea-level rise would have been included in the model calculations, which would have made the study more complex. Finally, the rise of the water table can generate unpleasant inconveniences such as the flooding of underground structures and geographically depressed areas and contact with underground infrastructures such as sewers (resulting in Escherichia Coli contamination problems) and power lines. Our results are consistent with the findings of Green and MacQuarrie [22], who performed a similar study on coastal sandstone aquifers, and Lemieux et al. [16], who studied an island system characterized by terrigenous sediments, carbonate and evaporite formations. They found that after performing transient saltwater intrusion modeling over 89 years and 83 years, respectively, sea-level rise had the least significant impact on the position of the freshwater-saltwater interface. A decrease in groundwater recharge was more important. Numerical simulations, performed on both confined and unconfined systems in which sea-level rise has no significant impact on saltwater intrusion, are also described by Chang et al. [56]. Moreover, the studies by Ferguson and Gleeson [20] and Hansen [57] on coastal aquifer vulnerability concluded that poor operating practices for drinking water supply are more important than sea-level changes with respect to saltwater intrusion intensification. This last point is particularly evident in the cross sections of Figure 8 where the interface, obtained with the worst-case conditions, reaches the highest positions in correspondence with the extraction wells. The interface position, under significant pumping, can approach the well screen where turbulent phenomena, caused by the high extraction rate, can produce interface dispersion with consequent well encroachment and aquifer contamination [58]. This could be a problem for the Maltese Islands since the groundwater body is thin in some places.
The simulation results in Figure 8 are based on a pessimistic scenario, as they are characterized by the combination of the worst conditions foreseen for the Mediterranean region. However, they elucidate certain aspects that should be taken into consideration in future planning of the available water resources management. Our calculations suggest that the anthropogenic impact, generated by the exploitation of the aquifer, is the most significant factor influencing the reduction in freshwater volume and its contamination by saltwater intrusion. This means that actions can be taken to mitigate this impact. The wells used for public water supply are concentrated in South-East Malta, and the most severe saltwater intrusion occurs here (Figure 8). The layout of the wells has been dictated by the local population distribution and the shallow location of groundwater bodies. The simulations show how the increase in water demand associated with current pumping well distributions causes a significant decrease in both the quantity and quality of the freshwater resource in South-East Malta. This phenomenon is much less pronounced in areas with a lower well density. All this suggests that a more homogeneous distribution of the pumping points across a larger surface area of the island could improve the resilience of the groundwater body against increased exploitation in the future.

Limitations of Groundwater Modelling Results
The predictions presented here are associated with some limitations due to a lack of available data and specific studies on the aspects that will be discussed below. Groundwater recharge fluctuations, due to the different seasonal regimes, or tidal effects, which can also affect saline water intrusion [59], have not been considered in the simulations. Furthermore, the methodology used to represent climate change impacts is based on average data derived for the entire Mediterranean zone and not for on-site specific observations. Including all these aspects in the modeling process is difficult in terms of data collection and their integration into the governing equations, but it would result in more realistic and sitespecific outcomes. The number and position of the irrigation wells in the Maltese Islands are not always available in the published literature. Thus, they were not included in the model, and their overall withdrawn annual volume has been deducted from the recharge one. This implies that the localized effect of each well, in the restricted portion of the aquifer in which it is located, is instead incorporated into a wider scale. This means that the overall impact generated by the increasing water demand required by the agricultural sector is included in the water balance of the model, but the local effects generated by the actual presence of wells are missing.

Conclusions
We assess the impact of predicted climate change on the freshwater resources of the Maltese islands by 2100. The study quantifies the effects of reduction in recharge, increase in water demand and sea-level on groundwater volume and saltwater intrusion. We have developed a three-dimensional density dependent and heterogeneous model, working in transient conditions, by using morphological and geological information. The model also incorporates the complex fault system intersecting the area, which plays a key role in the groundwater flow regime. The hydraulic conductivities of the faults have been estimated by running the model into an iterative optimization procedure. The adopted methodology has allowed us to reproduce the aquifer hydraulic head data monitored in three different time periods (1944, 1990 and 2014) and to define the manner in which the faults affect the balance between fresh and saltwater. The results show that when groundwater runs into the fault system, the flow in the horizontal direction is slowed down, since faults horizontal hydraulic conductivities were found to be a fraction of those of the geological formations in which they occur, while groundwater flow in the vertical direction is instead enhanced since vertical hydraulic conductivities result, instead, in a multiple. This is consistent with the goodness of the fit obtained for all the observed scenarios involved in the inversion procedure. The implemented approach is of broad applicability and can be transferred to other sites characterized by a non-negligible presence of faults. The obtained present-day conditions were then used in the calibrated model as initial conditions for the simulation of climate change scenarios. The simulation results show that among the three considered impacts, the most significant is the increase in water demand closely followed by the decrease in groundwater recharge. Sea-level rise plays a marginal role. The 80-year simulation period showed that these combined impacts would cause a loss of more than 16% of the volume of the groundwater body. Based on the modelling results, we suggest that a more homogeneous distribution of the pumping points serving the public supply facilities involving a larger surface area of the islands, together with a more organized use of the irrigation wells, could improve the resilience of the groundwater body against climate change in the future.   677898 (MARCAN)). The project has also received financial support from the University of Malta within the framework of the SMART project (Helmholtz European Partnering Initiative; PIE-0004).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Data are all already present in the paper. We choose to exclude this statement since no external data are available.