Implementation of a Far-Field Water Quality Model for the Simulation of Trace Elements in an Eastern Mediterranean Coastal Embayment Receiving High Anthropogenic Pressure

: Water quality modeling is a key element for the support of environmental protection and policymaking. The aim of this work is to describe the application of a far-field water quality model for the simulation of marine pollution occurring from heavy metals (cadmium, lead, nickel, copper, and zinc). The highly stressed marine area of the Saronikos Gulf (Aegean Sea, Eastern Mediterranean) was chosen for investigation. Major pollution sources were identified, loads were estimated, and the model was parameterized to reproduce the local seawater conditions. The distribution of the pollutants between the dissolved and particulate phases was examined. The performance of the model set-up was evaluated using field concentration measurements. The described implementation succeeded in reproducing the observed levels of pollution and therefore can be used as a baseline configuration to examine the cumulative impact of future pollution sources; for example, accidental pollution events.


Introduction
In a rapidly changing world, sustainable development requires careful planning and early impact prevision [1].To protect and preserve natural water resources and ecosystems it is important to develop and optimize tools for the evaluation and prediction of anthropogenic pressures [2,3].Water quality modeling is an implementation moving towards that direction, used in various applications and for all types of aquatic environments [4].Some fields in which water quality models are used are: human health protection from pathogen pollution [5]; prediction of eutrophication [6]; prediction and investigation of chemical pollution, including heavy metals [7], emerging contaminants, and antibiotics [8,9]; antibiotic resistance of aquatic environments [10]; and many others.The most crucial and valuable usage of water quality modeling as a tool is the provision of support for the development of policies and environmental management strategies [11,12].In this context, the development and validation of reliable models are of high importance, especially in areas where anthropogenic pressure is intense.
A marine environment receiving high anthropogenic pressure is the water body alongside the coastal area of Attica (Greece), an area populated by almost 4 million inhabitants [13].The Saronikos Gulf (Figure 1) is situated in the Eastern Mediterranean, in the Aegean Sea, and is the water body adjacent to the capital of Greece, Athens.The water quality of the Saronikos Gulf, over the last decades, has varied depending on human activities across different transitional eras.Anthropogenic pollution started in the early 1920s, with the installation of the first industries in the Thriassion Plain and the construction of central sewage pipes to serve the growing population of Attica.During the late 1970s and early 1980s, the pollution occurring from industrial activities close to Elefsis Bay, as well as sewage outflow of untreated wastewater of Athens, drew the attention of the scientific community [14][15][16].During the 1990s, the Saronikos Gulf reached high pollution levels, and at the same time, measures to mitigate further water quality deterioration were taken.The primary operation of the Psittalia wastewater treatment plant (WWTP) in 1994 partially alleviated the pollution levels until 2007, when the plant became fully operational [17][18][19].In the same time period, the industrial activities in the area of Elefsis started declining and, in parallel, industries gradually adopted measures for pollution reduction [20].Since the mid-1980s, the Saronikos Gulf has been monitored frequently under the aegis of the projects 'National Monitoring Program for the Assessment and Control of Marine Pollution in the Mediterranean' (MED-POL) MAP/UNEP and 'Monitoring of the Saronikos Gulf ecosystem under the influence of the Psittalia sewage outflow' [21], implemented by the Hellenic Centre for Marine Research (HCMR).During the last decade, the requirements of the European Water Framework Directive and the European Marine Strategy Framework Directive have driven the continuation of the monitoring effort.
Although a lot of studies focus on the state of the marine environment of the Saronikos Gulf, only a few relevant numerical modeling studies have been published for this specific area.Some of the published modeling studies focus on ecology, with an interest in the food web carbon flow and interactions of individual species [22][23][24].Only a few articles discuss the hydrodynamics of the Saronikos Gulf.Kontoyiannis [25] presented, for the first time, the actual three-dimensional flows in the Saronikos Gulf based on direct current measurements on a seasonal basis and under various wind conditions, while Belibassakis and Karathanasi simulated the nearshore hydrodynamics of Varkiza coast [26].Similarly, relatively few articles implement or develop models related to pollution for the case study area.Microplastic fate has been modeled in the past [27,28], and recently, Makatounis et al. [29] have modeled the oil spill of the Agia Zoni II tanker shipwreck, an accident that happened in 2017 and had a high environmental impact on the region.To the best of our knowledge, there are no published studies employing a model targeting chemical pollution.
The aim of this study is to describe the implementation of a far-field water quality model covering the Saronikos Gulf.In the following sections, we describe the pollutant sources and target compounds of this study, the set-up of the numerical models used, the simulated processes, the sources and sinks of pollutants, and the statistical tools that are used for the model assessment under the Material and Methods section.The following sections refer to the evaluation of the quality of the simulation, based mainly on comparison with available pollutant measurements, and to the spatiotemporal variability and partitioning of the selected compounds as reproduced by the model.The results and prospects are discussed at the end of this work.

Study Area and Target Compounds
The area of study is located in the Eastern Mediterranean, in the Aegean Sea.The waterbody of the Saronikos Gulf, enclosed by the western coast of the Attica peninsula and the Northeast coast of the Peloponnese, was selected for model implementation (Figure 1).The region under investigation has a total marine surface area of approximately 2800 km 2 and a water volume of approximately 3 × 10 12 m 3 .The study area is divided into the following subareas, based on bottom morphology: (i) Elefsis Bay, connected with the adjacent sub-basins through two narrow and shallow straits, located in the west (Megara strait) and the southeast (Keratsini strait), (ii) eastern Saronikos basin that is further divided into the inner Saronikos to the north and the outer Saronikos to the southeast, from where the gulf connects with the Aegean Sea, and (iii) western Saronikos Gulf (Table S1).The network of monitoring stations (S1, S2, S3, S7, S8, S11, S13, S16, S18, S43) as implemented in the framework of the action 'Monitoring of the Saronikos Gulf ecosystem under the influence of the Psittalia sewage outflow' [21], as well as the sources of pollution (1-9), are displayed in Figure 1b.
(Figure 1).The region under investigation has a total marine surface area of approximately 2800 km 2 and a water volume of approximately 3 × 10 12 m 3 .The study area is divided into the following subareas, based on bottom morphology: (i) Elefsis Bay, connected with the adjacent sub-basins through two narrow and shallow straits, located in the west (Megara strait) and the southeast (Keratsini strait), (ii) eastern Saronikos basin that is further divided into the inner Saronikos to the north and the outer Saronikos to the southeast, from where the gulf connects with the Aegean Sea, and (iii) western Saronikos Gulf (Table S1).The network of monitoring stations (S1, S2, S3, S7, S8, S11, S13, S16, S18, S43) as implemented in the framework of the action 'Monitoring of the Saronikos Gulf ecosystem under the influence of the Psittalia sewage outflow' [21], as well as the sources of pollution (1-9), are displayed in Figure 1b.Five trace elements were chosen to set up and validate a water quality model configuration: cadmium (Cd), lead (Pb), nickel (Ni), copper (Cu), and zinc (Zn).From the above-mentioned elements, Cd, Pb, and Ni are included in the list of priority substances of the European Water Framework Directive (2000/60/EC), while Cu and Zn are known to pose ecological risks in the aquatic environment, e.g., [30].The selection of the trace elements was made taking into consideration the capability of the selected model (Delft3D-WAQ) to explicitly simulate these compounds, as well as the available data measurements for validating the simulations, originating from Paraskevopoulou et al.'s [21] dataset (hereafter P-2014 data).

