Groundwater Flow Model Calibration Using Variable Density Modeling for Coastal Aquifer Management

: The paper investigates the mechanism of seawater intrusion and the performance of free and open-source codes for the simulation of variable density flow problems in coastal aquifers. For this purpose, the research focused on the Marathon Watershed, located in the northeastern tip of Attica, Greece. For the simulation of the groundwater system, MODFLOW, MT3DMS and SEAWAT codes were implemented, while sensitivity analysis and calibration processes were carried out with UCODE. Hydraulic head calibration was performed on the MODFLOW model


Introduction
The salinization of coastal aquifers due to seawater intrusion has been investigated for many years.The diversity and complexity of groundwater formations form different mechanisms and patterns of seawater intrusion, complicating the management and mitigation of the phenomenon [1].The mechanism of seawater intrusion (SWI) and the saltwater wedge has been investigated through several methods.The most common methods that are examined in the literature are geophysical investigation, e.g., refs.[2][3][4], laboratory experiments, e.g., refs.[5,6], mathematical models, e.g., refs.[7,8] or a combination of the above.
Different models and numerical solutions are utilized for the simulation of seawater intrusion.Dokou and Karatzas [9] used the finite element FEFLOW code [10] to model SWI in a karstic aquifer, while Stoeckl et al. [11] used the same code for investigating the effects of post-pumping in the SWI mechanism.Na et al. [12] employed the TOUGH2 multi-phase transport simulator [13] to model an experimental setup and test the influence of seawater density in a coastal aquifer under confined conditions.The SEAWAT [14] threedimensional finite difference code by the United States Geological Survey (USGS) is used by many researchers for the investigation of seawater intrusion dynamics, e.g., refs.[15][16][17].The simplest approach of the sharp interface with the SWI2 package of MODFLOW [18] is also utilized, while other researchers have applied the MODFLOW and MT3DMS [19] codes to simulate seawater encroachment on the mainland [20].
The sensitivity analysis and calibration process are very important when either simulating a natural system or in cases of experimental setups where observation data are obtained, e.g., refs.[21,22].In the field of seawater intrusion, the calibration process is much more complicated since the modeler has to deal with several issues, such as the use of concentration data, tidal data or head data that are affected by variable density [23].Due to the complexity and the high computational requirements of variable density flow codes along with the uncertainty of head measurements in aquifers affected by seawater intrusion [23], many researchers simplify this procedure by utilizing constant density models, e.g., refs.[24][25][26].
For the purposes of this research, the mechanism of seawater intrusion was investigated through numerical modeling in the case of the Marathon coastal aquifer in Greece.
The MODFLOW [18], MT3DMS [19] and SEAWAT codes [14] of the United States Geological Survey (USGS) were employed for this process.The implementation of the two groundwater flow codes (MODFLOW and SEAWAT) served two different purposes: (i) to explore seawater intrusion mechanism in the aquifer and (ii) to assess parameter transferability between MODFLOW and SEAWAT, i.e., the calibration of the MODFLOW parameters using observed head values and afterward, the use of the exact same dataset without calibration in the SEAWAT model.The main outcome of the research is the improvement of the groundwater flow model by comparing the performance of the SEAWAT code against MODFLOW in terms of the hydraulic heads.With respect to the comparison of the results of MODFLOW and SEAWAT in terms of the hydraulic heads, it is to be expected that SEAWAT will perform better after calibration as it entails a variable density factor.However, the parameter transferability, as studied herein, constitutes a research question.Furthermore, the SEAWAT variable density flow model revealed the irregular shape of the saltwater wedge at the eastern part of the plain, affected by the geological heterogeneity in the study area.

Materials and Methods
A graphical representation of the methodological framework utilized in the present research is presented in Figure 1.The following sections describe, in detail, each step of the methodology.
Hydrology 2024, 11, x FOR PEER REVIEW 2 The sensitivity analysis and calibration process are very important when either ulating a natural system or in cases of experimental setups where observation dat obtained, e.g., refs.[21,22].In the field of seawater intrusion, the calibration proc much more complicated since the modeler has to deal with several issues, such as th of concentration data, tidal data or head data that are affected by variable density Due to the complexity and the high computational requirements of variable density codes along with the uncertainty of head measurements in aquifers affected by seaw intrusion [23], many researchers simplify this procedure by utilizing constant de models, e.g., refs.[24][25][26].
For the purposes of this research, the mechanism of seawater intrusion was in gated through numerical modeling in the case of the Marathon coastal aquifer in Gr The MODFLOW [18], MT3DMS [19] and SEAWAT codes [14] of the United States Ge ical Survey (USGS) were employed for this process.The implementation of the groundwater flow codes (MODFLOW and SEAWAT) served two different purposes: explore seawater intrusion mechanism in the aquifer and (ii) to assess parameter tran ability between MODFLOW and SEAWAT, i.e., the calibration of the MODFLOW pa eters using observed head values and afterward, the use of the exact same dataset wi calibration in the SEAWAT model.The main outcome of the research is the improve of the groundwater flow model by comparing the performance of the SEAWAT against MODFLOW in terms of the hydraulic heads.With respect to the comparis the results of MODFLOW and SEAWAT in terms of the hydraulic heads, it is to b pected that SEAWAT will perform better after calibration as it entails a variable de factor.However, the parameter transferability, as studied herein, constitutes a res question.Furthermore, the SEAWAT variable density flow model revealed the irre shape of the saltwater wedge at the eastern part of the plain, affected by the geolo heterogeneity in the study area.