Software
The three-dimensional modeling suite Delft3D was used to investigate hydrodynamics, water quality, and pollution.The software has been developed by Deltares (Delft, The Netherlands) and has been used in many research applications over the past 15 years [31][32][33][34].The source code is freely available and distributed by Deltares (https://oss.deltares.nl(accessed on: 25 November 2021)) upon registration.

Domain and Simulation Period
A curvilinear grid (Figure 1) was developed with a horizontal resolution of 1400 m close to the open boundaries (outer Saronikos) and higher resolution in the inner Saronikos and Elefsis Bay (300-700 m).In the vertical direction, 15 σ layers were chosen to resolve the water column.The simulation period, starting on 01/11/2009 and terminating on 01/11/2010, was selected based on in situ data availability, which would be required for the validation of the model (Section 2.3).

Modeling Approach
In the current study, modules Delft3D-FLOW and Delft3D-WAQ were coupled offline to simulate the five selected trace elements (Cd, Pb, Ni, Cu, Zn).The hydrodynamics module Delft3D-FLOW has been previously developed and is described elsewhere [35,36].The fate of the selected heavy metals relies on partitioning between seawater and particulate matter and on transport [37].Their concentration change in each computational cell is calculated according to Equation (1).
where, C is the concentration of a trace element (g m −3 ), t is the time (day), and on the right hand side, Q loads is the time-dependent trace element loading from the various sources, the second term in the parentheses accounts for the contribution from the transport of material by currents (advective and turbulent diffusive fluxes), F res is the rate of change of concentration due to resuspension fluxes of material in the bottom cells, and the last term corresponds to the change due to settling of the various compounds.The three-dimensional current velocity is described by u (bold denoting a vector quantity), k is the eddy diffusivity, and w settling is the settling velocity.All terms on the right-hand side are expressed in units of g m −3 day −1 .The most critical component strongly influencing the concentration of target pollutants is "loads", which stands for the mass loads of pollutants entering the domain during the simulation period.The definition of pollution sources is of high difficulty and complexity and the adopted methodology is described in Section 2.2.5.As soon as pollutants enter the modeling area, they are subjected to partitioning processes.In this study, both the dissolved and particulate phases of trace elements were considered in the WAQ module.Regarding the "transport" of pollutants, advection and dispersion were computed via the FLOW module (Section 2.2.4), affecting both the particulate and the dissolved phase of each element.Trace elements are adsorbed to suspended matter, depending on the sorption affinity of each element and the particulate matter availability and quality, and are deposited on the sediment layer through the 'settling' process and subjected to sedimentation (see Sections 2.2.7 and 2.2.8).The concentration of trace elements in the water column is also influenced by the process of sediment 'resuspension', as fine-grained and loosely consolidated sediments are brought to resuspension [37].The process of sediment resuspension can be of high importance, mainly in shallow waters, e.g., [4]; therefore, this should be influencing only a very small part of the modeling area.
In the current modeling setup, the seabed was treated as a single sediment layer (1 cm thickness), where the flux of trace elements was only considered from the water column towards the sediment surface, due to 'settling'.Exchange fluxes due to the resuspension of sediments as well as direct 'adsorption to' and 'desorption from' the sediment layer were not considered.
With the use of Delft3D-WAQ, it is possible to consider multiple sediment layers using a computational grid of the active sediment to simulate the trace elements transport between sediment and overlying water, as well as between different sediment layers [38].However, to achieve this, a comprehensive set of information on sediment composition and quality of the modeling area is required, which unfortunately was not available.

Hydrodynamic Circulation (Delft3D-FLOW)
The module Delft3D-FLOW was previously applied to simulate the hydrodynamic conditions of the Saronikos Gulf during the Nov. 2009-Nov.2010 annual cycle [35,36].The model solves a system of non-linear equations consisting of the horizontal equations of motion, the continuity equation, and the transport equations for conservative constituents, derived from the three-dimensional Navier-Stokes equations for incompressible, free surface flow.For the application, relevant forcing fields (air temperature, air pressure, relative humidity, precipitation, cloudiness, net shortwave radiation) with hourly data were used from the ECMWF ERA5 database [39].Furthermore, freshwater sources were considered (Kifisos river, Mandra and Sarantapotamos intermittent water streams, and discharges from Psittalia and Thriassion WWTPs).Annual flows were determined with the use of measurements retrieved from databases (WWTPs online Monitoring Database) and prediction tools [40].Open boundary conditions were derived from a 30-year hindcast study of the coupled Eastern Mediterranean-Black Sea system employing the free surface, hydrostatic, primitive equation Regional Ocean Modelling System [41].

Process Description and Parameterization (Delft3D-WAQ)
For the reproduction of trace element fractionation between the dissolved and particulate phase, relevant water quality variables were simulated.With the use of the Delft3D-WAQ module, organic matter (POC, DOC, DON, DOP), dissolved inorganic nutrients (NH 4 + , NO 3 − , PO 4 3− , Si), and phytoplankton (diatoms, non-diatoms) were simulated.For these variables, the relevant processes were selected (Tables 1 and S2) and the parameter values were determined accordingly to be representative of the modeling area (Table S3).Similarly, the relevant nutrient loads entering the system were estimated and determined (Table S4).The inorganic matter was introduced as a time series, based on field observations [42] (Table S5).

Trace Element Loads
In the current study, pollution loads were restricted to land-based major pollution sources.These mass loads were introduced in the modeling area as "discharges" in the respective locations (Figure 1).In Table 2, the polluting source name, the water volume discharge, as well as trace element concentration data are displayed.For the river of Kifisos and the WWTPs of Psittalia and Thriassion, monitoring [43] and published data [17] were used.Regarding the concentration of pollutants in industrial wastewater and intermittent water streams, actual monitoring data were not available.Therefore, assumptions had to be made, taking into consideration permitted discharge limits for wastewater and drinking water, respectively, as defined by EU (Directive 2010/75/EU) and national legislation [44].Furthermore, the contribution of shipyard activities was taken into consideration, with the use of the MAMPEC software (Deltares V2.0) for the area of Perama, where a small-scale ship repairing zone exists (Figure S1 and description in Supplementary Materials).In this research, the atmospheric deposition of trace elements, as well as their entrance into the system through groundwater and other non-point sources were not examined.

Trace Element Partitioning
Adsorption was considered to take place in particulate inorganic matter (IM), non-living particulate organic matter (POC), and dissolved organic matter (DOC), and phytoplankton; in addition, the complexation of trace elements with dissolved organic matter (e.g., proteins, humic substances) is also taken into consideration [52].Although the adsorption/desorption is a dynamic process and is more efficiently modeled using kinetic rates [53], partitioning was simulated according to the equilibrium approach while the adsorption capacity of each element was described with the use of seasonal partitioning coefficients, relevant to the case study area.Partitioning coefficients were calculated according to the following equation (Equation (2), Table S6), using field concentrations of trace elements (P-2014 data) in both the dissolved and particulate phase (adsorbed to suspended solids).
where K p : the partition coefficient of a trace element (mg kg −1 /mg L −1 ); M p : the concentration of a trace element in suspended solids (mg kg −1 ); M d : the concentration of a trace element in water (mg L −1 ); C p : the particulate concentration of a trace element (mg L −1 ); C d : the dissolved concentration of a trace element (mg L −1 ); C ss : the concentration of suspended solids (mgSS L −1 ).
A dynamic modeling approach of this process could not be selected, as it was not possible to determine the necessary site-specific kinetic rate values that would include all possible reaction mechanisms influencing the adsorption/desorption process of trace elements.Conversely, the use of partitioning coefficients, calculated from suspended solids rather than sediments concentrations, eliminates the risk of underestimating this variable [54].By using the equilibrum modeling aproach, the following risks may be anticipated: (i) underestimation of trace elements in the sediments in the long term, (ii) overestimation of the mobility of pollutants by underestimating the interaction with sediments close to the release source, and (iii) omission of the role of the sediment, acting as a pollution source through desorption and release of pollutants [53].

Trace Element Settling
Trace elements were considered to settle with the carrier substances on which they were adsorbed (POC/IM/phytoplankton).The settling rate of each carrier substance, which is the fourth term on the right hand side of Equation ( 1), is computed by Delft3D-WAQ according to Equation (3), as described by Deltares [52].
where R set is the settling rate of a carrier substance (gDM m −3 d −1 ; DM: Dry Matter), f tau is the shear stress limitation function (dimensionless), F set is the settling flux of a carrier substance (gDM m −2 d −1 or gC m −2 d −1 ), and H is the water depth (m).This equation is applied when the depth (H) is greater than the minimum water depth for sedimentation (0.1 m).Otherwise, the settling rate is equal to zero.The settling flux (F set ) depends on the concentration of the substance and is calculated according to Equations ( 4) and (5).Finally, the shear stress limitation function (f tau ) is determined by the following equation, taking into consideration the shear stress (τ), (Equation ( 6)).
where τ: the shear stress (Pa); τc: the critical shear stress for the settling of a carrier substance (Pa).

Initial Conditions and Open Boundaries
To initiate the model, a single value was used for the concentration of each trace element in the modeling area, which was retrieved from the literature [21] (Table S7).After a one-year spin-up, the spatially varying concentrations recorded on the last time-step of the simulation were used as initial conditions.Regarding the concentration of trace elements at open boundaries (Table S7), a single value retrieved from the literature was used [21], representing the mean concentration of each trace element over the decade 2000-2010 in the outer Saronikos Gulf.

Validation Data and Statistics
The performance of the hydrodynamic model (Delft3D-FLOW) was quantitatively assessed and the model's ability to reproduce the hydrodynamics of the study area is discussed in other studies [35,36].For the purposes of completeness, the hyrodynamic model validation procedure is briefly addressed in the Supplementary Materials (Figures S2-S4, along with a short discussion).Regarding water quality and pollution, the model was validated with the use of monitoring values retrieved from sampling campaigns (discrete water samples) conducted in the past (Figure 1, P-2014 data).The concentration of trace elements at different stations and over various depths during the time period of 2009-2010 were compared to the predicted values.To assess model performance quantitatively, statistical analysis was performed, and the following indexes were calculated.
The percent bias (Pbias), as the sum of model error normalized by the data, was calculated according to Equation ( 7) [55].
The Cost Function (CF) was calculated to interpret the "goodness of fit" between model results and measurements data, according to Equation ( 8) [55].
The Willmott Model Skill (MSkill) reflects the degree to which the observed variable is accurately estimated by the simulated variate and is described by Equation (9) [56].
where τ: the shear stress (Pa); τc: the critical shear stress for the settling of a carrier substance (Pa).

Initial Conditions and Open Boundaries
To initiate the model, a single value was used for the concentration of each trace element in the modeling area, which was retrieved from the literature [21] (Table S7).After a one-year spin-up, the spatially varying concentrations recorded on the last time-step of the simulation were used as initial conditions.Regarding the concentration of trace elements at open boundaries (Table S7), a single value retrieved from the literature was used [21], representing the mean concentration of each trace element over the decade 2000-2010 in the outer Saronikos Gulf.

Validation Data and Statistics
The performance of the hydrodynamic model (Delft3D-FLOW) was quantitatively assessed and the model's ability to reproduce the hydrodynamics of the study area is discussed in other studies [35,36].For the purposes of completeness, the hyrodynamic model validation procedure is briefly addressed in the Supplementary Materials (Figures S2-S4, along with a short discussion).Regarding water quality and pollution, the model was validated with the use of monitoring values retrieved from sampling campaigns (discrete water samples) conducted in the past (Figure 1, P-2014 data).The concentration of trace elements at different stations and over various depths during the time period of 2009-2010 were compared to the predicted values.To assess model performance quantitatively, statistical analysis was performed, and the following indexes were calculated.
The percent bias (Pbias), as the sum of model error normalized by the data, was calculated according to Equation ( 7) [55].
The Cost Function (CF) was calculated to interpret the "goodness of fit" between model results and measurements data, according to Equation ( 8) [55].
The Willmott Model Skill (MSkill) reflects the degree to which the observed variable is accurately estimated by the simulated variate and is described by Equation (9) [56].
where τ: the shear stress (Pa); τc: the critical shear stress for the settling of a carrier substance (Pa).

Initial Conditions and Open Boundaries
To initiate the model, a single value was used for the concentration of each trace element in the modeling area, which was retrieved from the literature [21] (Table S7).After a one-year spin-up, the spatially varying concentrations recorded on the last time-step of the simulation were used as initial conditions.Regarding the concentration of trace elements at open boundaries (Table S7), a single value retrieved from the literature was used [21], representing the mean concentration of each trace element over the decade 2000-2010 in the outer Saronikos Gulf.

Validation Data and Statistics
The performance of the hydrodynamic model (Delft3D-FLOW) was quantitatively assessed and the model's ability to reproduce the hydrodynamics of the study area is discussed in other studies [35,36].For the purposes of completeness, the hyrodynamic model validation procedure is briefly addressed in the Supplementary Materials (Figures S2-S4, along with a short discussion).Regarding water quality and pollution, the model was validated with the use of monitoring values retrieved from sampling campaigns (discrete water samples) conducted in the past (Figure 1, P-2014 data).The concentration of trace elements at different stations and over various depths during the time period of 2009-2010 were compared to the predicted values.To assess model performance quantitatively, statistical analysis was performed, and the following indexes were calculated.
The percent bias (Pbias), as the sum of model error normalized by the data, was calculated according to Equation ( 7) [55].
The Cost Function (CF) was calculated to interpret the "goodness of fit" between model results and measurements data, according to Equation ( 8) [55].
The Willmott Model Skill (MSkill) reflects the degree to which the observed variable is accurately estimated by the simulated variate and is described by Equation (9) [56].); LWAQ (d).r stress limitation function (ftau) is determined by the following onsideration the shear stress (τ), (Equation ( 6)).
s and Open Boundaries odel, a single value was used for the concentration of each trace ing area, which was retrieved from the literature [21] (Table S7).
-up, the spatially varying concentrations recorded on the last ation were used as initial conditions.Regarding the concentration open boundaries (Table S7), a single value retrieved from the ], representing the mean concentration of each trace element over in the outer Saronikos Gulf.Statistics of the hydrodynamic model (Delft3D-FLOW) was quantitatively el's ability to reproduce the hydrodynamics of the study area is dies [35,36].For the purposes of completeness, the hyrodynamic edure is briefly addressed in the Supplementary Materials (Figures ort discussion).Regarding water quality and pollution, the model e use of monitoring values retrieved from sampling campaigns s) conducted in the past (Figure 1, P-2014 data).The concentration fferent stations and over various depths during the time period of pared to the predicted values.To assess model performance ical analysis was performed, and the following indexes were (Pbias), as the sum of model error normalized by the data, was o Equation ( 7) [55].
del Skill (MSkill) reflects the degree to which the observed variable by the simulated variate and is described by Equation (9) [56].
: the mean of observations.