Materials and Methods
A graphical representation of the methodological framework utilized in the pr research is presented in Figure 1.The following sections describe, in detail, each st the methodology.

Study Area
The Marathon Plain is located at the NE of Attica (Greece) and is formed within a hilly marble formation at the north, NE and SW, with the Mediterranean Sea at the south (Figure 2).The plain is approximately 40 km 2 , and it forms a typical coastal Mediterranean catchment.It is an area of great interest due to historical events (the battle of Marathon in 490 BC) and its coastal ecosystem (Schinias National Park), forming a resort available for several activities at the coastal zone of Attica.The climatic condition is semi-arid, with dry and warm summers and mild winters.According to weather data of the National Observatory of Athens meteorological station in the area (http://penteli.meteo.gr/stations/neamakri/, accessed on 22 January 2024), the average annual temperature is 18.3 • C, while a wide range of temperature changes has been observed during the year from −1.8 • C in January to 43.3 • C in August.The average annual precipitation is around 500 mm, occurring mainly between October and March, with the highest precipitation recorded in December.During the dry season, the monthly rainfall does not exceed 30 mm, while the driest month is August, with a range of 0 up to a few millimeters of precipitation.

Study Area
The Marathon Plain is located at the NE of Attica (Greece) and is formed within a hilly marble formation at the north, NE and SW, with the Mediterranean Sea at the south (Figure 2).The plain is approximately 40 km 2 , and it forms a typical coastal Mediterranean catchment.It is an area of great interest due to historical events (the battle of Marathon in 490 BC) and its coastal ecosystem (Schinias National Park), forming a resort available for several activities at the coastal zone of Attica.The climatic condition is semi-arid, with dry and warm summers and mild winters.According to weather data of the National Observatory of Athens meteorological station in the area (http://penteli.meteo.gr/stations/neamakri/,accessed on 22 January 2024), the average annual temperature is 18.3 °C, while a wide range of temperature changes has been observed during the year from −1.8 °C in January to 43.3 °C in August.The average annual precipitation is around 500 mm, occurring mainly between October and March, with the highest precipitation recorded in December.During the dry season, the monthly rainfall does not exceed 30 mm, while the driest month is August, with a range of 0 up to a few millimeters of precipitation.The groundwater system consists of two aquifer units: a karstic aquifer developed in marble and a granular aquifer that overlies the marble formation developed in terrestrial, alluvial deposits (Figure 2).The impermeable bedrock of the karstic aquifer is composed of a schist layer.
Along with the aquifer formations, the Marathon hydrosystem is formed by three major hydrologic compounds: the Haradros ephemeral stream, the karstic spring of Makaria that discharges in a drain canal and a naturally occurring wetland at the eastern part of the plain (Figure 2).
The land use in Marathon and the surrounding area is mainly agricultural, domestic and touristic.As a result of these activities, water demand highly varies all year round, with peak demand during summertime.The groundwater in the area is a valuable source The groundwater system consists of two aquifer units: a karstic aquifer developed in marble and a granular aquifer that overlies the marble formation developed in terrestrial, alluvial deposits (Figure 2).The impermeable bedrock of the karstic aquifer is composed of a schist layer.
Along with the aquifer formations, the Marathon hydrosystem is formed by three major hydrologic compounds: the Haradros ephemeral stream, the karstic spring of Makaria that discharges in a drain canal and a naturally occurring wetland at the eastern part of the plain (Figure 2).
The land use in Marathon and the surrounding area is mainly agricultural, domestic and touristic.As a result of these activities, water demand highly varies all year round, with peak demand during summertime.The groundwater in the area is a valuable source of freshwater, especially for the numerous farms and greenhouses that are irrigated entirely from groundwater resources.The uncontrolled abstraction-especially during the dry season-leads to the phenomena of aquifer quantity, quality deterioration and seawater encroachment on the mainland [27,28].Figure 3 presents the TDS concentration in the study area (October 2018).
of freshwater, especially for the numerous farms and greenhouses that are irrigated entirely from groundwater resources.The uncontrolled abstraction-especially during the dry season-leads to the phenomena of aquifer quantity, quality deterioration and seawater encroachment on the mainland [27,28].Figure 3 presents the TDS concentration in the study area (October 2018).

Modeling and Calibration Codes
Three different mathematical codes were utilized to investigate the seawater intrusion in the alluvial aquifer of Marathon.Firstly, the finite difference 3D code for groundwater movement in porous media, -MODFLOW-2005 [18], was applied.The partially differential equation that is executed by the code is defined as: where Kxx, Kyy and Kzz are values of hydraulic conductivity along the x, y and z coordinate axes (L/T); h is the potentiometric head (L); W is a volumetric flux per unit volume, representing sources and/or sinks of water, with W < 0.0 for flow out of the groundwater system and W > 0.0 for flow into the system (T − 1); SS is the specific storage of the porous material (L − 1); and t is time (T) [18].Finally, the SEAWAT v.4 code [14] was used in order to analyze groundwater movement under the effect of the variable density flow that is caused by seawater intrusion.The code integrates MODFLOW-2000 [29] and MT3DMS [19] codes simulating groundwater movement and contaminant transport, with modifications to address the problem of variable density flow in the aquifers.Further calibration for the hydraulic properties of the aquifer was not applied.The implementation of the SEAWAT code in this research serves two main purposes: (i) to explore the seawater intrusion mechanism in the several

Investigation of Seawater Intrusion Pattern and Mechanism through Modeling 2.2.1. Modeling and Calibration Codes
Three different mathematical codes were utilized to investigate the seawater intrusion in the alluvial aquifer of Marathon.Firstly, the finite difference 3D code for groundwater movement in porous media, -MODFLOW-2005 [18], was applied.The partially differential equation that is executed by the code is defined as: where K xx , K yy and K zz are values of hydraulic conductivity along the x, y and z coordinate axes (L/T); h is the potentiometric head (L); W is a volumetric flux per unit volume, representing sources and/or sinks of water, with W < 0.0 for flow out of the groundwater system and W > 0.0 for flow into the system (T − 1); SS is the specific storage of the porous material (L − 1); and t is time (T) [18].Finally, the SEAWAT v.4 code [14] was used in order to analyze groundwater movement under the effect of the variable density flow that is caused by seawater intrusion.The code integrates MODFLOW-2000 [29] and MT3DMS [19] codes simulating groundwater movement and contaminant transport, with modifications to address the problem of variable density flow in the aquifers.Further calibration for the hydraulic properties of the aquifer was not applied.The implementation of the SEAWAT code in this research serves two main purposes: (i) to explore the seawater intrusion mechanism in the several layers of the aquifer and (ii) to examine the performance of the model against the already calibrated MODFLOW model in terms of the hydraulic heads.
The ModelMuse platform of the USGS [30] was used to construct and simulate the model and visualize the model results.The SEAWAT code has not been integrated into ModelMuse yet, thus limiting the ability to simulate variable density flow problems using the finite difference method.Nevertheless, the generated files of both MODFLOW-2005 and MT3DMS can be used with minor modifications and additions so that it is possible to run the SEAWAT code externally (off-platform) [31].

Sensitivity Analysis and Calibration Code
For the purpose of this research, an existing model [26] was used, in which the sensitivity analysis and calibration were performed.The conceptualization of the model was not changed, but the model was updated with new hydraulic head measurements.For calibration purposes, we obtained hydraulic head data on a monthly timestep after a field campaign that lasted one hydrological year (October 2018-September 2019).The hydraulic head values were calculated by water level measurements (obtained with a water level meter) and absolute ground elevations (obtained with Differential GPS).For further information, interested readers are referred to [26].The sensitivity analysis took place in an automated way using the UCODE_2014 [32] coupled with MODFLOW [33].The code calculates parameters' sensitivity with fit-independent statistics in order to identify the importance of the observations for each parameter and the importance of a parameter for the simulation, as well as the importance of every observation for the simulation [34].The parameter sensitivity levels are calculated considering composite scaled sensitivities (CSS) that indicate the available information that observations offer to estimate a specific parameter.In conjunction with CSS, parameter correlation coefficients (PCC) are utilized to calculate the correlation between two parameters.In cases of extreme correlation (|PCC| > 0.95), the parameters cannot be calculated independently during the parameter estimation process [34].In the parameter estimation process, UCODE_2014 creates an objective function that is defined as a sum of the squared and weighted residuals of each observation minus the corresponding simulated value.The group of parameters that produces the minimum of the objective function forms the best prediction for the simulation [34].

Conceptual Model and Boundary Conditions
The conceptualization of the Marathon hydrosystem occurred through an extensive field investigation and review of past scientific reports [35][36][37].The research was focused on the eastern part of the plain-at the east of Haradros stream-due to the high chloride concentration in the area, as it was detected in previous field research and hydrochemical analysis [28].
The extent and boundaries of the model are formed within the boundaries of the alluvial deposits of the plain (Figure 2).The thickness of the deposits was determined by the geophysical explorations of Melissaris and Stavropoulos [35].Based on the data, the alluvial aquifer was divided into three sub-layers due to different grain size distribution with different hydraulic characteristics.The upper and lower layers comprise two different clay and sand-gravel zones, while the middle layer is composed of sand-gravel material.The extent of each zone is shown in Figure 4.The occurrence of different zones determines the hydrodynamic conditions and hydrochemical status of the aquifer.The thickness of each layer varies, although the formation does not exceed 60 m throughout the plain.Since no information is available for the extension of the aquifer offshore, the boundary of the system was defined along the coastline.
The coastal aquifer was simulated under transient conditions for 12 stress periods on a daily timestep.Each stress period represents a month of the hydrologic year, spanning from October 2018 to September 2019.The active area of the model is 19.67 km 2 , divided into a grid of 50 m × 50 m.This results in 7868 active cells per layer.

Neuman boundary conditions
• Abstraction wells and boreholes (Well package-WEL).Due to a lack of pumping rate data from the wells and drills, the pumping rate was determined based on data from the study of the Ministry of Agriculture [35] with an annual pumping rate of about 2 × 10 6 m 3 ; • Percolation from precipitation (Recharge package-RCH).The distribution of rainwater percolation throughout the year was calculated through a model by employing different hydrological methods [38].The model includes the Soil Conservation Service Curve Number (SCS-CN) approach [39], the Penman-Monteith method [40]) for Potential Evapotranspiration (PET), while the Soil Moisture Balance (SMB) [41] method was utilized to estimate the Actual Evapotranspiration and rainfall percolation in the aquifer.

Cauchy boundary conditions
• Inflows from the surrounding marble (General Head Boundary-GHB).For the boundary head of the GHB, several hydraulic head measurements were gathered during the field campaign (from October 2028 to September 2019).The hydraulic head of the karstic formation in the area varied from 2.5 m (during the dry season) up to 4 m (during the wet season); • Drain canal (Drain package-DRN).For the elevation of the drain canal in the area, a Digital Elevation Model (DEM) was utilized.

•
Haradros stream (River package-RIV).Since the stream has an ephemeral flow, a river stage of up to 1 m above the river bottom was utilized, according to the sparse observation data that were obtained during the field campaigns.The major hydrologic compounds of the area are represented by different boundary conditions (Figure 2), as listed below and presented in Figure 5:

Neuman boundary conditions
• Abstraction wells and boreholes (Well package-WEL).Due to a lack of pumping rate data from the wells and drills, the pumping rate was determined based on data from the study of the Ministry of Agriculture [35] with an annual pumping rate of about 2 × 10 6 m 3 ; • Percolation from precipitation (Recharge package-RCH).The distribution of rainwater percolation throughout the year was calculated through a model by employing different hydrological methods [38].The model includes the Soil Conservation Service Curve Number (SCS-CN) approach [39], the Penman-Monteith method [40]) for Potential Evapotranspiration (PET), while the Soil Moisture Balance (SMB) [41] method was utilized to estimate the Actual Evapotranspiration and rainfall percolation in the aquifer.

Cauchy boundary conditions
• Inflows from the surrounding marble (General Head Boundary-GHB).For the boundary head of the GHB, several hydraulic head measurements were gathered during the field campaign (from October 2028 to September 2019).The hydraulic head of the karstic formation in the area varied from 2.5 m (during the dry season) up to 4 m (during the wet season); • Drain canal (Drain package-DRN).For the elevation of the drain canal in the area, a Digital Elevation Model (DEM) was utilized.

•
Haradros stream (River package-RIV).Since the stream has an ephemeral flow, a river stage of up to 1 m above the river bottom was utilized, according to the sparse observation data that were obtained during the field campaigns.
The Head Observation package (HOB) was also utilized in order to assign the observed values to the model.
The Head Observation package (HOB) was also utilized in order to assign the observed values to the model.To develop the variable density flow model, the MODFLOW and MT3DMS models were employed.The contaminant transport model was integrated into the flow model with common spatial and temporal discretization, enabling the Sink & Source Mixing Package (SSM).For the simulation, a type of pollutant was defined at the coastline as a constant concentration boundary condition representing the Mediterranean Sea.Initial conditions for the contaminant transport model were set from the chemical analysis data in terms of total dissolved solids (TDSs).

Parametrization
An important task during the conceptualization of a flow model is the determination of the type and value of the parameters and their spatial distribution.For the model of Marathon, parametrization was applied in order to define the parameters of the aquifers' different zones and boundary conditions based on the collected data.Furthermore, the use of parameters is required for both the sensitivity analysis and calibration processes.In Table 1, the selected parameters for the aquifer of Marathon, as well as their initial values retrieved from past research [26,35,42], are presented.To develop the variable density flow model, the MODFLOW and MT3DMS models were employed.The contaminant transport model was integrated into the flow model with common spatial and temporal discretization, enabling the Sink & Source Mixing Package (SSM).For the simulation, a type of pollutant was defined at the coastline as a constant concentration boundary condition representing the Mediterranean Sea.Initial conditions for the contaminant transport model were set from the chemical analysis data in terms of total dissolved solids (TDSs).