Numerical Error Assessment
Delft3D-WAQ employs an artificial conservative tracer called 'Continuity' that is used to check the numerical correctness and stability of the simulation.According to the model's manual, a concentration of 1 g m −3 is assigned to all water sources-initial condition, boundary conditions, and discharges-so that the 'Continuity' concentration should remain 1 g m −3 during the whole simulation, as there are no processes that dilute or concentrate it.If the concentration deviates significantly from 1 g m −3 during a model run, the simulation is numerically unstable or/and a source of the 'continuity' conservative tracer has not been considered, i.e., a 'continuity' value of 1 g m −3 has not been assigned to a water source, leading to this deviation [37].During simulations of this study, deviations from this value of 1 g m −3 were quantified: in Elefsis Bay, deviation varies between 1 and 2%, while in the rest of the area, it is below 0.5%, indicating very good model performance.

Statistical Analysis
The modeled trace elements' concentrations were compared with the available in situ values (Table S8).The total concentration of trace elements (dissolved and particulate) was used for this analysis.In total, the concentrations of three sampling campaigns (7/12/2009, 11/03/2010, and 31/08/2010) were used, over nine sampling stations (S1, S2, S3, S7, S8, S11, S13, S16, and S18) and various depths from the surface (2 m), middle, and bottom of each station (20 m, 50 m, 74 m, 84 m, 89 m).For each sampling point (x, y, z), the respective modeled concentration was determined, and, in total, 87 values were compared (modeled vs. observed concentration).Regarding Model Bias (Figure 2a), zinc had the lowest %Bias (7.58%) while lead, nickel, and copper had good performance (%Bias less than 20%).Cadmium presented the highest %Bias (27.36%), indicating that there is higher uncertainty in the prediction of Cd concentrations.The Cost function, which quantifies the 'goodness of fit' between the two datasets, indicates that all element concentrations are reproduced in a very satisfactory way (0.61 < CF < 0.98, very good), with nickel exhibiting the least good score (CF = 1.07, good).The Willmott Skill indicator was very close to that for lead, nickel, and zinc (0.87, 0.95, 0.93, respectively) and lower for cadmium and copper (0.68, 0.84).
or concentrate it.If the concentration deviates significantly from 1 g m −3 during a model run, the simulation is numerically unstable or/and a source of the 'continuity' conservative tracer has not been considered, i.e., a 'continuity' value of 1 g m −3 has not been assigned to a water source, leading to this deviation [37].During simulations of this study, deviations from this value of 1 g m −3 were quantified: in Elefsis Bay, deviation varies between 1 and 2%, while in the rest of the area, it is below 0.5%, indicating very good model performance.