Parametrization
An important task during the conceptualization of a flow model is the determination of the type and value of the parameters and their spatial distribution.For the model of Marathon, parametrization was applied in order to define the parameters of the aquifers' different zones and boundary conditions based on the collected data.Furthermore, the use of parameters is required for both the sensitivity analysis and calibration processes.In Table 1, the selected parameters for the aquifer of Marathon, as well as their initial values retrieved from past research [26,35,42], are presented.

Sensitivity Analysis of MODFLOW Parameters
The results of the sensitivity analysis for the MODFLOW model were analyzed based on the CSS and PCC indices calculated by UCODE2014.
For the alluvial aquifer model, the results of the sensitivity analysis suggest that the model is highly influenced by the parameters related to the characteristics of the middle layer and sand-gravel zones, as well as the pumping parameters.The high values of the CSS index for the parameters HK2, SS_Par2, Q_KS and Q_M, as shown in Figure 6a, are partly due to the effect of pumping that takes place mainly in the middle layer formed entirely of sand-gravel.The specific yield (Sy) of the upper layer is less important for the simulation, while the parameters concerning the clay zones determine the result to an even lesser extent.The parameters concerning the boundary conditions and the model outputs (Cond_North, Cond_West and Cond_drn) have no significant sensitivity.A high correlation is observed between the Kato Souli pumping parameter (Q_KS) and the hydraulic properties of the sand-gravel (HK2, SS_Par2) and a smaller one with the Marathon pumping (Q_M) related to the effect of pumping on the middle layer, as described above and shown in Figure 6b.

Sensitivity Analysis of MODFLOW Parameters
The results of the sensitivity analysis for the MODFLOW model were analyzed based on the CSS and PCC indices calculated by UCODE2014.
For the alluvial aquifer model, the results of the sensitivity analysis suggest that the model is highly influenced by the parameters related to the characteristics of the middle layer and sand-gravel zones, as well as the pumping parameters.The high values of the CSS index for the parameters HK2, SS_Par2, Q_KS and Q_M, as shown in Figure 6a, are partly due to the effect of pumping that takes place mainly in the middle layer formed entirely of sand-gravel.The specific yield (Sy) of the upper layer is less important for the simulation, while the parameters concerning the clay zones determine the result to an even lesser extent.The parameters concerning the boundary conditions and the model outputs (Cond_North, Cond_West and Cond_drn) have no significant sensitivity.A high correlation is observed between the Kato Souli pumping parameter (Q_KS) and the hydraulic properties of the sand-gravel (HK2, SS_Par2) and a smaller one with the Marathon pumping (Q_M) related to the effect of pumping on the middle layer, as described above and shown in Figure 6b.

MODFLOW Model Calibration and Results
Parameter estimation for the MODFLOW model was performed with respect to the sensitivity analysis results.The parameters selected for the UCODE2014 parameter estimation process concern the aquifer characteristics HK2, SS_Par2, Sy, HK1 and SS_Par1.The parameters Q_M and Q_KS were not included in the model calibration, and the initial values were used for the entire simulation.These parameters were primarily utilized to investigate the impact of withdrawals to the aquifer.The calibrated parameters are presented in the following table (Table 2).The model improved significantly with respect to the measured hydraulic head after the calibration process.Furthermore, the statistics of the differences observed vs. the simulated hydraulic head, shown in Figure 7, present an acceptable agreement between the simulated and observed hydraulic head values.

MODFLOW Model Calibration and Results
Parameter estimation for the MODFLOW model was performed with respect to the sensitivity analysis results.The parameters selected for the UCODE2014 parameter estimation process concern the aquifer characteristics HK2, SS_Par2, Sy, HK1 and SS_Par1.The parameters Q_M and Q_KS were not included in the model calibration, and the initial values were used for the entire simulation.These parameters were primarily utilized to investigate the impact of withdrawals to the aquifer.The calibrated parameters are presented in the following table (Table 2).The model improved significantly with respect to the measured hydraulic head after the calibration process.Furthermore, the statistics of the differences observed vs. the simulated hydraulic head, shown in Figure 7, present an acceptable agreement between the simulated and observed hydraulic head values.

Parameter Fitting in SEAWAT Model
The parameters obtained from the MODFLOW-2005 calibration were used in the variable density flow model, SEAWAT.The Observation (OBS) Process, as described by Langevin et al. [43], was activated manually to test the fitting of the parameters in the SEAWAT model.The simulation results, in relation to the observed hydraulic head values, are presented in Figure 8.An acceptable agreement between the simulated and observed

Parameter Fitting in SEAWAT Model
The parameters obtained from the MODFLOW-2005 calibration were used in the variable density flow model, SEAWAT.The Observation (OBS) Process, as described by Langevin et al. [43], was activated manually to test the fitting of the parameters in the SEAWAT model.The simulation results, in relation to the observed hydraulic head values, are presented in Figure 8.An acceptable agreement between the simulated and observed hydraulic head values was achieved in the SEAWAT model by utilizing the hydraulic parameters that were previously calibrated in the MODFLOW model.In fact, there is a slightly better match between the observed and simulated values with the SEAWAT model, as indicated by the metrics in Figure 8.The use of the SEAWAT model against the MODFLOW model in terms of the hydraulic head produced a greater match between the simulated and measured values and improved the results of the already calibrated flow model (MODFLOW-2005) without further calibration.The main difference between the two models is reflected in the near-the-shoreline observed vs. simulated values where the root-mean-square error (RMSE) is equal to 0.58 for the MODFLOW model, while for the SEAWAT model, the RMSE decreases to 0.50.
hydraulic head values was achieved in the SEAWAT model by utilizing the hydraulic parameters that were previously calibrated in the MODFLOW model.In fact, there is a slightly better match between the observed and simulated values with the SEAWAT model, as indicated by the metrics in Figure 8.The use of the SEAWAT model against the MODFLOW model in terms of the hydraulic head produced a greater match between the simulated and measured values and improved the results of the already calibrated flow model (MODFLOW-2005) without further calibration.The main difference between the two models is reflected in the near-the-shoreline observed vs. simulated values where the root-mean-square error (RMSE) is equal to 0.58 for the MODFLOW model, while for the SEAWAT model, the RMSE decreases to 0.50.In addition, the simulation of the variable density flow model with the SEAWAT model produced different results in the distribution of the piezometric curves with respect to the MODFLOW model.The maps of iso-piezometric curves for the two models are presented in Figures 9 and 10, for the simulated hydrological year (2018-2019).The main difference between the two models lies in the near-shoreline piezometric contours that are shifted further inland.This shift is due to the effect of the transition zone that cannot be captured by a non-variable density model.The difference between the two simulations is more pronounced during the period of extensive pumping, i.e., between July 2019 and September 2019.
Regarding the transport model parameters, no automatic calibration was performed due to the limited data.However, the manual calibration of longitudinal dispersivity and porosity was applied by assigning different parameter values utilized in the existing research [44][45][46][47][48][49].The parameter calibration process revealed that no significant improvement in simulated-observed TDS values occurred.Nevertheless, there was an acceptable agreement between the observed vs. simulated TDS values (R 2 = 0.9).In addition, the simulation of the variable density flow model with the SEAWAT model produced different results in the distribution of the piezometric curves with respect to the MODFLOW model.The maps of iso-piezometric curves for the two models are presented in Figures 9 and 10, for the simulated hydrological year (2018-2019).The main difference between the two models lies in the near-shoreline piezometric contours that are shifted further inland.This shift is due to the effect of the transition zone that cannot be captured by a non-variable density model.The difference between the two simulations is more pronounced during the period of extensive pumping, i.e., between July 2019 and September 2019.
Regarding the transport model parameters, no automatic calibration was performed due to the limited data.However, the manual calibration of longitudinal dispersivity and porosity was applied by assigning different parameter values utilized in the existing research [44][45][46][47][48][49].The parameter calibration process revealed that no significant improvement in simulated-observed TDS values occurred.Nevertheless, there was an acceptable agreement between the observed vs. simulated TDS values (R 2 = 0.9).