Statistical Analysis
The modeled trace elements' concentrations were compared with the available in situ values (Table S8).The total concentration of trace elements (dissolved and particulate) was used for this analysis.In total, the concentrations of three sampling campaigns (7/12/2009, 11/03/2010, and 31/08/2010) were used, over nine sampling stations (S1, S2, S3, S7, S8, S11, S13, S16, and S18) and various depths from the surface (2 m), middle, and bottom of each station (20 m, 50 m, 74 m, 84 m, 89 m).For each sampling point (x, y, z), the respective modeled concentration was determined, and, in total, 87 values were compared (modeled vs. observed concentration).Regarding Model Bias (Figure 2a), zinc had the lowest %Bias (7.58%) while lead, nickel, and copper had good performance (%Bias less than 20%).Cadmium presented the highest %Bias (27.36%), indicating that there is higher uncertainty in the prediction of Cd concentrations.The Cost function, which quantifies the 'goodness of fit' between the two datasets, indicates that all element concentrations are reproduced in a very satisfactory way (0.61 < CF < 0.98, very good), with nickel exhibiting the least good score (CF = 1.07, good).The Willmott Skill indicator was very close to that for lead, nickel, and zinc (0.87, 0.95, 0.93, respectively) and lower for cadmium and copper (0.68, 0.84).The time variability of simulated Ni and Cu concentrations along with measured concentrations are plotted in Figures 3-5 for sampling stations S1, S3, S7, and S16, at depths where these field observations are available.These figures provide a visual comparison between the two series of data.The high variability of the modeled Cu concentration at the deeper layer of station S3 (Keratsini strait, Figure 3, bottom right) is attributed to the pollution source or Perama shipyard (Figure 1; no 9), in tandem with sedimentation processes acting upon the loads introduced and the hydrodynamic circulation advecting water masses and determining concentration levels.Respectively, the high variability of both Ni and Cu concentrations at the deeper layer of station S7 (Figure 4, bottom panels) is attributed to the loads discharged in the domain by the Psittalia WWTP.Conversely, at station S16, observed and modeled concentrations do not present fluctuations, which can be attributed to the fact that this station is far away from distinct pollution sources.Furthermore, station S16 can be characterized as a baseline location, with concentration levels being much lower than those observed and modeled at S1 (Elefsis Bay).

Spatial Seasonal Distribution of Trace Elements
The spatial seasonal distribution of trace elements was examined, focusing on the area of the Saronikos Gulf where the sampling stations network exists (Elefsis Bay and inner Saronikos, Figure 1b).For each trace element, concentration maps were created for the dissolved phase.The average concentrations over two periods were examined: (i) winter months (from December to April) and (ii) summer months (from May to November).This approach takes into account the pycnocline formation and circulation patterns both observed [25] and modeled [35,36], as well as similar approaches in the literature [21].At the same time, concentrations are averaged both in the upper 20% of the water column (top three layers of the sigma grid model) and the rest of the water column, dividing it into a 'surface' and a 'subsurface', deeper layer.The thickness of the surface layer may vary from approximately 6 to 7 m in the deep parts of Elefsis Bay to approximately 20 m where the bathymetry is around 100 m (see Section 3.3).These averages were adopted to compare our findings with the previously published distribution maps of trace elements for the time period of 2000-2010 [21].
concentrations are plotted in Figures 3-5 for sampling stations S1, S3, S7, and S16, at depths where these field observations are available.These figures provide a visual comparison between the two series of data.The high variability of the modeled Cu concentration at the deeper layer of station S3 (Keratsini strait, Figure 3, bottom right) is attributed to the pollution source or Perama shipyard (Figure 1; no 9), in tandem with sedimentation processes acting upon the loads introduced and the hydrodynamic circulation advecting water masses and determining concentration levels.Respectively, the high variability of both Ni and Cu concentrations at the deeper layer of station S7 (Figure 4, bottom panels) is attributed to the loads discharged in the domain by the Psittalia WWTP.Conversely, at station S16, observed and modeled concentrations do not present fluctuations, which can be attributed to the fact that this station is far away from distinct pollution sources.Furthermore, station S16 can be characterized as a baseline location, with concentration levels being much lower than those observed and modeled at S1 (Elefsis Bay).