The Mechanism of Seawater Intrusion
In order to investigate further the mechanism of seawater intrusion in Marathon, the model was simulated for a longer period to achieve the natural mixing between the aquifer concentration and the seawater.The main goal was to investigate and determine the shape and the extent of the saltwater wedge in the coastal granular formation.According to Mabrouk et al. [50], two approaches are commonly used to determine the distribution of salinity in the aquifer: (a) assigning the initial salinity distribution spatially based on existing salinity data or (b) starting the simulation with a fully saline or fully freshwater

The Mechanism of Seawater Intrusion
In order to investigate further the mechanism of seawater intrusion in Marathon, the model was simulated for a longer period to achieve the natural mixing between the aquifer concentration and the seawater.The main goal was to investigate and determine the shape and the extent of the saltwater wedge in the coastal granular formation.According to Mabrouk et al. [50], two approaches are commonly used to determine the distribution of salinity in the aquifer: (a) assigning the initial salinity distribution spatially based on existing salinity data or (b) starting the simulation with a fully saline or fully freshwater aquifer and simulating for an extended period until the groundwater system reaches a dynamic equilibrium of a freshwater-seawater interface.Gomaa et al. [51], in their research in a coastal granular aquifer of Egypt, concluded that in the second case, freshwater occupied the development area, and the result did not match the salinity data.In contrast, the results of the first approach were in agreement with the observed values.Based on this conclusion and considering the existing salinity data, the first approach was adopted for the case of Marathon.
A reliable steady-state SEAWAT model was initially developed.To achieve steadystate conditions, the SEAWAT model was run for a "warm-up" period of 3000 timesteps.It must also be noted that the model was also run for a "warm-up" period of 10,000 timesteps, but the results revealed that steady-state conditions for the study area were achieved after 2000 timesteps.To this end, in order to reduce the computational burden, it was chosen to proceed with the modeling exercise using the initial model with the 3000 timesteps.The variable density flow was then simulated for 1827 timesteps under transient flow conditions for 5 hydrological years (2015-2020).
In Figure 11, cross-sections with the distribution of total dissolved solids (TDS) are plotted in the most affected areas of the aquifer as a result of the SEAWAT simulation.
Hydrology 2024, 11, x FOR PEER REVIEW 13 of 18 aquifer and simulating for an extended period until the groundwater system reaches a dynamic equilibrium of a freshwater-seawater interface.Gomaa et al. [51], in their research in a coastal granular aquifer of Egypt, concluded that in the second case, freshwater occupied the development area, and the result did not match the salinity data.In contrast, the results of the first approach were in agreement with the observed values.Based on this conclusion and considering the existing salinity data, the first approach was adopted for the case of Marathon.
A reliable steady-state SEAWAT model was initially developed.To achieve steadystate conditions, the SEAWAT model was run for a "warm-up" period of 3000 timesteps.It must also be noted that the model was also run for a "warm-up" period of 10,000 timesteps, but the results revealed that steady-state conditions for the study area were achieved after 2000 timesteps.To this end, in order to reduce the computational burden, it was chosen to proceed with the modeling exercise using the initial model with the 3000 timesteps.The variable density flow was then simulated for 1827 timesteps under transient flow conditions for 5 hydrological years (2015-2020).
In Figure 11, cross-sections with the distribution of total dissolved solids (TDS) are plotted in the most affected areas of the aquifer as a result of the SEAWAT simulation.In cross-section A-A', a typical seawater intrusion wedge is formed.At this location, the extent of the freshwater-saltwater interface increases with depth toward inland.In the central part of the plain, the TDS concentration increases with depth up to 4 g/m 3 .The main cause of high TDS concentration is attributed to the cone of depression that appears during and after the dry period, a fact that is also validated by the chemical analysis and simulation results.During that season, extensive pumping led to upward movement of the freshwater-saltwater interface.The cone of depression is also validated by the hydraulic head in the area during the period of extensive pumping, i.e., from July to October (Figure 10).The upward movement of the interface is reflected in cross-section A-A as the In cross-section A-A', a typical seawater intrusion wedge is formed.At this location, the extent of the freshwater-saltwater interface increases with depth toward inland.In the central part of the plain, the TDS concentration increases with depth up to 4 g/m 3 .The main cause of high TDS concentration is attributed to the cone of depression that appears during and after the dry period, a fact that is also validated by the chemical analysis and simulation results.During that season, extensive pumping led to upward movement of the freshwater-saltwater interface.The cone of depression is also validated by the hydraulic head in the area during the period of extensive pumping, i.e., from July to October (Figure 10).The upward movement of the interface is reflected in cross-section A-A as the TDS concentration is higher in the bottom layer (4 g/m 3 ), and it is gradually decreased in the top layer (2.5 g/m 3 ).Salinity is reduced toward the north where the karstic marble formation of lower salinity recharges the alluvial aquifer.
In cross-section B-B, the interface is moved further inland across the coastline, while in the northern part, the TDS concertation is lower than in cross-section A-A.The decrease in salinity is due to the seizing of groundwater withdrawals since the specific area is no longer cultivated.
Cross-section C-C presents the TDS concentration in the area of the wetland, where seawater intrusion extends further inland than in the western part of the plain.The extended transition zone is caused by the low topography of the wetland, as also stated by Margoni [37], in conjunction with the negative values of the hydraulic head that lower the hydraulic pressure in the aquifer, leading to extended seawater intrusion.In addition, the aquifer geological heterogeneity near the coastline (Figure 4) could have an effect on the extent and irregular shape of the transition zone.