Spatial Seasonal Distribution of Trace Elements
The spatial seasonal distribution of trace elements was examined, focusing on the area of the Saronikos Gulf where the sampling stations network exists (Elefsis Bay and inner Saronikos, Figure 1b).For each trace element, concentration maps were created for the dissolved phase.The average concentrations over two periods were examined: (i) For all trace elements, the simulated concentrations were higher in Elefsis Bay (S1, S2) and close to Psittalia (S7, S3).This is in line with the observed concentrations, and was expected, as major sources of pollution exist and were set close to those areas.On the other hand, simulated concentrations in the central inner Saronikos Gulf (close to S11) are lower than those recorded during field measurements.
Mean dissolved Cd concentrations were higher in Elefsis Bay, especially during winter (Figure 6a,b), reaching up to 0.1 µg L −1 .This can be attributed to the seasonal variation of mass loads, as water streams (Figure 1b) provide higher Cd loads during winter.In the area of inner Saronikos, mean concentrations were close to 0.02 µg L −1 .On the other hand, Cd did not present any notable difference between surface (Figure 6a,c) and subsurface (Figure 6b,d) concentrations.The uniform vertical distribution in the water column can be attributed to the very low adsorption affinity of Cd to particulate matter, e.g., [57], as the fraction of trace elements adsorbed to particulate matter is subject to sedimentation processes, thus altering the vertical distribution.In the case of Cd, the larger part of the total concentration remains in the dissolved phase.
Considering Pb, a different trend was observed (Figure 7), with higher mean concentrations being recorded close to S3 (Keratsini strait), reaching 0.25 µg L −1 .It is probable that the point source pollution close to S3 (Figure 1b), together with the small water volume of the strait, contributed to these increased concentrations.In the rest of the area (Elefsis Bay, inner Saronikos), the mean concentrations varied between 0.05 and 0.1 µg L −1 .
Regarding Ni, the mean modeled concentrations (Figure 8) were significantly higher in Elefsis Bay (close to 1 µg L −1 ) compared to the rest of the modeling area, (0.25 µg L −1 ) which is in line with the monitoring studies [21].This is due to the high contribution of pollution sources in Elefsis Bay, with more than 30% of the total mass load being discharged in that area.The actual pollution sources in Elefsis Bay can be identified on the distribution maps of Figure 8, where higher concentrations are modeled (red areas).Vertical distribution did not present important differences, with slightly higher concentrations being observed in the subsurface.
winter.In the area of inner Saronikos, mean concentrations were close to 0.02 μg L −1 .On the other hand, Cd did not present any notable difference between surface (Figure 6a,c) and subsurface (Figure 6b,d) concentrations.The uniform vertical distribution in the water column can be attributed to the very low adsorption affinity of Cd to particulate matter, e.g., [57], as the fraction of trace elements adsorbed to particulate matter is subject to sedimentation processes, thus altering the vertical distribution.In the case of Cd, the larger part of the total concentration remains in the dissolved phase.Considering Pb, a different trend was observed (Figure 7), with higher mean concentrations being recorded close to S3 (Keratsini strait), reaching 0.25 μg L −1 .It is probable that the point source pollution close to S3 (Figure 1b), together with the small water volume of the strait, contributed to these increased concentrations.In the rest of the area (Elefsis Bay, inner Saronikos), the mean concentrations varied between 0.05 and 0.1 μg L −1 .Regarding Ni, the mean modeled concentrations (Figure 8) were significantly higher in Elefsis Bay (close to 1 μg L −1 ) compared to the rest of the modeling area, (0.25 μg L −1 ) which is in line with the monitoring studies [21].This is due to the high contribution of pollution sources in Elefsis Bay, with more than 30% of the total mass load being discharged in that area.The actual pollution sources in Elefsis Bay can be identified on the distribution maps of Figure 8, where higher concentrations are modeled (red areas).Vertical distribution did not present important differences, with slightly higher concentrations being observed in the subsurface.
Regarding the mean dissolved concentrations of Cu, values were higher in the vicinity of Elefsis Bay, especially during the summer period (Figure 9c,d).Furthermore, higher concentrations were recorded in the subsurface layer (Figure 9b,d).Both these trends were in line with the findings of Paraskevopoulou et al. [21]; however, the P-2014 dataset reports a higher concentration of Cu close to S8 and S11, which could be attributed to shipping activity in the anchorage area and the shipping lane, respectively.With the current model set-up, these higher levels close to S8 and S11 (inner Saronikos)   For the case of Zn, the mean dissolved concentrations presented no significant spatial or seasonal fluctuations (Figure 10).The model managed to reproduce mean baseline concentrations varying from 3 to 5 μg L −1 , which are similar to those reported in monitoring studies [21].Regarding the mean dissolved concentrations of Cu, values were higher in the vicinity of Elefsis Bay, especially during the summer period (Figure 9c,d).Furthermore, higher concentrations were recorded in the subsurface layer (Figure 9b,d).Both these trends were in line with the findings of Paraskevopoulou et al. [21]; however, the P-2014 dataset reports a higher concentration of Cu close to S8 and S11, which could be attributed to shipping activity in the anchorage area and the shipping lane, respectively.With the current model set-up, these higher levels close to S8 and S11 (inner Saronikos) are not reproduced.For the case of Zn, the mean dissolved concentrations presented no significant spatial or seasonal fluctuations (Figure 10).The model managed to reproduce mean baseline concentrations varying from 3 to 5 μg L −1 , which are similar to those reported in monitoring studies [21].For the case of Zn, the mean dissolved concentrations presented no significant spatial or seasonal fluctuations (Figure 10).The model managed to reproduce mean baseline concentrations varying from 3 to 5 µg L −1 , which are similar to those reported in monitoring studies [21].

Vertical Distribution
To investigate the vertical distribution of trace elements, figures of cross sections were created.In Figure 11, the modeled concentration of each element (at a precise day) is presented along a cross-section starting from Elefsis Bay and passing over stations S1, S3, S7, S11, and S16 (see Figure 1).This section was chosen in order to visualize the following different areas: the highly polluted Elefsis Bay (S1), the straits connecting Elefsis Bay with the inner Saronikos (S3), the area where the Psittalia WWTP outflow is located (S7), the inner Saronikos (S11), the beginning of the outer Saronikos deeper basin (S16), and the open boundaries to the Aegean Sea (Figure 1b).To compare the modeled concentrations with observed ones (P-2014 data), dots representing the measured concentrations were added at the relevant stations and depths (Figure 11).
Overall, the model managed to reproduce the trend observed from field values.In all cases, the average daily modeled concentrations are presented for the day for which the field measurement was available.Indicatively, some of those time frames are displayed.Regarding Cd during March (Figure 11a), the concentration is higher in Elefsis Bay, where oil refineries and water streams (active in winter) are the major sources of pollution.Pb exhibits higher values (both modeled and measured) in Elefsis Bay's southeast straits (where station S3 is located) (Figure 11b).This was also reflected in distribution maps (Figure 7).Regarding Ni, during August, the concentrations were adequately reproduced in Elefsis Bay (S1), but overall, Ni presented higher field concentrations in the area of the inner Saronikos (Figure 11c).For Cu, during March (Figure 11d), the higher concentrations are observed in Elefsis Bay straits (S3), where the Ship Repairing Zone of Perama is estimated to be an important source of pollution (via antifouling paints), releasing pollutants in a small, confined volume of water.Regarding Zn, Psittalia is an important source of pollution and the effect of the Ship Repairing Zone of Perama (S3) can be seen in Figure 11e.The high measured concentration of Zn on top layers may indicate that atmospheric deposition or antifouling paints leaching from ship hulls-not examined in this study-may be an important source of pollution.

Vertical Distribution
To investigate the vertical distribution of trace elements, figures of cross sections were created.In Figure 11, the modeled concentration of each element (at a precise day) is presented along a cross-section starting from Elefsis Bay and passing over stations S1, S3, S7, S11, and S16 (see Figure 1).This section was chosen in order to visualize the following different areas: the highly polluted Elefsis Bay (S1), the straits connecting Elefsis Bay with the inner Saronikos (S3), the area where the Psittalia WWTP outflow is located (S7), the inner Saronikos (S11), the beginning of the outer Saronikos deeper basin (S16), and the open boundaries to the Aegean Sea (Figure 1b).To compare the modeled concentrations with observed ones (P-2014 data), dots representing the measured concentrations were added at the relevant stations and depths (Figure 11).
Overall, the model managed to reproduce the trend observed from field values.In all cases, the average daily modeled concentrations are presented for the day for which the field measurement was available.Indicatively, some of those time frames are displayed.Regarding Cd during March (Figure 11a), the concentration is higher in Elefsis Bay, where oil refineries and water streams (active in winter) are the major sources of pollution.Pb exhibits higher values (both modeled and measured) in Elefsis Bay's southeast straits (where station S3 is located) (Figure 11b).This was also reflected in distribution maps (Figure 7).Regarding Ni, during August, the concentrations were adequately reproduced in Elefsis Bay (S1), but overall, Ni presented higher field concentrations in the area of the inner Saronikos (Figure 11c).For Cu, during March (Figure 11d), the higher concentrations are observed in Elefsis Bay straits (S3), where the Ship Repairing Zone of Perama is estimated to be an important source of pollution (via antifouling paints), releasing pollutants in a small, confined volume of water.Regarding Zn, Psittalia is an important source of pollution and the effect of the Ship Repairing Zone of Perama (S3) can be seen in Figure 11e.The high measured concentration of Zn on top layers may indicate that atmospheric deposition or antifouling paints leaching from ship hulls-not examined in this study-may be an important source of pollution.

Partitioning
Regarding the partitioning of trace elements, the distribution between the dissolved and particulate phases was compared with the trend observed from field monitoring values (P-2014 data).As the actual monitoring values during the annual cycle Nov. 2009-Nov.2010 describe only three instances (December 2009, March, and August 2010), monitoring values from the previous years (2008) were used to examine the distribution as a percentage of the total concentration.Indicatively, the partitioning of Cu and Ni at stations S1, S7, and S16 as computed by the model are compared with the recorded partitioning at the same station and depth (Figure 12).same day are superimposed as dots at the respective stations and depths.The dashed line indicate the separation of the surface/subsurface layers.

Partitioning
Regarding the partitioning of trace elements, the distribution between the dissolved and particulate phases was compared with the trend observed from field monitoring values (P-2014 data).As the actual monitoring values during the annual cycle Nov 2009-Nov.2010 describe only three instances (December 2009, March, and August 2010) monitoring values from the previous years (2008) were used to examine the distribution as a percentage of the total concentration.Indicatively, the partitioning of Cu and Ni a stations S1, S7, and S16 as computed by the model are compared with the recorded partitioning at the same station and depth (Figure 12).Mean monthly modeled distribution of Cu and Ni, according to WAQ terminology, referring to the dissolved phase: free dissolved metal in water column (blue), adsorbed to DOC (light blue); and the particulate phase: adsorbed to POC (grey), adsorbed to phytoplankton (green), adsorbed to inorganic matter (yellow) for stations S1, S7, and S16, in comparison with the observed distribution pattern at the same month (dashed line).When monitoring values were not available (white marker), interpolation was applied between the two available distributions.On the x-axis, the months are displayed, beginning from November 2009 (N) to October 2010 (O), and on the y-axis, is the fraction of each trace element encountered in each phase.
In Figure 12, the dashed line indicates the fraction of each trace element encountered in the dissolved phase, as determined by field measurements.Similarly, in the same graph, the upper limit of the light blue bar indicates the model-predicted fraction of each trace element distributed in the dissolved phase.Ideally, the dashed line should meet the upper limit of the light blue bar, to reach a perfect agreement of measured and modeled distribution of trace elements.An example of very good agreement is the case of Ni distribution at S7 (Figure 12).The variability of Cu partitioning at station S7 during the summer months deviates from this general trend.The deviation of the model may be related to the partition coefficient (K d ) value that was attributed to Cu and is uniform for all the modeling area, but may not be representative of the area of S7 during summer.Furthermore, the selected equilibrium approach for modeling the distribution is expected to perform poorly close to pollution sources [53,54].Marine field studies highlight the difficulty in disclosing the complexity of processes related to binding mechanisms of trace elements, depending on the availability of organic ligands and the nature and origin of dissolved organic matter [58].It is possible that the different local conditions, due to WWTP outflow, such as particulate matter quality, DOC, pH, and salinity, influence the adsorption-desorption behavior of Cu [59].Another explanation may be that during summer months, a source of particulate Cu is present, and was not introduced into the model.For example, in-water hull cleaning activities of ships, anchored in the marine area southwest of S7, may lead to particle release, as in many cases the wastewater occurring from this cleaning process is not collected for treatment [60].In a similar monitoring study, investigating the partitioning of trace elements in the same marine area during the year 2011-2012, no intense seasonal differences were observed close to S7 [61].