Discussion
During the present research, an important issue arose related to the calibration process of the model since three different codes were implemented for the simulation of the coastal aquifer.Generally, for hydraulic parameter calibration, the MODFLOW model is utilized without taking into account the effects of the variable density flow that the SEAWAT code incorporates.Abd-Elaty and Zelenakova [47], in their research on the Nile Delta aquifer, calibrated the SEAWAT model with the trial-and-error method with respect to the TDS concentration.Ansarifar et al. [25] simulated the effect of sea level fluctuation with a calibrated MODFLOW model.Lyra et al. [24] utilized the MODFLOW model and PEST analysis to simulate and calibrate the hydraulic heads of a coastal aquifer, while solute transport parameters were calibrated through trial and error in the SEAWAT model.Nasiri et al. [52] in their research on seawater intrusion mitigation in coastal aquifers, calibrated hydraulic heads in the MODFLOW model and TDS concentration in the MT3DMS model before simulating the variable density flow with SEAWAT.In our research, the hydraulic heads were calibrated in the MODFLOW model and the TDS concentration in the SEAWAT model.In addition to methods used in the aforementioned studies, the OBS process was activated in the SEAWAT model [43].This procedure revealed that the hydraulic parameters calibrated in MODFLOW produced better results for the observed vs. simulated head values of the SEAWAT model without further calibration by capturing the hydraulic head differences that the constant density flow model (MODFLOW) could not capture.Furthermore, this research reveals that the distribution of the hydraulic heads with SEAWAT differs from the piezometry produced by the MODFLOW model, especially in the areas of high TDS concentration (near the shoreline and at the wetland), where variable density effects are more pronounced.SEAWAT, the density-variant model, can be compared with MODFLOW, the density-invariant model, in which concentration or temperature changes do not affect the density.The results of the two cases are very similar or identical (Figures 7 and 8) for the upstream part, where it can be observed that density is likely an insignificant part of the system behavior in that part of the two models.However, there is a significant difference between the two cases within the coastal zone, where it is proved that density is an important factor.As Post and Simmons [53] quote, the convective groundwater systems are mixed convection systems where forced convection (flow driven by pressure or hydraulic head gradients) and free convection (flow driven by density variations) co-exist.
As for the solute transport model parameters, in our case study, longitudinal dispersivity was of low importance for the TDS concentration, a result that is in accordance with Shoemaker's research.Further validation of the model results can be carried out by obtaining hydraulic head and flow data near the coastline since effective locations for the calibration of SEAWAT with hydraulic head observations are the toe of the transition zone and areas near the coastline [54].Unfortunately, during field research in the area, it was not possible to locate wells near the coastline, possibly due to the seizing of groundwater abstractions as a result of seawater intrusion.Furthermore, for the specific case study, the vertical discretization was decided with respect to the different hydrogeological zones of the alluvial aquifer, as indicated in Figure 4. Further vertical discretization would also improve the model.In order to apply vertical discretization and validate the produced results, an advanced research infrastructure (e.g., multilevel piezometers) is essential.As a result, this topic is reserved for future work.
Another outcome of the present research is the investigation of the seawater intrusion mechanism in the alluvial aquifer of Marathon.Considering the distribution and the shape of the hydraulic head in the aquifer (Figures 9 and 10), the results of the simulation reveal that the aquifer is subjected to passive-active seawater intrusion, as analyzed by Werner [55].In this situation, groundwater flows both inland and toward the sea.The phenomenon is more pronounced in the vicinity of the low-topography wetland, while at the central part of the plain, there is a temporal variability of the phenomenon.During the wet season, the hydraulic gradient is seaward, and a typical case of passive seawater intrusion occurs, while during the period of extensive pumping (July-October), the negative values of the hydraulic head increase the landward movement of seawater.
The extent and shape of the transition zone, as affected by the heterogeneity that characterizes the coastal formations, has been investigated through mathematical models by many researchers (e.g., refs.[56][57][58]), with the main conclusion that in heterogeneous porous formations, the freshwater-saltwater transition zone intensifies as a result of irregular flow conditions [1].Across the wetland area, the transition zone does not form the typical saltwater wedge.In the middle layer of the aquifer, a much wider transition zone develops than in the lower layer, as it is affected by the different hydraulic properties of the aquifer's zones.At the bottom layer, the zone of lower hydraulic conductivity (Figure 4c) prevents a faster development of seawater encroachment.The results are in accordance with Michael et al.'s [59] conclusion that geological heterogeneity has such a strong influence on the flow conditions that the salinity distribution and the shape of the freshwater-saltwater interface have an irregular and complex form, with no resemblance to the wedge-shaped interface consideration.