Sedimentation
Sedimentation of trace elements to the surface sediment layer, as a consequence of the elements partitioning to the particulate phase and settling of particulate matter, was modeled.The mass of each trace element deposited to the sediment top layer (1 cm), after a one-year simulation time period, is indicated in Figure 13 (in g of trace element per m 2 ).Among the five trace elements examined, Zn attained the highest rates of sedimentation, with up to 12 g m −2 per year being deposited around the area of Psittalia WWTP outflow (Figure 13e).Furthermore, both Ni and Cu presented high deposition rates in Elefsis Bay, reaching 0.5 and 1.2 g m −2 per year, respectively, (Figure 13c,d).This is explained by the high discharge fluxes as well as the contribution of pollution sources in Elefsis Bay (approximately 30% of the total mass load of each trace element is discharged in Elefsis Bay).On the other hand, the deposition of Pb is lower (Figure 13b), due to lower mass loads entering the domain (5 times lower than Ni).Due to low sorption affinity (low K d value, Table S6; [57]) as well as low total mass loads, Cd has the lowest deposition rates, in the range of µg m −2 y −1 (Figure 13a).Among the other trace elements, Ni and Zn present lower sorption affinity compared to Pb and Cu, which cannot be directly related to the deposition rates, as the total mass loads of each trace element entering the domain are different, following the order Zn > Cu > Ni > Pb > Cd.
To compare the predicted values with the actual deposition rates, a comparison was made at stations where field data are available.In Table 3, the modeled sedimentation rates at two stations (S1, S7) are compared with the fluxes calculated based on observed concentrations of trace elements in the sediment [62] and sedimentation rates [63].Observed concentrations included measurements in nearby areas of each station (west: W; east: E; north: N).For the calculation of fluxes, the sediment mass was converted to volume (Tables S9-S11; Methodology described in Supplementary Materials).To compare the predicted values with the actual deposition rates, a comparison was made at stations where field data are available.In Table 3, the modeled sedimentation rates at two stations (S1, S7) are compared with the fluxes calculated based on observed concentrations of trace elements in the sediment [62] and sedimentation rates [63].Observed concentrations included measurements in nearby areas of each station (west: W; east: E; north: N).For the calculation of fluxes, the sediment mass was converted to volume (Tables S9-S11; Methodology described in Supplementary Materials).
Table 3. Deposition rates of trace elements at stations S1, S7, and nearby areas: simulated (columns A) and calculated mass fluxes (columns B) based on measured concentrations [62], Table S10, and reported sedimentation rates [63].Units in mg m −2 year −1 .Table 3. Deposition rates of trace elements at stations S1, S7, and nearby areas: simulated (columns A) and calculated mass fluxes (columns B) based on measured concentrations [62], Table S10, and reported sedimentation rates [63].Units in mg m −2 year −1 .The modeled rates (column A of Table 3) are within the same order of magnitude as the observed/calculated fluxes, indicating that the process of "Sedimentation" is efficiently parametrized.Between the two stations compared in Table 3, higher deposition rates are observed for Cd, Pb, Ni, and Zn close to S7.This is expected as the contribution of Psittalia WWTP as a pollution source is important in the vicinity of S7.Regarding Cu, the model seems to efficiently reproduce the sedimentation rate close to S7, but overestimates the deposition at S1, in Elefsis Bay.Similarly, the Zn deposition rate seems to be overestimated by the model.Taking into consideration that there are no exact field measurements to compare directly with our findings, this evaluation should remain at this basic level.

Discussion
The implementation of a far-field water quality model was herein described, aiming to serve as a tool for the assessment of environmental impacts on Saronikos Gulf water quality due to human activities (planned or accidental) in its coastal zone.Furthermore, the presented methodology can serve as a guide to other researchers for the application of the Delft3D-WAQ model in other case study areas.The evaluation of trace elements simulations through the comparison with field measurement data proves that the baseline concentrations were successfully reproduced.Furthermore, the predicted distribution of trace elements in the particulate and dissolved phases is efficiently simulated.Despite the different characteristics of each sub-area (Elefsis Bay, inner Saronikos, outer Saronikos), the model was able to reproduce, to a satisfactory degree, the spatiotemporal variability within the modeling area.
The detailed comparison of the predicted concentrations with the measured ones, on the exact same day and geographical area (station and depth), leads to a quantitative evaluation of the model's performance.Despite the challenges that such a 'point-to-point' comparison poses, the outcome revealed a very good reproduction of the field values in total for the year 2010.A more general comparison with the seasonal variations observed over a longer period (2000-2010) by Paraskevopoulou et al. [21] was attempted (Figures 5-9).Observations of Cd over that decade demonstrated higher concentrations close to Psittalia, while the model reproduced higher concentrations in Elefsis Bay.Regarding Pb, observations demonstrate that there is an accumulation of this trace element in the western part of Elefsis Bay, which was not reproduced by the model (Figure 6).Concentrations of Ni presented similarity, with both observed and simulated levels being higher in Elefsis Bay.On the other hand, over the decade of 2000-2010, the levels of Cu were higher in Elefsis Bay, a feature not well-reflected in the model results (Figure 8), which refer only to the year 2010.Similar monitoring studies, examining the time period 2012-2019, also report high concentrations of Cu in Elefsis Bay [64].It is highly probable that the mass loads of Cu entering the modeling area are underestimated, as various pollutant sources were not included; for example, leaching of antifouling paints from vessels' hulls, which can be an important source in the anchorage area south of Salamina island and in the marine traffic lane of the inner Saronikos Gulf.
Through this seasonal comparison, we aimed to obtain a better understanding of the model behavior in relation to the trends observed in the field.Between seasons, three factors may be influencing the differences observed in modeled concentrations: (i) the seasonal variability of pollution fluxes by land sources (rivers and streams), (ii) the seasonal circulation patterns and mixing intensity, and (iii) the distribution between dissolved and particulate phases due to sorption affinity and sorbent availability.Considering the land sources of pollution, for all trace elements, the mass loads occurring in the marine area through the user-defined riverine inputs had a high contribution to the total mass load (from October to February).It is worth mentioning that during December, more than 60% of the waste load was due to the river and water stream discharges.This indicates that it is of high importance to adequately and precisely compute the water volume discharged and the pollutant concentrations observed in freshwater sources.In our case, this estimation could be improved only if more detailed field values were available.Otherwise, the parallel use of a hydrological model predicting river heavy metal concentrations [65] could be used to obtain more accurate predictions of discharges in the marine area.
Overall, the most crucial point greatly influencing the accuracy of predictions is the determination of anthropogenic coastal pollution sources.Unfortunately, this is also the most challenging part [66], as it is almost impossible to precisely determine every humaninduced activity.Some examples of not-easily identifiable pollution sources comprise the following activities: (i) hull cleanup (in water), releasing Cu and Zn [60], (ii) leaching of metals from various parts of the ships (hull, anodes) [67,68], (iii) illegal discharges of wastewater (bilge wastewater, gray wastewater) or non-recorded accidental spills of hazardous materials [69], (iv) dredging and disposal of dredged sediment [70].Furthermore, the quantification of fluxes from pollution sources is also difficult to determine precisely [71].Pollutant fluxes from point sources (e.g., the outfall of a WWTP) can be estimated precisely only in the case where analytical measurements and recordings of discharge quality are available (e.g., daily measurement of pollutants).Indirect pollution sources also exist, such as leaching of pollutants through soil and groundwater, surface runoff, non-recorded or small-scale industrial activities, and many others.Thus, it is of high importance to successfully identify the major pollution sources, quantify their contributions, and optimize this methodology.In the current study, all available tools and information were used in that direction, and future applications will review and update this methodology to incorporate new developments and data availability (e.g., European Industrial Emissions Portal).
Similarly, natural atmospheric events in combination with the coastal and maritime air emissions may contribute seasonally or occasionally to the mass loads directly deposited to the surface layer [72].For example, Remoundaki et al. [73] reported atmospheric Zn concentrations that were 15 times higher than usual during a Saharan dust transport event in 2009 in Athens.Such phenomena may locally influence the concentration of trace elements observed.In the current study, atmospheric loads were not computed but it is within the scope of future research.
Modeling tools have been used for the determination of the fate of heavy metals in the aquatic environment in studies similar to ours.The simulation of the hydrodynamic and biogeochemical processes in similar case studies and with the use of the same software (Delft3D) is described by Vaz et al. [74] and Mendes et al. [75].Regarding pollution simulation, in some cases, models are used to identify pollution sources through back-tracking [76].Most applications focus on riverine environments, where high levels of pollution are identified [77,78].In this case, the flow velocity, freshwater characteristics, and sediment behavior dominate the distribution mechanism of pollutants.In other studies, models are applied to investigate the distribution of heavy metals in estuarine waters, where tidal cycles, salinity, and suspended sediments are the main drivers [79][80][81][82].Our study could better be compared with others, focusing on sea areas where salinity levels are relatively stable and suspended particulate matter is of low concentration.Such a case is the study of Periáñez [83], investigating the distribution of Cu, Ni, and Zn in the Gulf of Cadiz.In this study, sediment transport and tidal cycles are examined, while riverine input of metals is dealt with as a steady state condition, and the background concentrations of examined metals are successfully reproduced.In our study, we have estimated the fluctuation of freshwater inputs due to rainfall events, although we had to use an average concentration of pollutants in river water.This value could be improved if more real-time measurements were available and under different conditions (after or during a rain event).Similarly, Fang et al. [84] have modeled the distribution of Cu, Cd, Zn, Pb, and Ni in Hangzhou Bay (China).In that study, the role of the sediment layer is investigated, and the environmental impact of a nuclear power plant located onshore is discussed.In a similar approach, Ma et al. [85] attempt to predict the hydraulic residence time of pollutants occurring from point and nonpoint pollution sources in Xiamen Bay (China).In almost all of the abovementioned studies, the difficulties encountered and the uncertainty of predictions, related to various factors, are highlighted.
In all the above-mentioned cases, the models were developed to serve as tools assisting the investigation of marine processes in areas presenting ecological interest or facing intense pressure.The current study contributes in the same direction, providing one more tool for better management and protection of the marine environment of the Saronikos Gulf.With further development of such tools, and coupled with monitoring activities, better results in water quality surveillance could be achieved [86,87].