Conclusions
The present study involves the evaluation of variable density flow and assessment of the seawater intrusion mechanism in a coastal aquifer through modeling.The study area is a typical coastal Mediterranean hydrosystem located in Marathon, Greece.
The main conclusions derived from the use of MODFLOW and SEAWAT models to simulate groundwater flow are as follows: (i) the sensitivity analysis performed on the MODFLOW model with the UCODE_2014 code indicates that the most sensitive parameters are the hydraulic conductivity (HK) of the main aquifer layer (middle layer of the model), as well as the parameter of pumping (Q), concluding that the parameters related to the particular characteristics of the aquifer and the total water use (withdrawals) in the watershed have the greatest influence on the simulation of the aquifer, (ii) the use of the SEAWAT code to simulate flow under variable density conditions produced a greater match between the simulated and observed values of the hydraulic heads and, in this sense, further improved the results of the already calibrated flow model (MODFLOW-2005) and (iii) the distribution of hydraulic heads in the area produced with the MODFLOW model differ significantly with those produced with the SEAWAT model in the areas of high TDS concentration.Additionally, it appears that longitudinal dispersivity is not sensitive to the simulation.Regarding the "porosity" parameter, it appears that this parameter is more important for the simulation results in the variable density flow model than the "longitudinal dispersivity" parameter.However, the model results seem to be affected at very small porosity values (i.e., 0.05).
Parameter transferability seems to work, as even without calibration, the results of the SEAWAT model indicate a slightly better performance; thus, parameter transferability can take place between the two models.The limitations of this study lie in the lack of observed hydraulic head data near the coastline that would further validate the proposed framework since the most effective observation locations for the calibration of seawater intrusion models are those that are located near the coastline.
Finally, the mechanism of seawater intrusion in the area was assessed through the TDS concentration results produced from the SEAWAT model.The results indicate that the cone of depression at the central part of the plain is responsible for the high TDS concentration throughout the year.Furthermore, in the eastern part of the plain, the heterogeneity that characterizes the coastal formation in terms of the hydraulic parameters creates an irregular shape of the saltwater wedge in the specific area.

Figure 4 .
Figure 4. Layer zonation in the (a) upper aquifer layer, (b) middle aquifer layer and (c) lower aquifer layer.The coastal aquifer was simulated under transient conditions for 12 stress periods on a daily timestep.Each stress period represents a month of the hydrologic year, spanning from October 2018 to September 2019.The active area of the model is 19.67 km 2 , divided into a grid of 50 m × 50 m.This results in 7868 active cells per layer.The major hydrologic compounds of the area are represented by different boundary conditions (Figure 2), as listed below and presented in Figure 5: 1. Dirichlet boundary conditions • Coastline (Time-Variant Specified-Head package-CHD) at 0 m; • Makaria Spring (Time-Variant Specified-Head package-CHD) at 2 m; • Coastal Wetland (Time-Variant Specified-Head package-CHD) ranging from −1.2 m to 0.6 m to obtain observations in the vicinity of the wetland.

Figure 4 .
Figure 4. Layer zonation in the (a) upper aquifer layer, (b) middle aquifer layer and (c) lower aquifer layer.

1 .
Dirichlet boundary conditions • Coastline (Time-Variant Specified-Head package-CHD) at 0 m; • Makaria Spring (Time-Variant Specified-Head package-CHD) at 2 m; • Coastal Wetland (Time-Variant Specified-Head package-CHD) ranging from −1.2 m to 0.6 m to obtain observations in the vicinity of the wetland.

Figure 7 .
Figure 7. Observed vs. simulated hydraulic head values for the MODFLOW model.

Figure 7 .
Figure 7. Observed vs. simulated hydraulic head values for the MODFLOW model.

Figure 8 .
Figure 8. Observed vs. simulated hydraulic head values for the SEAWAT model.

Figure 8 .
Figure 8. Observed vs. simulated hydraulic head values for the SEAWAT model.

Figure 11 .
Figure 11.Aquifer TDS concentration in cross-section as a result of SEAWAT model.

Figure 11 .
Figure 11.Aquifer TDS concentration in cross-section as a result of SEAWAT model.