Conclusions
The described model set-up successfully reproduces the concentration levels of five trace elements (Cd, Pb, Ni, Cu, and Zn) in the marine coastal area of the Saronikos Gulf, Eastern Mediterranean.Major land-based pollution sources were identified, and pollutant mass loads were estimated and used as inputs to force simulations.The partitioning of the examined metals between the dissolved and particulate phases was efficiently simulated.The modeled accumulation rate of pollutants in the sediment layer is within the observed accumulation rates.
To improve model performance, it is essential to develop methodologies to quantify certain pollution fluxes characterized by a large degree of uncertainty, thus not included in the present study.This task was beyond the scope of this application.For example, additional pollution sources may include fluxes from industries such as 'Skaramangas shipyards', atmospheric deposition, shipping discharges (e.g., grey, bilge, and scrubber water), and antifouling paints.Concurrently, more marine ecosystem processes could be modeled (e.g., simulation of sediment layer interactions and resuspension), extending simulation capabilities.Lastly, to further develop and optimize the model, it is necessary to obtain larger datasets of field measurements, covering more areas of the gulf and more recent time periods.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jmse12050797/s1,Table S1: Sampling Stations network information; Table S2: Analytical table of processes activated in Delft3D-WAQ module for the simulation of ecosystem dynamics; Table S3: Table of processes values, activated in Delft3D-WAQ module, for the simulation of ecosystem dynamics; Table S4: Nutrient loads introduced in the modeling area as discharges; Table S5: Inorganic Matter (IM) average monthly concentration used in the simulations; Figure S1: Perama small-scale industries ship repairing zone, (a) google earth capture, and (b) marine traffic capture; Table S6: Kd time-series table-seasonal values, calculated according to monitoring data [21] and average SS concentrations [42]; Table S7 S8: Statistical analysis of observed vs. modeled concentrations of trace elements; Table S9: Predicted mass of heavy metals reaching the top sediment layer after 365 days of simulation (in mg m−2 year−1); Table S10: Measured concentrations of heavy metals in the top sediment layer (in mg kg −1 ), year 2010, as reported by Karageorgis et al., 2020 [62]; Table S11: Mass fluxes of trace elements (in mg m −2 year −1 ), based on measured concentrations [62] (Table S10) and reported sedimentation rates [63].
where C x : the concentration of a carrier substance (gDM m −3 or gC m −3 ); F set0 : the zero-order settling flux of a carrier substance (gDM m −2 d −1 or gC m −2 d −1 ); s: the settling velocity of a carrier substance (m d −1 ); ∆t: the timestep in DELWAQ (d).
i : each observation (concentration, field monitoring value); M i : modeled concentration; n: the number observations/modeled concentrations; σ O : the standard deviation of observations; f a carrier substance (gDM m −3 or gC m −3 ); tling flux of a carrier substance (gDM m −2 d −1 or gC m −2 d −1 ); of a carrier substance (m d −1

Figure 5 .
Figure 5. Simulated (line) and measured (dots, P-2014 data) concentrations of nickel and copper at station S16 (outer Saronikos Gulf) and various depths.

Figure 5 .
Figure 5. Simulated (line) and measured (dots, P-2014 data) concentrations of nickel and copper at station S16 (outer Saronikos Gulf) and various depths.

J 27 Figure 11 .Figure 11 .
Figure 11.Daily average vertical distribution of modeled (a) cadmium, (b) lead, (c) nickel, (d) copper, and (e) zinc concentration along the section shown in Figure 1 as a red line (Elefsis Bay-Inner Saronikos).Concentrations of trace elements measured in the field (P-2014 data) on the Figure 11.Daily average vertical distribution of modeled (a) cadmium, (b) lead, (c) nickel, (d) copper, and (e) zinc concentration along the section shown in Figure 1 as a red line (Elefsis Bay-Inner Saronikos).Concentrations of trace elements measured in the field (P-2014 data) on the same day are superimposed as dots at the respective stations and depths.The dashed line indicates the separation of the surface/subsurface layers.

Figure 12 .
Figure 12.Mean monthly modeled distribution of Cu and Ni, according to WAQ terminology referring to the dissolved phase: free dissolved metal in water column (blue), adsorbed to DOC Figure 12.Mean monthly modeled distribution of Cu and Ni, according to WAQ terminology, referring to the dissolved phase: free dissolved metal in water column (blue), adsorbed to DOC (light blue); and the particulate phase: adsorbed to POC (grey), adsorbed to phytoplankton (green), adsorbed to inorganic matter (yellow) for stations S1, S7, and S16, in comparison with the observed distribution pattern at the same month (dashed line).When monitoring values were not available (white marker), interpolation was applied between the two available distributions.On the x-axis, the months are displayed, beginning from November 2009 (N) to October 2010 (O), and on the y-axis, is the fraction of each trace element encountered in each phase.

:
Initial concentrations of trace elements used for the start-up of the model and concentrations at open boundaries; Figure S2.Comparison of Sea Surface Temperature (SST, • C, left) and SST anomaly (monthly average subtracted, right) derived from Delft3D-FLOW simulation (red line) and Copernicus satellite data (blue line), averaged over the inner Saronikos Gulf basin.The Root Mean Square Error (RMSE) and the correlation coefficient (r) in each case are reported.;Figure S3.Simulated (line) and measured with CTD (dots) sea water temperature (oC) at various depths of station S16 (inner Saronikos Gulf) during the study period (November 2009-October 2010).Simulation mean error = −0.44 • C, statistically significant correlation coefficient r = 0.93 (p-value < 0.001).; Figure S4.Simulated (line) and measured with CTD (dots) salinity (psu) at various depths of station S16 (inner Saronikos Gulf) during the study period (November 2009-October 2010).Simulation mean error = 0.26 psu, statistically significant correlation coefficient r = 0.61 (p-value < 0.001).;Table

Table 1 .
Table of processes activated in the Delft3D-WAQ module for the simulation of ecosystem dynamics.

Table 2 .
Pollution loads introduced in the modeling area as discharges.