HR3DHG version 1: modeling the spatiotemporal dynamics of mercury in the Augusta Bay (southern Italy)

Abstract

The biogeochemical dynamics of $\mathrm{Hg}$, and specifically of its three species ${\mathrm{Hg}}^{\mathrm{0}}$, ${\mathrm{Hg}}^{\mathrm{II}}$, and $\mathrm{MeHg}$ (elemental, inorganic, and organic, respectively), in the marine coastal area of Augusta Bay (southern Italy) have been explored by the high-resolution 3D Hg (HR3DHG) model, namely an advection–diffusion–reaction model for dissolved mercury in the seawater compartment coupled with a diffusion–reaction model for dissolved mercury in the pore water of sediments in which the desorption process for the sediment total mercury is taken into account. The spatiotemporal variability of the mercury concentration in both seawater ($\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$) and the first layers of bottom sediments ($\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ and$\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$), as well as the $\mathrm{Hg}$ fluxes at the boundaries of the 3D model domain, have been theoretically reproduced, showing acceptable agreement with the experimental data collected in multiple field observations during six different oceanographic cruises. Also, the spatiotemporal dynamics of the total mercury concentration in seawater have been obtained by using both model results and field observations. The mass balance of the total $\mathrm{Hg}$ species in seawater has been calculated for the Augusta Harbour, improving previous estimations. The HR3DHG model could be used as an effective tool to predict the spatiotemporal distributions of dissolved and total mercury concentrations, while contributing to better assessing hazards for the environment and therefore for human health in highly polluted areas.

1 Introduction

The investigation of the biogeochemical dynamics of $\mathrm{Hg}$ species in the marine environment addresses the need to accurately model sources and pathways of this priority contaminant within and among the different abiotic and biotic compartments of the aquatic ecosystem . Over the last few years some theoretical studies have offered innovative tools to reproduce the mass balance and the dynamics of $\left[\mathrm{Hg}\right]$ in the marine environment by means of biogeochemical models based on interconnected zero-dimensional boxes representing water or sediment compartments: among these are the River MERLIN-Expo model  and the WASP (Water Analysis Simulation Program) model . In particular, the River MERLIN-Expo model  has been used to reproduce the spatiotemporal distribution of inorganic and organic contaminants in the 1D domain of rivers and to calculate the mass balance for each of them. Although the model is able to describe many of the physical and chemical processes involved in fresh water and sediment, it specifically targets environments characterized by (i) nearly homogeneous water bodies and (ii) limited variations in landscape geometry. The WASP models have been used to simulate the $\mathrm{Hg}$ cycle within aquatic ecosystems characterized by well-mixed water layers and homogeneous sediment layers coupled through the boundary conditions at the water–sediment interface . In particular, a WASP model applied to a 1D domain and calibrated by using experimental data for dissolved $\mathrm{Hg}$ and $\mathrm{MeHg}$ allowed us to explore $\left[\mathrm{Hg}\right]$ dynamics in the Black Sea . Similarly, the box model approach has been adopted in a 2D configuration  to calculate $\mathrm{Hg}$ mass balance in the coastal areas of the Marano–Grado lagoon (northern Italy), where heterogeneous spatial distributions of $\mathrm{Hg}$ species have been observed experimentally. In general, models based on zero-dimensional boxes do not deliver reliable concentration values of contaminants in highly heterogeneous environments unless they provide high spatial resolution and a proper parameterization of the biogeochemical system.

For these reasons, in a recent work  the biochemistry of $\mathrm{Hg}$ in aquatic ecosystems has been studied using a high-resolution (HR) 1D advection–reaction–diffusion model, in which a mercury module has been integrated with the Bottom RedOx Model (BROM)  to reproduce the vertical dynamics of the total dissolved $\mathrm{Hg}$ and $\mathrm{MeHg}$ in the marine coastal areas of the Étang de Berre lagoon (France) . However, even this model includes some criticalities in the estimation of mercury dynamics. For example, the temporal variations of mercury benthic fluxes, due to reaction and diffusion processes that involve mercury species present in sediments, are not taken into account in the boundary conditions of this model. On the other hand, sediment chemistry and diffusion were investigated recently by , who devised a high-resolution 1D model for $\mathrm{Hg}$ species present in water and sediments of the Baltic Sea . In both HR models, however, the strong impact of the horizontal velocity field on the spatiotemporal distribution of $\left[\mathrm{Hg}\right]$ could not be considered since 1D modeling was used.

In general, the appropriate modeling to reproduce the spatial and temporal variability of $\mathrm{Hg}$ species in highly heterogeneous marine ecosystems, such as Augusta Harbour, requires the use of a hydrodynamics model integrated with a biogeochemical model . For this aim, introduced a PCFLOW3D model upgraded with the biogeochemical module for simultaneously simulating the velocity field of marine currents, suspended particle transport, and mercury biogeochemical transformations for the whole Mediterranean Sea. The modified PCFLOW3D is a nonstationary 3D model, which consists of four real-time integrated modules: (i) a hydrodynamic module and (ii) transport-dispersion module, both based on the finite-volume method and implemented for obtaining the velocity field of marine currents and the turbulent diffusivities; (iii) a sediment–transport module used to simulate the transport, sedimentation, and resuspension of solid particles; and (iv) a biogeochemical module to reproduce the advection, diffusion, and reaction processes of $\mathrm{Hg}$ species. Although the grid used did not guarantee a high spatial resolution, the modified PCFLOW3D model allowed for obtaining, for all the $\mathrm{Hg}$ species, theoretical vertical profiles of $\left[\mathrm{Hg}\right]$ in acceptable agreement with experimental data for most of the Mediterranean Sea  and to improve the flux estimations of $\mathrm{Hg}$ mass balance for the whole Mediterranean basin . By following the same modeling approach of , over the last decade several authors have used 3D advection–diffusion–reaction models to simulate the spatiotemporal dynamics of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ in oceans, lakes, and estuaries . However, the $\mathrm{Hg}$ partition mechanisms between the liquid phase and (biotic and abiotic) particulate organic matter (POM) were explicitly included in only a few studies. Among these, reproduced the $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ in oceans and calculated $\mathrm{Hg}$ mass balance by using a 3D ocean tracer model (OFFTRAC) coupled with a general circulation model (GEOS-Chem) . Here, the sinking flux of $\mathrm{Hg}$ bound to POM was calculated by exploiting remote sensing data for net primary production (NPP) and chlorophyll concentration, which are associated with phytoplankton abundance.

All these approaches do not allow for a fine representation of the spatial variability by approximating the model domain as a set of interconnected boxes or by detailing only in the seawater compartment the spatiotemporal dynamics of the investigated chemical species. For these reasons, we developed a new model to reproduce the spatiotemporal dynamics of $\left[\mathrm{Hg}\right]$ in polluted marine sites characterized by very high spatial heterogeneity, such as the Augusta Harbour. In the present work we report results obtained using a 3D advection–diffusion–reaction biogeochemical model for three $\mathrm{Hg}$ species in seawater (${\mathrm{Hg}}^{\mathrm{0}}$, ${\mathrm{Hg}}^{\mathrm{II}}$, and $\mathrm{MeHg}$) coupled with a diffusion–reaction model for dissolved $\mathrm{Hg}$ in the pore water of sediments. The model, named HR3DHG, has been applied to the investigation of mercury dynamics in Augusta Bay (southern Italy; see Fig. 1), specifically in its harbor, a highly polluted coastal site. In this area, a substantial experimental dataset has been collected and improved upon in recent years : oceanographic cruises and data on key physical and chemical parameters from the atmosphere, seawater, and sediments are used to verify and validate the modules of HR3DHG for a reliable and accurate high-resolution investigation of the spatiotemporal dynamics of $\mathrm{Hg}$ at highly contaminated coastal marine sites. The HR3DHG model has been designed to predict the biogeochemical behavior of $\mathrm{Hg}$ in seawater and sediments, specifically in confined and highly polluted marine coastal areas. It offers the opportunity to explore the effects of both the desorption dynamics of total mercury (${\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}$) in sediments and of ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}$ diffusion dynamics in pore water in nearly steady conditions. To this aim, in the model we consider both the sediment–pore water distribution coefficients and the desorption rate for the total mercury concentration in the sediment. The former described the ratio between adsorption and desorption rate constants at the steady state without considering perturbations induced by mercury concentration reduction in pore water. The latter reproduced the effects of these perturbations on the solid phase of the sediments.

Moreover, the role played by the spatiotemporal behavior of phytoplankton  and the mechanisms responsible for the uptake of $\mathrm{Hg}$ within cells  are taken into account as specific contributions to the scavenging process and $\mathrm{Hg}$ release process by POM, respectively. Also, seasonal oscillations of key environmental variables (velocity of marine currents, amount of precipitation, elemental and inorganic mercury concentration in the atmosphere, etc.) are taken into account.

The main objectives of the HR3DHG model can be synthesized as follows: (i) to accurately reproduce and localize the peaks of $\left[\mathrm{Hg}\right]$ within the 3D domain; (ii) to estimate the $\mathrm{Hg}$ fluxes at domain boundaries; and (iii) to predict the evolution of mercury in the sediment of polluted sites. Moreover, the HR3DHG model offers the possibility to describe the $\mathrm{MeHg}$ and ${\mathrm{Hg}}^{\mathrm{II}}$ partition between the dissolved phase (both seawater and pore water) and the particulate phase (suspended particulate matter and sediment particles). Specifically, in the dissolved phase the model describes the overall behavior of $\mathrm{Hg}$ in ionic form and complexed with dissolved organic carbon (DOC). Finally, the HR3DHG model can be a useful tool to predict and prevent the risks for human health in marine areas close to industrial sites affected by $\mathrm{Hg}$ pollution extended for very long time intervals.

The paper is organized as follows: a brief overview of the study site is provided in Sect. 2. The HR3DHG model and the model simulation setup are described in Sect. 3, referring to the Supplement for further details. In Sect. 4 the obtained results are reported and compared with experimental data. In Sect. 5 the model and the results are discussed, and, finally, conclusions are drawn in Sect. 6.

2 The study area

The Augusta Bay (Fig. 1) is a semi-closed marine area that occupies a surface of about $\mathrm{30}$ km${}^{\mathrm{2}}$ on the eastern coast of Sicily (southern Italy). The location of one of the most important harbors of the Mediterranean over time since the early 1960s, the Augusta site also hosts several industrial plants, which have adversely affected the whole area with the diffusion of several priority pollutants. In particular, a huge amount of $\mathrm{Hg}$ from one of the largest European chlor-alkali plants (Syndial Priolo Gargallo) was discharged into the sea without any treatment until the 1970s, when waste treatment became operational . Although discharge activities were definitively stopped in 2005, the $\mathrm{Hg}$ contamination from the chlor-alkali plant remains a critical environmental threat, with extremely high $\left[\mathrm{Hg}\right]$ in the bottom sediments , significant $\mathrm{Hg}$ evasion fluxes from sediments to seawater and to the atmosphere , and evident and recently documented risks for the ecosystem  and for human health . The geographical position, together with its geological and oceanographic features, assigns this area a key role in the $\mathrm{Hg}$ inventory at Mediterranean scale. The estimate of the Hg export from Augusta Bay to the open sea (0.54 kmol yr${}^{-\mathrm{1}}$; ) corresponds to about 4 % of total input from coastal point and diffuse sources to the Mediterranean Sea (12.5 kmol yr${}^{-\mathrm{1}}$; ). A very narrow shelf develops down to $\mathrm{100}$–$\mathrm{130}$ m with a mean gradient of about $\mathrm{1.0}$${}^{\circ }$ and a next steep slope characterized by a dense net of canyons dropping to the deep Ionian basin . The Augusta Harbour covers a surface of $\mathrm{23.5}$ km${}^{\mathrm{2}}$ with two main inlets connecting with the open sea: the Scirocco ($\mathrm{300}$ m wide and $\mathrm{13}$ m deep) and the Levante inlets ($\mathrm{400}$ m wide and $\mathrm{40}$ m deep). The bottom is mainly flat with an average depth of $\mathrm{15}$ m, with the exception of a deeper channel about $\mathrm{30}$ m deep connecting the inner part of the harbor with the Levante inlet. Water circulation inside the port and the exchanges through the inlets are mainly ruled by the wind and tidal forcing. Tidal fluctuations are generally low, with amplitudes ranging between $\mathrm{10}$ and $\mathrm{20}$ cm, and the winds are generally from the northwest and northeast with an average speed of around $\mathrm{3}$ m s${}^{-\mathrm{1}}$ . Water circulation in the outer coastal areas is also mainly affected by wind and tidal forcing and only weakly influenced by the outer baroclinic ocean circulation, which takes place mainly from the shelf break area offshore.

Figure 1

Map of the area under investigation including the Augusta Bay and the eponymous harbor. The sampling sites of each oceanographic survey are indicated with different symbols.

[Figure omitted. See PDF]

3 Model description

The HR3DHG model has been designed and implemented to reconstruct, at high spatiotemporal resolution, the behavior of $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$. The model consists of an advection–diffusion–reaction model for the seawater compartment coupled with a diffusion–reaction sub-model for pore water, in which the dynamics of the desorption of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ between the solid (sediments) and liquid phase (pore water) are considered.

As in the PCFLOW3D model of , the module of the biogeochemical model for the seawater compartment is integrated with a hydrodynamics module (see Fig. 2). Specifically, SHYFEM (Shallow-water HYdrodynamic Finite-Element Model) is used to calculate the spatiotemporal behavior of the horizontal components of the velocity field in the seawater compartment , fixing to zero the vertical velocity according to the experimental data (see Sect. S3 of the Supplement for details). In the HR3DHG model, the mercury exchange between the abiotic and biotic compartments is also taken into account. For this purpose, the spatiotemporal behavior of picoeukaryote abundance is reproduced by the Nutrient–Phytoplankton (NP) model  (see Sect. S4 of the Supplement for details). By using the curve of the mean vertical profile obtained by the picoeukaryote abundances are converted into the chlorophyll concentration, which allows for the reproduction of the spatiotemporal distribution of NPP. This is used in our model to calculate both the biological rate constants and the sinking flux of $\mathrm{Hg}$ adsorbed by POM. The amount of Hg absorbed and released by each picoeukaryote cell in seawater is calculated by using the Phytoplankton MERLIN-Expo model  (see Sect. S5 of the Supplement for details). The two modules are coupled with the advection–diffusion–reaction sub-model in order to reproduce the spatiotemporal behavior of the load of dissolved $\mathrm{Hg}$ released by dead picoeukaryote cells in the seawater compartment (see Fig. 2).

Figure 2

Basic structure of the HR3DHG model.

[Figure omitted. See PDF]

Figure 3

Basic scheme used for implementing the HR3DHG model. The scheme describes the transformation processes ($k_{\mathrm{1}}$ – photooxidation, $k_{\mathrm{2}}$ – photoreduction, $k_{\mathrm{3}}$ – biological oxidation, $k_{\mathrm{4}}$ – biological reduction, $k_{\text{Ph-de}}$ – photodemethylation, $k_{\text{deme}}$ – demethylation, $k_{\text{me}}$ – methylation, $K_{\text{demeth}}$ – demethylation, $K_{\text{meth}}$ – methylation) and the main transport processes ($A$ – anthropogenic input, AD – atmospheric deposition, $B$ – benthic flux, $D$ – net flux due to particulate deposition and settling, $O$ – net outflow from basin, $V$ – atmospheric evasion), which involve dissolved and particulate-bound $\mathrm{Hg}$ species in seawater (${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{II}}$, ${\mathrm{MeHg}}_{\mathrm{D}}$, ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{0}}$, ${\mathrm{Hg}}_{\mathrm{P}}^{\mathrm{II}}$, ${\mathrm{MeHg}}_{\mathrm{P}}$) and sediments (${\mathrm{Hg}}_{\mathrm{pw}}^{\mathrm{II}}$, ${\mathrm{MeHg}}_{\mathrm{pw}}$, ${\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}$).

[Figure omitted. See PDF]

3.1

The advection–diffusion–reaction model for the $\mathrm{Hg}$ species in seawater

The dynamics of the $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ in the Augusta Bay have been reproduced using an advection–diffusion–reaction model. Specifically, the model equations are solved to identify the behavior of the three main $\mathrm{Hg}$ species in seawater, indicated by ${\mathrm{Hg}}^{\mathrm{0}}(x,y,z,t)$, ${\mathrm{Hg}}^{\mathrm{II}}(x,y,z,t)$, and $\mathrm{MeHg}(x,y,z,t)$, which denote the concentrations of each $\mathrm{Hg}$ species in the position ($x$,$y$,$z$) within the three-dimensional domain at a specific time $t$ and whose reciprocal interactions are modeled with the reaction terms of the partial differential equations (PDEs). Since the experimental data indicate that the ${\mathrm{Me}}_{\mathrm{2}}\mathrm{Hg}$ concentration is very low in the Augusta Harbour , the behavior of this $\mathrm{Hg}$ species is not reproduced in our model. The spatial domain is composed of the sum of several sub-domains (regular parallelepipeds), which cover the bathymetric map of the Augusta Bay . Specifically, $z$ represents the depth of the barycenter of each sub-domain, localized between the surface ($z=\mathrm{0}$) and the bottom ($z=z_{\mathrm{b}}$), while $x$ and $y$ indicate the distance in meters measured from a reference point (lat. $\mathrm{37}{}^{\circ }\mathrm{14.618}$' N, long. $\mathrm{15}{}^{\circ }{\mathrm{11.069}}^{\prime }$ E) located northwest of the town of Augusta.

The model for both compartments is coded in C++ and adopts a finite-volume scheme in explicit form with spatial and temporal discretizations treated separately. The approach followed allows for the combination of various types of discretization procedures for solving the diffusion, advection, and reaction terms. Specifically, the differential equations are solved by performing centered-in-space differencing for the diffusion terms and first-order upwind-biased differencing for the advection terms.

The model domain in seawater is constituted by a mesh of $\mathrm{10}$ and $\mathrm{18}$ elements regularly spaced at $\mathrm{454.6}$ m in both the $x$ and $y$ direction and with a variable number of vertical layers of 5 m depth in the $z$ direction. The mesh covers the whole Augusta Harbour and part of the adjacent coastal area. In Fig. 1 the model domain is shown along with the location of the open boundaries corresponding to the two port inlets. In both compartments (seawater and sediment), a fixed time step of 300 s has been chosen to satisfy the stability conditions and constraints associated with the numerical method adopted . The stability analysis, performed according to previously published methods , indicates that the convergence of our algorithm is guaranteed.

The dynamics of the $\mathrm{Hg}$ species in seawater are represented through six processes : (i) photochemical and biological redox transformations (reaction terms); (ii) methylation–demethylation reactions (reaction terms); (iii) movement due to turbulence (diffusion terms); (iv) passive drift due to marine currents (advection terms); (v) organic and inorganic particle scavenging; and (vi) organic particle remineralization.

The photochemical and biological redox transformations between ${\mathrm{Hg}}^{\mathrm{0}}$ and ${\mathrm{Hg}}^{\mathrm{II}}$ have been described as reaction terms with a first-order kinetic . In particular, the rate constants of photochemical redox reactions are directly proportional to the shortwave radiation flux at the sea surface attenuated along the water column due to dissolved organic carbon (DOC) and suspended particulate matter (SPM) . At the same time, the rate constants of biological redox reactions are proportional to the organic carbon remineralization rate (OCRR), which depends on the net primary production at the sea surface (NPP), surface chlorophyll concentration, and surface atmospheric temperature .

The model includes three reaction terms regulated by first-order kinetics, which describe the photodemethylation of $\mathrm{MeHg}$, the methylation of ${\mathrm{Hg}}^{\mathrm{II}}$, and the biotic demethylation of $\mathrm{MeHg}$ . The first is the amount of ${\mathrm{Hg}}^{\mathrm{0}}$ produced by $\mathrm{MeHg}$ through photochemical reactions. The second is the amount of $\mathrm{MeHg}$ obtained by ${\mathrm{Hg}}^{\mathrm{II}}$ through biotic and abiotic pathways in seawater. The third is the amount of ${\mathrm{Hg}}^{\mathrm{II}}$ produced by $\mathrm{MeHg}$ through reductive demethylation processes caused by the activity of bacteria in contaminated environments. The rate constants of the three reaction terms are fixed according to previous works .

The PDEs include terms of advection and diffusion for each dimension of the 3D domain. In particular, the diffusion terms reproduce the effects of turbulence on the 3D distribution of ${\mathrm{Hg}}_{\mathrm{D}}$ through horizontal ($D_x$ and $D_y$) and vertical ($D_z$) turbulent diffusivities, which are fixed as constant (see Sect. S1 of the Supplement).

The advection terms describe the effects on the $\mathrm{Hg}$ distributions induced by (i) the horizontal velocity components ($v_x(x,y,z,t)$ and $v_y(x,y,z,t)$) of the marine currents along the $x$ and $y$ directions and (ii) the vertical velocity component ($v_z(x,y,z)$) along the $z$ direction. The horizontal velocities are calculated using results achieved by applying a hydrodynamic model to the area  (see Sect. S3 of the Supplement) and change as a function of space and time. The vertical velocity is fixed to zero according to available experimental data.

Moreover, we estimated the dynamics of the dissolved ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ species, also considering effects due to (i) adsorption by SPM (scavenging process) and (ii) release by particulate organic matter. The scavenging process for both ${\mathrm{Hg}}_{\mathrm{D}}$ species is regulated by the sinking flux of particle-bound mercury along the water column , which depends on variables calculated by using the NP model. The amount of $\mathrm{Hg}$ released by particulate organic matter is primarily estimated through variables defined in the NP model and Phytoplankton MERLIN-Expo model  (see Sects. S4 and S5 of the Supplement). Specifically, the NP model provides the spatiotemporal distribution of picoeukaryote abundance, which is used to obtain the chlorophyll concentration and the net primary production through suitable conversion functions  (see Sects. S1 and S4 of the Supplement). These two variables are then exploited to calculate the contribution of the sinking flux for POM-bound $\mathrm{Hg}$ within the suspended particulate matter. The Phytoplankton MERLIN-Expo model gives the spatiotemporal dynamics of the ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ contents within the picoeukaryote cells . These two variables are then used, together with the picoeukaryote abundance, to obtain the amount of ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ released by the dead picoeukaryote cells (see Sect. S1.2 and S1.3 of the Supplement).

Thus, the advection–diffusion–reaction model for the $\mathrm{Hg}$ species in seawater is defined by the following coupled partial differential equations. $$\begin{array}{}\text{1}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial t}}={\displaystyle \frac{\partial }{\partial x}}\left[D_x{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial x}}\right]-{\displaystyle \frac{\partial }{\partial x}}\left(v_x{\mathrm{Hg}}^{\mathrm{0}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left[D_y{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial y}}\right]\\ & -{\displaystyle \frac{\partial }{\partial y}}\left(v_y{\mathrm{Hg}}^{\mathrm{0}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial z}}\right]-{\displaystyle \frac{\partial }{\partial z}}\left(v_z{\mathrm{Hg}}^{\mathrm{0}}\right)\\ & +k_{\text{Ph-de}}\cdot \mathrm{MeHg}-(k_{\mathrm{1}}+k_{\mathrm{3}})\cdot {\mathrm{Hg}}^{\mathrm{0}}+(k_{\mathrm{2}}+k_{\mathrm{4}})\\ & \cdot {\mathrm{Hg}}^{\mathrm{II}}+S_{\mathrm{L}}^{\mathrm{0}}\end{array}\text{2}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial t}}=+{\displaystyle \frac{\partial }{\partial x}}\left[D_x{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial x}}\right]-{\displaystyle \frac{\partial }{\partial x}}\left(v_x{\mathrm{Hg}}^{\mathrm{II}}\right)\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[D_y{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial y}}\right]-{\displaystyle \frac{\partial }{\partial y}}\left(v_y{\mathrm{Hg}}^{\mathrm{II}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial z}}\right]\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(v_z{\mathrm{Hg}}^{\mathrm{II}}\right)+(k_{\mathrm{1}}+k_{\mathrm{3}})\cdot {\mathrm{Hg}}^{\mathrm{0}}-(k_{\mathrm{2}}+k_{\mathrm{4}})\cdot {\mathrm{Hg}}^{\mathrm{II}}\\ & -k_{\text{me}}\cdot {\mathrm{Hg}}^{\mathrm{II}}+k_{\text{deme}}\cdot \mathrm{MeHg}+S_{\mathrm{L}}^{\mathrm{II}}+S_{\mathrm{DOM}}^{\mathrm{II}}-S_{\mathrm{SPM}}^{\mathrm{II}}\end{array}\text{3}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \mathrm{MeHg}}{\partial t}}=+{\displaystyle \frac{\partial }{\partial x}}\left[D_x{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial x}}\right]-{\displaystyle \frac{\partial }{\partial x}}\left(v_x\mathrm{MeHg}\right)\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[D_y{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial y}}\right]-{\displaystyle \frac{\partial }{\partial y}}\left(v_y\mathrm{MeHg}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[D_z{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial z}}\right]-{\displaystyle \frac{\partial }{\partial z}}\left(v_z\mathrm{MeHg}\right)-K_{\text{Ph-de}}\\ & \cdot \mathrm{MeHg}+k_{\text{me}}\cdot {\mathrm{Hg}}^{\mathrm{II}}-k_{\text{deme}}\cdot \mathrm{MeHg}+S_{\mathrm{L}}^{\mathrm{MM}}\\ & +S_{\mathrm{DOM}}^{\mathrm{MM}}-S_{\mathrm{SPM}}^{\mathrm{MM}}\end{array}\end{array}$$ Here, $k_{\mathrm{1}}$, $k_{\mathrm{2}}$, $k_{\mathrm{3}}$, and $k_{\mathrm{4}}$ are the rate constants for the photooxidation of ${\mathrm{Hg}}^{\mathrm{0}}$, the photoreduction of ${\mathrm{Hg}}^{\mathrm{II}}$, the biological oxidation of ${\mathrm{Hg}}^{\mathrm{0}}$, and the biological reduction of ${\mathrm{Hg}}^{\mathrm{II}}$, respectively (h${}^{-\mathrm{1}}$); $k_{\text{Ph-de}}$ is the rate constant for the photodemethylation of $\mathrm{MeHg}$ (h${}^{-\mathrm{1}}$); $k_{\text{deme}}$ and $k_{\text{me}}$ are the rate constants for the biotic demethylation of $\mathrm{MeHg}$ and the methylation of ${\mathrm{Hg}}^{\mathrm{II}}$, respectively (h${}^{-\mathrm{1}}$); $S_{\mathrm{L}}^{\mathrm{0}}$, $S_{\mathrm{L}}^{\mathrm{II}}$, and $S_{\mathrm{L}}^{\mathrm{MM}}$ are the direct loads for ${\mathrm{Hg}}^{\mathrm{0}}$, ${\mathrm{Hg}}^{\mathrm{II}}$, and $\mathrm{MeHg}$, respectively ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}$); $S_{\mathrm{DOM}}^{\mathrm{II}}$ and $S_{\mathrm{DOM}}^{\mathrm{MM}}$ are the loads of ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{II}}$ and ${\mathrm{MeHg}}_{\mathrm{D}}$, respectively, released by POM ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}$); $S_{\mathrm{SPM}}^{\mathrm{II}}$ and $S_{\mathrm{SPM}}^{\mathrm{MM}}$ are the sinking fluxes of the SPM-bound mercury for ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$, respectively $\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}\right)$.

The photochemical rate constants ($k_{\mathrm{1}}$ and $k_{\mathrm{2}}$) are directly proportional to the shortwave radiation flux (RAD) at the water–atmosphere interface , while the biological rate constants ($k_{\mathrm{3}}$ and $k_{\mathrm{4}}$) are calculated by the organic carbon remineralization rate (OCRR) of the microbial reactions  (see Sect. S1 of the Supplement). The $k_{\text{me}}$ and $k_{\text{deme}}$ values are fixed according to , while $k_{\text{Ph-de}}$ is set according to .

The two sinking fluxes ($S_{\mathrm{SPM}}^{\mathrm{II}}$ and $S_{\mathrm{SPM}}^{\mathrm{MM}}$) are obtained according to previous works , as follows. $$\begin{array}{}\text{4}& {\displaystyle }\begin{array}{rl}& S_{\mathrm{SPM}}^{\mathrm{II}}=S_{\mathrm{POM}}^{\mathrm{II}}+S_{\mathrm{silt}}^{\mathrm{II}}=\\ & -{\displaystyle \frac{\partial }{\partial z}}\left[\mathrm{NPP}\cdot \text{(pe ratio)}\cdot {\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-\mathrm{0.9}}\cdot \left({\displaystyle \frac{k_{\mathrm{D}}}{f_{\mathrm{oc}}}}\right)\cdot {\mathrm{Hg}}^{\mathrm{II}}\left(z\right)\right]\\ & -v_{\mathrm{silt}}\cdot k_{\text{D silt}}^{\mathrm{II}}\cdot \text{SPIM}\cdot {\mathrm{Hg}}^{\mathrm{II}}\left(z\right)\end{array}\text{5}& {\displaystyle }\begin{array}{rl}& S_{\mathrm{SPM}}^{\mathrm{MM}}=S_{\mathrm{POM}}^{\mathrm{MM}}+S_{\mathrm{silt}}^{\mathrm{MM}}=\\ & -{\displaystyle \frac{\partial }{\partial z}}\left[\text{NPP}\cdot \text{(pe ratio)}\cdot {\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-\mathrm{0.9}}\cdot \left({\displaystyle \frac{k_{\mathrm{D}}}{f_{\mathrm{oc}}}}\right)\cdot \mathrm{MeHg}\left(z\right)\right]\\ & -v_{\mathrm{silt}}\cdot k_{\text{D silt}}^{\mathrm{MM}}\cdot \text{SPIM}\cdot \mathrm{MeHg}\left(z\right)\end{array}\end{array}$$ Here, $S_{\mathrm{POM}}^{\mathrm{II}}$ and $S_{\mathrm{POM}}^{\mathrm{MM}}$ are the sinking fluxes of POM-bound $\mathrm{Hg}$ for ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ $\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}\right)$, respectively; $S_{\mathrm{silt}}^{\mathrm{II}}$ and $S_{\mathrm{silt}}^{\mathrm{MM}}$ are the sinking fluxes of silt-bound $\mathrm{Hg}$ for ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ $\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}\right)$, respectively; NPP is the net primary production (g C m${}^{-\mathrm{2}}$ h${}^{-\mathrm{1}}$); the pe ratio is the ratio of particulate organic carbon (POC) export to NPP out of the euphotic zone (dimensionless); $z_{\mathrm{0}}$ is the depth of the euphotic zone (m); $k_{\mathrm{D}}$ is the seawater–SPM partition coefficient for ${\mathrm{Hg}}_{\mathrm{D}}$ (L kg${}^{-\mathrm{1}}$); $f_{\mathrm{oc}}$ is the fraction of suspended particulate matter as organic carbon (dimensionless); $v_{\mathrm{silt}}$ is the silt settling velocity (m h${}^{-\mathrm{1}}$); $k_{\text{D silt}}^{\mathrm{II}}$ is the partition coefficient of ${\mathrm{Hg}}^{\mathrm{II}}$ to silt (L kg${}^{-\mathrm{1}}$); $k_{\text{D silt}}^{\mathrm{MM}}$ is the partition coefficient of $\mathrm{MeHg}$ to silt (L kg${}^{-\mathrm{1}}$); SPIM is the suspended particulate inorganic matter (kg L${}^{-\mathrm{1}}$). The NPP is obtained by using the conversion equation for the chl $a$ concentration . This is calculated by the picoeukaryote abundance using the conversion curve of . The pe ratio is calculated by the surface atmospheric temperature measured from remote sensing and the chl $a$ concentration obtained by the NP model . The $k_{\mathrm{D}}$ value was measured within the Augusta Harbour during the last oceanographic survey, while the spatial distributions of $f_{\mathrm{oc}}$ and SPIM at the steady state have been reproduced using experimental findings for suspended particulate matter (see Sect. 2.1 of the Supplement). The $v_{\mathrm{silt}}$, $k_{\text{D silt}}^{\mathrm{II}}$, and $k_{\text{D silt}}^{\mathrm{MM}}$ values for marine environments with silty SPIM are fixed according to .

The loads of ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{II}}$ and ${\mathrm{MeHg}}_{\mathrm{D}}$ released by POM are calculated by using the following equations: $$\begin{array}{}\text{6}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{DOM}}^{\mathrm{II}}=\mathit{\lambda }\cdot m\cdot b\cdot {\mathrm{PHg}}^{\mathrm{II}},\text{7}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{DOM}}^{\mathrm{MM}}=\mathit{\lambda }\cdot m\cdot b\cdot \mathrm{PMeHg},\end{array}$$ where ${\mathrm{PHg}}^{\mathrm{II}}$ and PMeHg are the ${\mathrm{Hg}}^{\mathrm{II}}$ content and $\mathrm{MeHg}$ content, respectively, in each cell of picoeukaryotes ($\mathrm{\mu }$g per cell); $b$ is the picoeukaryote abundance (cells per cubic meter); $m$ is the mortality of picoeukaryotes (h${}^{-\mathrm{1}}$); $\mathit{\lambda }$ is the $\mathrm{Hg}$ recycling coefficient for picoeukaryotes (dimensionless). The spatiotemporal dynamics of ${\mathrm{PHg}}^{\mathrm{II}}$ and PMeHg are obtained by solving the ODEs of the Phytoplankton MERLIN-Expo model for ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$, respectively  (see Sect. S5 of the Supplement). The spatiotemporal distribution of $b$ is reproduced by using the NP model  (see Sect. S4 of the Supplement). The parameter $m$ is set according to , while $\mathit{\lambda }$ is fixed as equal to the nutrient recycling coefficient for picoeukaryotes .

The concentrations $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ are calculated as a function of position $(x,y,z)$ and time $t$, as follows: $$\begin{array}{}\text{8}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{D}}& {\displaystyle }={\mathrm{Hg}}^{\mathrm{0}}+{\mathrm{Hg}}^{\mathrm{II}}+\mathrm{MeHg},\text{9}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{T}}& {\displaystyle }={\mathrm{Hg}}_{\mathrm{D}}+k_{\mathrm{D}}\cdot \text{SPM}\cdot ({\mathrm{Hg}}^{\mathrm{II}}+\mathrm{MeHg}).\end{array}$$

Here, $k_{\mathrm{D}}$ is the seawater–SPM partition coefficient for ${\mathrm{Hg}}_{\mathrm{D}}$ (only ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$) (L kg${}^{-\mathrm{1}}$), and SPM is the suspended particulate matter concentration (kg L${}^{-\mathrm{1}}$). The partition coefficient $k_{\mathrm{D}}$ is set to the value recently experimentally observed in seawater samples collected within the Augusta Bay. The spatial distribution of SPM was set according to the experimental information collected during the oceanographic cruise of October 2017 and is assumed constant for the whole simulation time.

The advection–diffusion–reaction model is completed by a set of ordinary differential equations (ODEs), which describe the mercury fluxes at the boundaries of Augusta Harbour. Specifically, we take the following into account for the three mercury species: (i) the evasion and the deposition of ${\mathrm{Hg}}^{\mathrm{0}}$ at the water–atmosphere interface ; (ii) the lack of ${\mathrm{Hg}}^{\mathrm{0}}$ diffusion at the water–sediment interface ; (iii) the wet and dry deposition of ${\mathrm{Hg}}^{\mathrm{II}}$ at the water–atmosphere interface ; (iv) the wet and dry deposition of $\mathrm{MeHg}$ at the water–atmosphere interface ; (v) the diffusion of ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ at the water–sediment interface; (vi) the constant fixed value of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ out of Augusta Bay (Ionian Sea) ; and (vii) the exchange of elemental mercury, ${\mathrm{Hg}}^{\mathrm{II}}$, and $\mathrm{MeHg}$ between the Augusta basin and the Ionian Sea through the two inlets . Since the Augusta Bay is considered a semi-closed basin, the lateral fluxes at the boundaries of the domain are set to zero except for the two inlets . Here, the lateral fluxes depend on the direction of horizontal velocities and therefore change as a function of depth and time (see Sect. S1.1.2 of the Supplement). The boundary conditions for the three mercury species are defined by the following equations. $$\begin{array}{}\text{10}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial z}}-v_z{\mathrm{Hg}}^{\mathrm{0}}\right]\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathrm{Hg}}_{\text{gas-atm}}\cdot \mathrm{Pr}}{\mathrm{\Delta }t}}\\ & +{\mathrm{MTC}}_{\text{water-atm}}\cdot \left({\mathrm{Hg}}_{\text{gas-atm}}-H\cdot {\mathrm{Hg}}^{\mathrm{0}}{|}_{z=\mathrm{0}}\right)\end{array}\text{11}& {\displaystyle }\begin{array}{rl}& \left[D_x{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial x}}-v_x{\mathrm{Hg}}^{\mathrm{0}}\right]=\left[D_y{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial y}}-v_y{\mathrm{Hg}}^{\mathrm{0}}\right]\\ & ={\left.\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{0}}}{\partial z}}-v_z{\mathrm{Hg}}^{\mathrm{0}}\right]\right|}_{z=z_{\mathrm{b}}}=\mathrm{0}\end{array}\text{12}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial z}}-v_z{\mathrm{Hg}}^{\mathrm{II}}\right]\right|}_{z=\mathrm{0}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{II}}\cdot \mathrm{Pr}}{\mathrm{\Delta }t}}+{\text{Drydep}}_{{\mathrm{Hg}}^{\mathrm{II}}}\\ & {\left.\left[D_z{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial z}}-v_z\mathrm{MeHg}\right]\right|}_{z=\mathrm{0}}=\\ & \mathrm{0.005}\cdot {\left.\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial z}}-v_z{\mathrm{Hg}}^{\mathrm{II}}\right]\right|}_{z=\mathrm{0}}\end{array}\text{13}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_z{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial z}}-v_z{\mathrm{Hg}}^{\mathrm{II}}\right]\right|}_{z=z_{\mathrm{b}}}={\text{MTC}}_{\text{sed-water}}^{\mathrm{II}}\\ & \cdot \left({{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}-{\mathrm{Hg}}^{\mathrm{II}}{|}_{z=z_{\mathrm{b}}}\right)\end{array}\text{14}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_z{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial z}}-v_z\mathrm{MeHg}\right]\right|}_{z=z_{\mathrm{b}}}={\text{MTC}}_{\text{sed-water}}^{\mathrm{MM}}\\ & \cdot \left({\mathrm{MeHg}}_{\text{pore water}}-\mathrm{MeHg}{|}_{z=z_{\mathrm{b}}}\right)\end{array}\text{15}& {\displaystyle }\left[D_x{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial x}}-v_x{\mathrm{Hg}}^{\mathrm{II}}\right]=\left[D_y{\displaystyle \frac{\partial {\mathrm{Hg}}^{\mathrm{II}}}{\partial y}}-v_y{\mathrm{Hg}}^{\mathrm{II}}\right]=\mathrm{0}\text{16}& {\displaystyle }\begin{array}{rl}& \left[D_x{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial x}}-v_x\mathrm{MeHg}\right]\\ & =\left[D_y{\displaystyle \frac{\partial \mathrm{MeHg}}{\partial y}}-v_y\mathrm{MeHg}\right]=\mathrm{0}\end{array}\text{17}& {\displaystyle }\begin{array}{rl}& {\mathrm{Hg}}^{\mathrm{0}}(x_{\mathrm{inlet}},y_{\mathrm{inlet}},z)={\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{0}},\\ & {\mathrm{Hg}}^{\mathrm{II}}(x_{\mathrm{inlet}},y_{\mathrm{inlet}},z)={\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{II}},\\ & \mathrm{MeHg}(x_{\mathrm{inlet}},y_{\mathrm{inlet}},z)={\mathrm{MeHg}}_{\mathrm{ext}}\end{array}\end{array}$$ Here, ${\mathrm{Hg}}_{\text{gas-atm}}$ is the gaseous elemental mercury (GEM) concentration in the atmosphere (ng m${}^{-\mathrm{3}}$); Pr is the amount of precipitation (mm); $\mathrm{\Delta }t$ is the exposition time to precipitation (h); Drydep${}_{{\mathrm{Hg}}^{\mathrm{II}}}$ is the atmospheric dry deposition of ${\mathrm{Hg}}^{\mathrm{II}}$ (ng m${}^{-\mathrm{2}}$ h${}^{-\mathrm{1}}$); MTC${}_{\text{water-atm}}$ is the gas-phase overall mass transfer coefficient (m h${}^{-\mathrm{1}}$); $H$ is the Henry's law constant (dimensionless); ${\mathrm{Hg}}^{\mathrm{0}}{|}_{z=\mathrm{0}}$ is the $\left[{\mathrm{Hg}}^{\mathrm{0}}\right]$ at the sea surface ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$); ${\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{II}}$ is the $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ in the atmosphere (ng m${}^{-\mathrm{3}}$); ${\mathrm{MTC}}_{\text{sed-water}}^{\mathrm{II}}$ is the mass transfer coefficient for ${\mathrm{Hg}}^{\mathrm{II}}$ at the water–sediment interface (m h${}^{-\mathrm{1}}$); ${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$ is the $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ in the pore water of the surface layer (upper 10 cm) of the sediments ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$); ${\mathrm{Hg}}^{\mathrm{II}}{|}_{z=z_{\mathrm{b}}}$ is the dissolved $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ at the deepest layer of the water column ($z=z_{\mathrm{b}}$) ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$); MTC${}_{\text{sed-water}}^{\mathrm{MM}}$ is the mass transfer coefficient for $\mathrm{MeHg}$ at the water–sediment interface (m h${}^{-\mathrm{1}}$); MeHg${}_{\text{pore water}}$ is the [MeHg] in the pore water in the surface layer (upper 10 cm) of the sediments ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$); MeHg${|}_{z=z_{\mathrm{b}}}$ is the dissolved [MeHg] in the deepest layer of the water column ($z=z_{\mathrm{b}}$) ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$); ${\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{0}}$, ${\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{II}}$, and MeHg${}_{\text{ext}}$ are the average $\left[{\mathrm{Hg}}^{\mathrm{0}}\right]$, $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$, and [MeHg], respectively, reported from the Ionian Sea ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$). The dynamics of the GEM and ${\mathrm{Hg}}^{\mathrm{II}}$ concentrations in the atmosphere (${\mathrm{Hg}}_{\text{gas-atm}}$ and ${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{atm}}$) are reproduced using the experimental data collected in the Augusta Bay between August 2011 and June 2012 , whereas rainfall is derived from remote sensing (see http://eosweb.larc.nasa.gov, last access: March 2020). The spatiotemporal dynamics of pore water mercury concentrations (${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$ and MeHg${}_{\text{pore water}}$) at the sediment surface layer are obtained with the diffusion–reaction model for the sediment compartment, while the mass transfer coefficients (MTC${}_{\text{sed-water}}^{\mathrm{II}}$ and MTC${}_{\text{sed-water}}^{\mathrm{MM}}$) at the water–sediment interface are calculated by sediment porosity, the molecular diffusion coefficient, and boundary layer thickness above and below the sediment according to previous works  (see Sect. S1.2.2 and S1.3.2 of the Supplement).

3.1.1

Mass balance of $\mathrm{Hg}$ in Augusta Bay

The annual mass balance for the total Hg in the seawater compartment can be estimated using the boundary conditions given in Eqs. ()–(), according to the following equation : 18$$A+\mathrm{AD}+B=O+D+V,$$ where $A$ is the input of ${\mathrm{Hg}}_{\mathrm{D}}$ from anthropogenic activities; AD is the atmospheric mercury deposition; $B$ is the mercury flux from sediments to seawater due to diffusion processes; $O$ is the net mercury outflow from the Augusta Harbour to the Ionian Sea; $D$ is the amount of mercury recycled in the Augusta Bay (or the net mercury deposition for settling and burial); $V$ is the GEM evasion from the Augusta Bay to the atmosphere.

By integrating Eqs. ()–(), we obtain the terms of the annual mass balance referring to the mercury fluxes exchanged at the interfaces (AD, $B$, $V$) and the net mercury outflow from the Augusta Bay to the Ionian Sea ($O$), while the input of anthropogenic activities ($A$) is set to zero according to literature sources .

Finally, we estimate the total amount of mercury recycled ($D$) from the other terms and compare it with the amount of mercury recycled by scavenging ($S$). A simple scheme of the fluxes exchanged in the mercury biogeochemical cycle of the Augusta Bay is shown in Fig. .

3.2 The diffusion–reaction model for Hg species in pore water

The dynamics of $\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ in the Augusta sediments (average thickness of $\mathrm{1.9}$ m) have been studied using a diffusion–reaction model. In particular, we investigated the behavior of the two mercury species dissolved in pore water, i.e., ${\mathrm{Hg}}^{\mathrm{II}}$ (${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$) and $\mathrm{MeHg}$ (MeHg${}_{\text{pore water}}$), which interact with each other directly through the reaction terms of the two PDEs. Moreover, in the model we took into account the variations of mercury concentrations in pore water due to the slow desorption of the fraction bound to particulate sediments. In order to better reproduce the experimental findings, we describe mercury desorption using an exponential equation, which accounts, in the absence of external sources, for the loss of mercury through the desorption mechanism. Since mercury desorption has to depend on its instantaneous concentration, the mechanism is regulated by a first-order kinetic.

The model provides solutions for the spatiotemporal behavior of the mercury concentration for the two species dissolved in pore water, i.e., inorganic mercury (${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}(x^{\prime },y^{\prime },z^{\prime },t)$) and methyl-mercury (MeHg${}_{\text{pore water}}(x^{\prime },y^{\prime },z^{\prime },t)$), and the total mercury concentration in the sediments (${{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}(x^{\prime },y^{\prime },z^{\prime },t)$). Here, the coordinates ($x^{\prime }$, $y^{\prime }$, $z^{\prime }$) indicate the position within the irregular three-dimensional domain of the sediment compartment. Since the surface sediment slope is very low for the whole basin, the domain is approximated as the sum of several sub-domains shaped as regular parallelepipeds, which reproduce the sediment columns in each position ($x^{\prime },y^{\prime },z^{\prime }$) of the Augusta Bay. Specifically, $z^{\prime }$ represents the depth of the barycenter of each sub-domain, localized between the top ($z^{\prime }=\mathrm{0}$) and the bottom ($z^{\prime }=\mathrm{1.9}$ m) of the surface sediment layer, while the other coordinates ($x^{\prime }=x$ and $y^{\prime }=y$) indicate the distance in meters measured from the same reference point used for the seawater compartment.

In the sediment compartment, the PDEs of the model are solved by performing a centered-in-space differencing for the diffusion terms. The sea bottom is discretized in the horizontal plane using the same regular mesh adopted for simulating the dissolved mercury distribution in the seawater compartment (see Fig. 1) with $\mathrm{454.6}$ m regularly spaced elements. In this case, the vertical discretization is constituted by equally spaced layers of 0.2 m depth, with the exception of the interface layer between the water and sediment, whose depth is set at 0.1 m. This choice has been made in order to best adapt the 3D grid of the model to the scheme used to interpolate available experimental data.

The initial conditions for $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ are fixed on the basis of experimental findings. As a first step we reproduced the spatial distribution of ${{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}$ at time $t=\mathrm{0}$ by interpolating the experimental data collected by ICRAM in $\mathrm{2005}$  (see Sect. S1.2.4 of the Supplement). We then calculated both $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ and [MeHg] in pore water using the following equations:

19$$\begin{array}{rl}& {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\left(\mathrm{0}\right)=(\mathrm{1}-f_{\mathrm{MeHg}})\cdot {\displaystyle \frac{{{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(\mathrm{0}\right)}{K_{\mathrm{d}}^{\mathrm{II}}}},\\ & {\mathrm{MeHg}}_{\text{pore water}}\left(\mathrm{0}\right)=f_{\mathrm{MeHg}}\cdot {\displaystyle \frac{{{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(\mathrm{0}\right)}{K_{\mathrm{d}}^{\mathrm{MM}}}},\end{array}$$ where ${{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(\mathrm{0}\right)$ represents the spatial distribution of $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ in the sediments at an initial time (mg kg${}^{-\mathrm{1}}$), $f_{\mathrm{MeHg}}$ is the fraction of $\mathrm{MeHg}$ in the sediments (dimensionless), $K_{\mathrm{d}}^{\mathrm{II}}$ is the sediment–pore water distribution coefficient for ${\mathrm{Hg}}^{\mathrm{II}}$ (L kg${}^{-\mathrm{1}}$), and $K_{\mathrm{d}}^{\mathrm{MM}}$ is the sediment–pore water distribution coefficient for $\mathrm{MeHg}$ (L kg${}^{-\mathrm{1}}$).

In pore water, the dynamics of $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ and [MeHg] are modeled by considering three chemical–physical processes : (i) methylation and demethylation (reaction terms); (ii) passive movement due to the Brownian motion of each chemical species (diffusion terms); and (iii) the desorption of mercury bound to sediment particles (desorption term).

The methylation and demethylation processes involved in the dynamics of ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ are considered in the model through reaction terms describing first-order kinetics. The rate constants of these reactions are fixed according to previous works .

The diffusion terms reproduce the effects of Brownian motion on the spatial distribution of $\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ in pore water. In particular, the magnitude of the Brownian motion is described by the molecular diffusion coefficients for ${\mathrm{Hg}}^{\mathrm{II}}$ ($D_{\mathrm{sed}}^{\mathrm{in}}(x^{\prime },y^{\prime },z^{\prime })$) and $\mathrm{MeHg}$ ($D_{\mathrm{sed}}^{\mathrm{or}}(x^{\prime },y^{\prime },z^{\prime })$), which change in each position of the domain as a function of porosity and tortuosity (see Sect. S1.2.2 and S1.3.2 of the Supplement). The molecular diffusion coefficients are assumed isotropic in all directions and are set as constant functions of time according to previous works .

The desorption term estimates the increase in ${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$ and MeHg${}_{\text{pore water}}$ due to the mercury release from the sediment particles to pore water. The desorption process is regulated by the temporal gradient of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ ($\partial {\mathrm{dHg}}_{\mathrm{T}}^{\mathrm{sed}}/\partial \mathrm{d}t$), which changes as a function of position and time (see Sect. S1.2.2 and S1.3.2 of the Supplement).

Thus, the module for the sediment compartment is defined by the following coupled partial differential equations. $$\begin{array}{}\text{20}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathrm{dHg}}_{\mathrm{pore}\mathrm{water}}^{\mathrm{II}}}{\mathrm{d}t}}=+K_{\text{demeth}}\cdot {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}-K_{\text{meth}}\\ & \cdot {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}+{\displaystyle \frac{\partial }{\partial x}}\left[D_{\mathrm{sed}}^{\mathrm{in}}\cdot {\displaystyle \frac{\partial {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}}{\partial x}}\right]\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[D_{\mathrm{sed}}^{\mathrm{in}}\cdot {\displaystyle \frac{\partial {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}}{\partial y}}\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{sed}}^{\mathrm{in}}\cdot {\displaystyle \frac{\partial {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}}{\partial z}}\right]-{\displaystyle \frac{(\mathrm{1}-f_{\mathrm{MeHg}})}{K_{\mathrm{d}}^{\mathrm{II}}}}\\ & \cdot {\displaystyle \frac{{\mathrm{dHg}}_{\mathrm{T}}^{\mathrm{sed}}}{\mathrm{d}t}}\end{array}\text{21}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathrm{dMeHg}}_{\mathrm{pore}\mathrm{water}}}{\mathrm{d}t}}=-K_{\text{demeth}}\cdot {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}\\ & +K_{\text{meth}}\cdot {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}+{\displaystyle \frac{\partial }{\partial x}}\left[D_{\mathrm{sed}}^{\mathrm{or}}\cdot {\displaystyle \frac{\partial {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}}{\partial x}}\right]\\ & +{\displaystyle \frac{\partial }{\partial y}}\left[D_{\mathrm{sed}}^{\mathrm{or}}\cdot {\displaystyle \frac{\partial {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}}{\partial y}}\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{sed}}^{\mathrm{or}}\cdot {\displaystyle \frac{\partial {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}}{\partial z}}\right]-{\displaystyle \frac{f_{\mathrm{MeHg}}}{K_{\mathrm{d}}^{\mathrm{MM}}}}\cdot {\displaystyle \frac{{\mathrm{dHg}}_{\mathrm{T}}^{\mathrm{sed}}}{\mathrm{d}t}}\end{array}\text{22}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathrm{dHg}}_{\mathrm{T}}^{\mathrm{sed}}}{\mathrm{d}t}}=-\mathit{\alpha }\cdot {{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\Rightarrow {{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(t\right)={{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(\mathrm{0}\right)\\ & \cdot \exp (-\mathit{\alpha }\cdot t)\end{array}\end{array}$$ Here, $K_{\text{demeth}}$ is the rate constant for the demethylation of $\mathrm{MeHg}$ (h${}^{-\mathrm{1}}$); $K_{\text{meth}}$ is the rate constant for the methylation of ${\mathrm{Hg}}^{\mathrm{II}}$ (h${}^{-\mathrm{1}}$); $\mathit{\alpha }$ is the desorption rate for the $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ bound to the sediment particles (h${}^{-\mathrm{1}}$).

As boundary conditions, we assume a null value of mercury flux at the bottom of the sediment column ($\mathrm{1.9}$ m of depth), mainly due to the measured very low porosity, while the vertical gradients of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ are set to zero at the water–sediment interface according to field observations. The mercury concentration in sediments is fixed to zero at the lateral boundaries $(x_{\mathrm{b}}^{\prime },y_{\mathrm{b}}^{\prime })$ of the 3D domain. The boundary conditions for dissolved and total mercury concentrations in sediments are described by the following equations. $$\begin{array}{}\text{23}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_{\mathrm{sed}}^{\mathrm{in}}{\displaystyle \frac{\partial {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{0}}\\ & ={\left.\left[D_{\mathrm{sed}}^{\mathrm{in}}{\displaystyle \frac{\partial {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{1.9}\mathrm{m}}=\mathrm{0}\end{array}\text{24}& {\displaystyle }\begin{array}{rl}& {\left.\left[D_{\mathrm{sed}}^{\mathrm{or}}{\displaystyle \frac{\partial {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{0}}=\\ & {\left.\left[D_{\mathrm{sed}}^{\mathrm{or}}{\displaystyle \frac{\partial {\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{1.9}\mathrm{m}}=\mathrm{0}\end{array}\text{25}& {\displaystyle }\begin{array}{rl}& {\left.\left[{\displaystyle \frac{\partial {{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{0}}=\\ & {\left.\left[{\displaystyle \frac{\partial {{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}}{\partial z}}\right]\right|}_{z^{\prime }=\mathrm{1.9}\mathrm{m}}=\mathrm{0}\end{array}\text{26}& {\displaystyle }\begin{array}{rl}& {{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}{|}_{(x_{\mathrm{b}}^{\prime },y_{\mathrm{b}}^{\prime })}=\mathrm{0},{\mathrm{MeHg}}_{\mathrm{pore}\mathrm{water}}{|}_{(x_{\mathrm{b}}^{\prime },y_{\mathrm{b}}^{\prime })}=\mathrm{0},\\ & {{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}{|}_{(x_{\mathrm{b}}^{\prime },y_{\mathrm{b}}^{\prime })}=\mathrm{0}\end{array}\end{array}$$ Equations ()–() represent the three-dimensional diffusion–reaction model used to describe and reproduce the spatiotemporal dynamics of $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ and [MeHg] in pore water, as well as of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ in sediments. It should be noted that Eqs. ()–(), ()–(), and (), which reproduce the spatiotemporal distributions of the mercury concentrations in both compartments (seawater and sediment), strongly depend on the initial condition for the total mercury concentration observed in the sediments.

3.3 Model and simulation setup

In our model, as initial conditions we assumed a uniformly distributed concentration of ${\mathrm{Hg}}_{\mathrm{D}}$ and ${\mathrm{Hg}}_{\mathrm{T}}$, set to $\mathrm{1.9}$ ng L${}^{-\mathrm{1}}$, corresponding to the experimental detection limit. Specifically, the initial concentration of each $\mathrm{Hg}$ species was fixed according to the percentage observed in seawater samples from the Ionian Sea  in such a way as to respect the detection limit for total $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$. The numerical results were not affected by the chosen initial conditions; indeed, the same spatial distribution of $\left[\mathrm{Hg}\right]$ at nearly steady state was obtained when initial $\mathrm{Hg}$ concentrations higher than the detection limit were fixed.

The model results were obtained by running a single long simulation. To reproduce the spatial mercury distributions at nearly steady state, we integrated the model equations over a time interval ($t_{max}>\mathrm{7}$ years) long enough to reach an annual decrease in the mercury concentration of less than 2 %. This percentage value progressively declines for longer time intervals down to an annual decrease of 0.12 % for $t_{max}=\mathrm{250}$ years.

All environmental parameters and variables used in the model are reported in Tables S1–S3 of the Supplement. Most of the environmental parameters were set to values experimentally observed at sites contaminated by mercury , while other parameters, including those most sensitive for the model, were calibrated to correctly reproduce the experimental data collected during the six oceanographic surveys . Furthermore, the photochemical and biological rate constants of the redox reactions were calculated using both the outputs of the NP model and the data from remote sensing (see Sect. S1 of the Supplement).

Experimental measurements were carried out during the period between $\mathrm{2005}$ and $\mathrm{2017}$ at several stations inside and outside Augusta Harbour (see Fig. 1). The mercury concentration and mercury fluxes were measured both in sediments and seawater (see Tables S6–S10 and S12 of the Supplement). We refer to , , and for a detailed description of the measured parameters, the related dynamics, and the analytical methods used . These experimental data were used to identify the most sensitive parameters for the model and to compare them with the theoretical results in order to estimate the model accuracy in reproducing Hg dynamics.

Concerning the calibration procedure, we first focused on the best values of the parameters for the sediment compartment (i.e., sediment–pore water distribution coefficients, desorption rate, and boundary layer thickness above the sediment) in such a way as to optimize the match between theoretical results and experimental observations. Specifically, in Eqs. ()–() the sediment–pore water distribution coefficients were calibrated to guarantee the best theoretical $\left[\mathrm{Hg}\right]$ in pore water in agreement with the value ranges experimentally observed in a previous work , whereas the fraction of methyl-mercury in sediments for the whole spatial domain was set to that obtained by field observations during the oceanographic survey of October $\mathrm{2017}$ (see Table S1). In Eq. (), the desorption rate $\mathit{\alpha }$ was calibrated to obtain the best fit between the theoretical results and experimental observations for $\left[\mathrm{Hg}\right]$ in pore water. Before calculating the mass transfer coefficients at the water–sediment interface, the boundary layer thickness above the sediment was optimized to better reproduce the spatial distribution of the mercury benthic flux observed experimentally. Unlike the boundary layer thickness above the sediment, the other parameters used to obtain ${\mathrm{MTC}}_{\text{sed-water}}^{\mathrm{II}}$ and ${\mathrm{MTC}}_{\text{sed-water}}^{\mathrm{MM}}$ were not calibrated. In fact, the boundary layer thickness below the sediment was estimated by using the relationship between this parameter and the average velocity of marine currents defined by , while the spatial distribution of the sediment porosity within Augusta Harbour was reproduced according to previous works  by exploiting the measurements on the sediment samples performed by ICRAM in $\mathrm{2005}$. Also, the molecular diffusion coefficient was that reported by .

As a second step, we calibrated model parameters for the seawater compartment (i.e., vertical and horizontal diffusivities) in order to better reproduce the spatiotemporal dynamics of the dissolved mercury concentration. The vertical turbulent diffusivity was calibrated according to experimental data, which indicate weakly mixed water column conditions within the Augusta Bay during the whole year. Specifically, the vertical turbulent diffusivity was set in such a way as to obtain the best match with the experimentally observed dissolved mercury concentration at the surface layer of the water column. The calibrated vertical diffusivity was in good agreement with previously reported values  under the condition of weakly mixed waters. The horizontal turbulent diffusivity was assumed isotropic in the horizontal water plane ($D_x=D_y$) and calibrated by considering the values obtained in . In particular, the horizontal turbulent diffusivities were optimized to obtain the best possible match with the observed mercury evasion flux. The calibrated horizontal diffusivities were in accordance with the values estimated by other authors  for basins similar in size to those of the Augusta Bay.

In our analysis, no comparison between the calibrated desorption rate and experimental data was possible. However, the other calibrated environmental parameters were in good agreement with those obtained experimentally in the Augusta Bay and at other sites contaminated by mercury .

Finally, the calibrated model was run by considering the seasonal oscillations of the environmental data (water currents, wind, etc.) provided by hydrodynamic modeling (see Sect. S3 of the Supplement). The model results were validated using the other experimental findings acquired in the Augusta Bay: $\left[{\mathrm{MeHg}}_{\mathrm{D}}\right]$, $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$, and $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ measured in seawater and all annual $\mathrm{Hg}$ fluxes estimated for the mass balance.

4 Results

In the following the simulation results obtained for the seawater and sediment compartments are described and compared with experimental data.

4.1 Mercury in seawater

The spatial distribution of the three mercury species dissolved in seawater is obtained by solving Eqs. ()–(), together with the equation system ()–(), for the sediment compartment.

Figure 4

Comparison between the experimental data (red points) and the theoretical results (black lines) for the dissolved mercury concentration at six stations in the Augusta Bay. The vertical profiles of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ are obtained by the model for the sites closest to stations $\mathrm{1}$, $\mathrm{10}$, $\mathrm{12}$, and $\mathrm{20}$ (sampling May 2011), station $\mathrm{15}$ (sampling June 2012), and station $\mathrm{26}$ (sampling February 2012).

[Figure omitted. See PDF]

The theoretical concentrations of the three mercury species dissolved in seawater are reported in Table S5 of the Supplement. Here, we observe that the average concentration ratios among the three mercury species dissolved in seawater ($\left[{\mathrm{Hg}}^{\mathrm{II}}\right]/\left[\right[{\mathrm{Hg}}_{\mathrm{D}}]=\mathrm{0.790}$, $\left[\mathrm{MMHg}\right]/\left[\right[{\mathrm{Hg}}_{\mathrm{D}}]=\mathrm{0.022}$, and $\left[{\mathrm{Hg}}^{\mathrm{0}}\right]/\left[\right[{\mathrm{Hg}}_{\mathrm{D}}]=\mathrm{0.188}$) are in good agreement with both the experimental and theoretical values reported in recent publications . Moreover, the theoretical results for the vertical profiles of the mercury concentration show a similar shape for the whole simulated period (2005–2254), while the magnitude of the concentrations in the whole water column decreases slowly as a function of time (see Fig. S1 of the Supplement).

The model results indicate that the dissolved mercury concentration is usually maximal at the seawater–sediment interface (see Fig. 4) where the main sources of ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ are localized. These numerical results are in reasonable agreement with the field observations (see Tables S6–S7 of the Supplement). Moreover, taking into account the redox conditions of sediments in the area, we speculate that maxima in $\mathrm{MeHg}$ production are confined to the seawater–sediment interface.

Conversely, at some sites of the calculation grid (see Fig. 1) we observe that the peaks of the mercury concentration occur at mid-depth in the water column, possibly due to the distribution of the marine current velocity field within the Augusta basin, which sometimes determines the presence of an $\left[\mathrm{Hg}\right]$ maximum in the intermediate layers of seawater. In general, in our model the dynamics of the mercury concentration in seawater are strictly connected to the behavior of the benthic mercury fluxes, which decrease slowly as a function of time due to the slow molecular diffusion process of mercury within the pore waters of the sediments.

A quantitative analysis, based on the reduced ${\mathit{\chi }}^{\mathrm{2}}$ test, indicates good agreement between the model results and experimental findings for [MeHg] at stations $A\mathrm{3}$ (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0005}$) and $A\mathrm{7}$ (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0005}$), while differences can be observed at stations $A\mathrm{9}$ (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0955}$) and $A\mathrm{11}$ (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.1065}$), where the theoretical concentrations appear overestimated at the bottom layer (see Table S6 of the Supplement). This result is probably due to the overestimation of the $\mathrm{MeHg}$ benthic fluxes at these two stations.

Figure 5

Distribution of $\mathrm{MeHg}$ and ${\mathrm{Hg}}_{\mathrm{D}}$ fluxes calculated at the seawater–sediment interface. The maps reproduce the spatial distribution of the benthic flux in the Augusta Bay during the two sampling periods, i.e., 19–21 September 2011 (a, b) and 23–26 June 2012 (c, d).

[Figure omitted. See PDF]

In our analysis, the spatiotemporal behavior of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ is obtained as the sum of the three dissolved mercury species. On the other hand, the dynamics of the spatial distribution of $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ are estimated according to Eq. (), assuming a linear correlation between the modeled $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ and the experimental SPM concentrations. The spatial distributions of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ and $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ are reported for May 2011 in Figs. S2–S5 of the Supplement.

The numerical results for $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ are in good agreement with the experimental findings for the four investigated periods (see Table S7 of the Supplement). Specifically, the difference between the model result and field observation for $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ is less than the experimental error ($\mathit{\sigma }=\mathrm{3.2}$ ng L${}^{-\mathrm{1}}$) in 59 % of sampling points, while it exceeds $\mathrm{2}\mathit{\sigma }$ in only 17 % of sampling sites. As a conclusion, the comparison between experimental data and theoretical results for $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ shows mostly small discrepancies except in some of the most contaminated areas, where concentration hot spots are hard to capture due to the resolution grid used in the present work.

The model results for $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ show some discrepancies with experimental data at most of the sites investigated during the first sampling period (May 2011), while in general they evidence acceptable agreement for the other sampling periods (see Table S8 of the Supplement). As a whole, the discrepancy for $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ is less than $\mathit{\sigma }$ in 44 % of cases, while it exceeds $\mathrm{2}\mathit{\sigma }$ at 32 % of sampling sites. The differences (larger than $\mathit{\sigma }=\mathrm{3.2}$ ng L${}^{-\mathrm{1}}$) can mainly be explained by the significant distance between the sampling sites and the model calculation grid nodes (see Fig. 1). Additionally, we cannot neglect the role played by the theoretical spatial distribution of the SPM concentration (see Eq. ), which could significantly affect the spatiotemporal dynamics of the total mercury concentration in seawater. In particular, the spatial distribution of SPM concentrations used in the model is probably not appropriate for the first sampling period investigated (May 2011), while it produces good agreement for the other three sampling periods.

Figure 6

Distribution of the ${\mathrm{Hg}}^{\mathrm{0}}$ flux calculated at the seawater–atmosphere interface. The maps reproduce the spatial distribution of the evasion flux in the Augusta Bay during the two sampling periods, i.e., 29–30 November 2011 (a) and 23–25 June 2012 (b).

[Figure omitted. See PDF]

The theoretical distributions of the benthic mercury fluxes simulated by the model for the two investigated periods (September 2011 and June 2012) are shown in Fig. . Here, very high benthic ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$ fluxes are documented in the southwest sector of Augusta Harbour, where the chlor-alkali plant discharged high amounts of contaminants until the late 1970s. The model also reliably reproduces the high benthic mercury fluxes in the part of the southeast sector close to the inlets of the Augusta Bay. The benthic mercury fluxes are very low in the coastal zones at the north of the basin, while intermediate values have been calculated in the central part of the bay. As a whole, the estimated benthic mercury fluxes are in good agreement with the experimental data collected during the two sampling periods (see Table S9 of the Supplement). It should be noted that the model results suggest that the benthic ${\mathrm{Hg}}_{\mathrm{D}}$ fluxes are mainly generated by the diffusion process at the seawater–sediment interface and that the amount of ${\mathrm{Hg}}_{\mathrm{D}}$ release from the resuspended particulate matter is negligible. Moreover, the model results confirm that the spatial heterogeneity of the benthic fluxes observed experimentally is strictly connected to that of the ${\mathrm{Hg}}_{\mathrm{T}}$ concentration in sediments.

In general, the theoretical distribution of the mercury evasion fluxes is in acceptable agreement with the experimental results for the investigated periods (see Table S10 of the Supplement). Specifically, small discrepancies are observed at most stations (four out of six), while larger difference emerge at stations $\mathrm{3}$ (November $\mathrm{2011}$) and $\mathrm{5}$ (June $\mathrm{2012}$). From a qualitative point of view, the model results for elemental mercury evasion confirm that a high flux is present in the coastal zones at the southwest of the Augusta Bay , while a reduced evasion flux is observed at the northern sector of the basin (see Fig. ).

Figure 7

Mercury benthic fluxes (a), evasion flux to the atmosphere (b), net outflows at the inlets (c), and recycling fluxes (d).

[Figure omitted. See PDF]

In this work, we calculate the annual mass balance of the Augusta Bay to study the fate of $\mathrm{Hg}$ coming from sediments and to estimate the $\mathrm{Hg}$ outflows at the inlets of basin. In Fig. , we show the temporal behavior of the annual mercury fluxes used for the mass balance calculation (see also Table S11 of the Supplement). The results of the annual benthic mercury fluxes ($B$) show that most of the mercury coming up from sediments is in inorganic form (see Fig. 7a), while the benthic $\mathrm{MeHg}$ flux appears to be 1 to 2 orders of magnitude lower. The model results are compared with experimental information reported by for three different sampling sites and in two different periods (September 2011 and June 2012). The modeled ${\mathrm{Hg}}_{\mathrm{D}}$ benthic fluxes (2.65 kmol yr${}^{-\mathrm{1}}$ for the year 2011 and 2.61 kmol yr${}^{-\mathrm{1}}$ for the year 2012) are significantly larger than those estimated for both sampling periods on the basis of the field observations (1.1 kmol yr${}^{-\mathrm{1}}$ in September 2011 and 1.4 kmol yr${}^{-\mathrm{1}}$ in June 2012) . This probably depends on the limited number of sampling sites available in the experimental work with a consequent limited capacity to capture reliable estimates of benthic fluxes within a basin, such as Augusta Bay, where the spatial distribution of sediment mercury is highly heterogeneous. Moreover, the higher resolution of the grid used in our model guarantees a better estimation of the annual benthic mercury fluxes once the spatiotemporal integration is performed.

The model results for the dynamics of the annual mercury evasion fluxes are shown in Fig. 7b. The comparison with experimental findings indicates that the mercury evasion fluxes ($V$) obtained from the model ($\mathrm{1.93}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr${}^{-\mathrm{1}}$ for the year 2011 and $\mathrm{1.90}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr${}^{-\mathrm{1}}$ for the year 2012) are in good agreement with those estimated by for each year ($\mathrm{1.70}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr${}^{-\mathrm{1}}$). Conversely, a significant discrepancy is observed between the annual atmospheric mercury deposition (AD) obtained by our model ($\mathrm{0.22}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr${}^{-\mathrm{1}}$) and that estimated in the experimental work ($\mathrm{0.42}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr${}^{-\mathrm{1}}$) . This discrepancy is due to the different calculation methods used in the two works. Specifically, in our model the AD is calculated by using both the atmospheric mercury concentrations and the average precipitation measured for all months of the year. In contrast, in the AD is calculated by averaging the experimental data acquired during a time-limited sampling period (from $\mathrm{29}$ August $\mathrm{2011}$ to $\mathrm{23}$ April $\mathrm{2012}$), namely without considering the period in which the amount of precipitation is very low. In this way, the AD obtained by is much higher than that of our model, even if it is probably overestimated due to the calculation method used. In general, the contribution of AD is negligible in the mercury mass balance of the Augusta Bay. Indeed, the simulations indicate that a strong increase in atmospheric mercury deposition caused by environmental changes (dust fall increase and/or rainfall increase) would not affect the numerical results of our model significantly.

The annual net mercury inflows ($A$) from rivers and sewerage to the basin are assumed to be negligible, in agreement with field observations. Specifically, the flow rate of the Marcellino River is equal to zero for most of the year, while the inflow from sewerage is low. Moreover, it is fair to speculate that the $\mathrm{Hg}$ concentration in fresh waters discharged into the Augusta Bay decreased significantly after the chlor-alkali plant closure.

The dynamics of the annual net mercury outflow ($O$) at the Levante and Scirocco inlets are described in Fig. 7c. The results encompass both the inflow and outflow of the water mass in each inlet for the whole year and the associated ${\mathrm{Hg}}_{\mathrm{T}}$ contribution. In Fig. 7c, we show the annual ${\mathrm{Hg}}_{\mathrm{T}}$ outflow from the Augusta Bay towards the open sea. This has been estimated to be 0.13 kmol yr${}^{-\mathrm{1}}$ for the year 2012 and appears significantly lower than the 0.51 kmol yr${}^{-\mathrm{1}}$ calculated by for the same year. Our hypothesis to explain this discrepancy is that the previous study does not consider the dynamics of $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ at the inlets (the ${\mathrm{Hg}}_{\mathrm{T}}$ outflow is calculated only on the basis of the mercury concentration measured in February 2012) and that the approach used in the previous paper does not take into account the dynamics of the inflow and outflow of the water mass at the two inlets.

In this work, the annual recycled mercury flux ($D$) is calculated by subtraction using the mass balance Eq. (), as was done in previous works on the Augusta Bay . The model results for the recycled mercury flux are shown in Fig. 7d. Here, the values calculated by our model (2.50 kmol yr${}^{-\mathrm{1}}$ for the year 2011 and 2.46 kmol yr${}^{-\mathrm{1}}$ for the year 2012) are larger and probably more realistic than those estimated in (0.84 kmol yr${}^{-\mathrm{1}}$).   Indeed, the former are obtained by considering the seasonal oscillations of all other mercury fluxes during the year, while the latter are calculated without considering the seasonal changes in mercury fluxes .

In order to reproduce the effects induced by the scavenging process on mercury dynamics, our model calculates the annual sinking mercury flux, whose results are shown in Fig. 6d. Here, a significant gap between the recycled flux (2.50 kmol yr${}^{-\mathrm{1}}$ for the year 2011) and the sinking flux (0.07 kmol yr${}^{-\mathrm{1}}$ for the year 2011) is observed, probably due to the underestimation of the amount of mercury captured by POM (see Eqs. –). More specifically, this behavior could be caused by the underestimation of NPP, which is calculated by using a conversion equation calibrated for oceans  rather than for coastal zones.

In contrast, very high values of the annual ${\mathrm{Hg}}_{\mathrm{T}}$ accumulation rate in the surface sediment layer (12.07 kmol yr${}^{-\mathrm{1}}$ for the year 2011), with respect to those of the annual recycled flux (2.50 kmol yr${}^{-\mathrm{1}}$ for the year 2011), are obtained by our model. This result is caused by the high sedimentation rate (11.7 mm yr${}^{-\mathrm{1}}$) estimated experimentally  and used in our calculations for the annual ${\mathrm{Hg}}_{\mathrm{T}}$ accumulation rate . However, the sedimentation rate could be overestimated due to the sampling methods used. In fact, the results obtained by the sediment transport model indicate a low average sedimentation rate for the Augusta Bay.

4.2 Mercury in sediments

The spatiotemporal dynamics of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ in the sediments of Augusta Bay and the mercury concentration of the two species (${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$) dissolved in pore water have been obtained by solving Eqs. ()–(). All environmental parameters and variables used for the sediment compartment are reported in Tables S1–S2 of the Supplement.

Figure 8

Dynamics of the vertical profiles of $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ in sediments (a, d) as well as $\left[{\mathrm{Hg}}^{\mathrm{II}}\right]$ and [MeHg] in pore waters (b, c, e, f) at stations 8 and 16 (sampling May 2011).

[Figure omitted. See PDF]

In Fig. , the vertical profiles of the mercury concentration in the sediments indicate that $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$, $\left[{{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\right]$, and $\left[{\mathrm{MeHg}}_{\text{pore water}}\right]$ always reach their maximum value within the surface layer of the sediments ($<\mathrm{0.5}$ m of depth). However, the shape of the vertical profiles for $\left[{{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\right]$ and $\left[{\mathrm{MeHg}}_{\text{pore water}}\right]$ in pore water changes as a function of time. Also, the magnitude of the concentration peaks decreases over the whole 3D domain during the period studied. In particular, the pore water mercury concentration assumes a nearly uniform distribution along the whole sediment column after several years of model simulation, even if the highest mercury concentrations are always observed in the shallowest layer of the sediments.

The highest $\left[{{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\right]$ and $\left[{\mathrm{MeHg}}_{\text{pore water}}\right]$ in the sediment surface layer support the high benthic mercury fluxes measured even several years after the chlor-alkali plant closure. Moreover, the results of $\left[{{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\right]$ and $\left[{\mathrm{MeHg}}_{\text{pore water}}\right]$ also indicate that the benthic mercury fluxes will remain elevated until the beginning of the 23rd century.

Finally, the comparison performed for $\left[{\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{sed}}\right]$ in pore water indicates good agreement between the theoretical results and the experimental data (see Table S12 of the Supplement).

5 Discussion

In this work we introduced the innovative HR3DHG biogeochemical model, which has been verified and validated for all its modules with the rich database acquired for the Augusta Bay. The model is an advection–diffusion–reaction model  that reproduces the spatiotemporal dynamics of the mercury concentration in seawater. The advection–diffusion–reaction model was coupled with (i) a diffusion–reaction model, which estimates the mercury concentration in the pore waters of the sediment compartment, and (ii) the equation that reproduces the mechanism responsible for the desorption of the two mercury species from the solid to the liquid phase of the sediments. This “integrated” model, which allows us to give a description of the mercury dynamics at highly polluted marine sites, introduces some novelties in the landscape of the mathematical modeling of spatiotemporal dynamics in a biogeochemical context.

This integrated model also estimates the total amount of mercury present in biological species that occupy the lowest trophic level of the food chain, i.e., phytoplankton populations. For this purpose, we incorporated the Phytoplankton MERLIN-Expo model  to describe the mechanism of mercury uptake in phytoplankton cells. Moreover, we reproduced the spatiotemporal dynamics of phytoplankton communities in seawater using the Nutrient–Phytoplankton model . This integrated model, together with the Nutrient–Phytoplankton model and the Phytoplankton MERLIN-Expo model, constitutes a new global biogeochemical (HR3DHG) model describing mercury dynamics and their effects on the lowest level of the trophic chain.

The HR3DHG model simultaneously provides the high-resolution spatiotemporal dynamics of $\left[\mathrm{Hg}\right]$ in seawater and sediment, as well as $\mathrm{Hg}$ fluxes at the boundaries of the 3D domain. The former are useful to locate the most polluted areas within the investigated basin. The latter are necessary to obtain the annual mercury mass balance of the basin in the quasi-stationary condition and to predict the mercury outflow towards the open sea, even after a very long time.

For comparison, the different approach used in the WASP models did not allow us to reproduce the dynamics of the mercury concentration distribution at 3D high resolution at polluted sites characterized by elevated spatial heterogeneity. Similar criticalities came out from the study of HR 1D models , in which the effects of the horizontal velocity field on the mercury dynamics could not be taken into account. Moreover, both the mechanism of the desorption of the total mercury in sediments and the processes involved in dissolved mercury dynamics in pore water were not considered in most advection–diffusion–reaction models, such as the BROM. In general, only a few models  were able to make forecasts about the mercury depletion time in the sediment compartment of highly polluted sites, such as Augusta Bay.

Finally, the biogeochemical models introduced in previous works  provided neither the NPP from the Nutrient–Phytoplankton model  nor the load of POM-released ${\mathrm{Hg}}_{\mathrm{D}}$ obtained using the Phytoplankton MERLIN-Expo model  (see Sect. 3.1).

All the aforementioned aspects are therefore elements of novelty in the context of 3D biogeochemical modeling. The HR3DHG model considers the effects of seasonal changes in the environmental variables on the mercury outflows towards the atmosphere and the open sea.

Application of the HR3DHG model to the case study of Augusta Bay provides crucial information for that environment, helping us to revise our view of the mercury dynamics at the highly contaminated coastal marine sites of the Mediterranean Sea. Firstly, the mass transfer coefficients at the water–sediment interface are highly sensitive to the layer thickness above the sediment. Specifically, for each mercury species in sediments, a small decrease in this parameter causes a great increase in benthic fluxes, with a consequent strong enhancement of the dissolved mercury concentration in seawater.

The model framework for the sediment compartment causes the spatiotemporal dynamics of the benthic mercury flux to strongly depend on the spatial distribution of the sediment porosity and the initial total mercury concentration in the top sediments, which were fixed using the experimental data.

A sensitivity analysis performed on the environmental parameters and variables used in the seawater compartment indicates that the spatiotemporal dynamics of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ primarily depend on the velocity field of the marine currents obtained from the hydrodynamic model , even if the role played by the vertical and horizontal diffusivities  cannot be neglected. Specifically, the spatiotemporal behavior of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ changed significantly when alternative velocity fields for the Augusta Bay were used in the biogeochemical module, confirming a feature already observed in previous models . Conversely, limited changes in the spatial distribution of $\left[{\mathrm{Hg}}_{\mathrm{D}}\right]$ were observed when different values of the vertical and horizontal diffusivities were set in our model.

The magnitude of the elemental mercury concentration is tightly connected to the values assigned to the rate constants of the photochemical redox reactions, while the role played by the other reaction rates appears negligible for this mercury species.

According to the available experimental data, the theoretical results obtained with the HR3DHG model suggest that the amount of mercury bound to particulate matter is quite high in the seawater compartment (about $\mathrm{47}$ % of the ${\mathrm{Hg}}_{\mathrm{T}}$ on average). Because of the exponential decay of $\left[{\mathrm{Hg}}_{\mathrm{T}}\right]$ in sediments, the concentration of the three mercury species dissolved in seawater decreases slowly as a function of time, whereas their concentration ratios remain approximately constant. Specifically, the mean concentrations of mercury are partitioned as $\mathrm{79.0}$ % of ${\mathrm{Hg}}^{\mathrm{II}}$, $\mathrm{18.8}$ % of elemental mercury, and $\mathrm{2.2}$ % of $\mathrm{MeHg}$, namely values very similar to those observed experimentally in other contaminated sites . The same ratio is observed for mercury that outflows from the inlets of Augusta Bay to the open sea. Here, the theoretical results of the HR3DHG model show a progressive decrease in annual mercury outflow from the bay.

On the whole, the mercury dissolved in seawater is derived from sediments through the benthic flux of ${\mathrm{Hg}}^{\mathrm{II}}$ and $\mathrm{MeHg}$. In particular, these two mercury species are released directly by the sediments, while elemental mercury is generated by redox reactions that involve the other two species. The elemental mercury concentration at the water surface contributes to the mercury evasion flux, even if only a small part of the elemental mercury in the seawater is released in the atmosphere.

Notably, the theoretical results of the HR3DHG model demonstrate the pivotal role played by the recycling process in the mercury mass balance of the Augusta Bay. Estimates for the annual recycled mercury flux indicate that most ($\mathrm{94}$ %) of the mercury released by sediments remains within the Augusta basin, while the mercury outflows at the boundaries of the basin are negligible with respect to the annual benthic mercury fluxes. More specifically, in the quasi-stationary condition, the model results indicate that most of the recycled mercury returns to the sediments where it is reburied, and the amount of mercury absorbed by POM (0.008 kmol yr${}^{-\mathrm{1}}$ for the year 2011) and recycled in seawater is negligible. In this last respect, however, it is important to underscore that even a reduced amount of $\mathrm{MeHg}$ entering living phytoplankton cells can be very dangerous for the health of human beings due to the bioaccumulation processes that occur throughout the food chain .

The theoretical results show that the recycled mercury flux in the Augusta Bay is only partially described by the scavenging process. In particular, an underestimation of the sinking flux for POM-bound mercury is observed when the NPP from the NP model is used in Eqs. ()–(). This behavior is probably due to the chl $a$ concentration conversion equation of , which has been calibrated for oceans instead of coastal zones. For this reason, the NPP estimation would need further experimental and theoretical investigations. Moreover, deeper knowledge of the scavenging process, which determines the particulate $\mathrm{Hg}$ dynamics, would be necessary from a theoretical point of view to obtain a better estimation of the ${\mathrm{Hg}}_{\mathrm{T}}$ removed from the water column.

The theoretical results from the HR3DHG model show that, without specific and appropriate recovery actions, the mercury benthic flux could remain high for a very long time, representing a threat for this environment, for its ecosystems, and for human health.

Finally, with its features, the HR3DHG model may represent a useful tool to explore and predict the effects of environmental changes on mercury dynamics for several possible forthcoming scenarios.

6 Conclusions

A novel biogeochemical integrated model, HR3DHG, has been designed and implemented to reproduce the spatiotemporal dynamics of three species of mercury in the highly contaminated Augusta Bay. The model consistently reproduces the biogeochemical dynamics of mercury fluxes at the boundaries of the 3D domain, which is necessary for an accurate and reliable approximation of the annual mass balance for the whole basin. Direct comparison of model and experimental data suggests a good capacity of HR3DHG to capture the crucial processes dominating the dynamics of $\mathrm{Hg}$ species in the different marine compartments and at their interfaces, with reliable estimations of benthic fluxes and evasion towards the atmosphere. The model provides robust information on the recycling of the $\mathrm{Hg}$ species in a confined coastal area and can be considered a reliable numerical tool to describe the high-resolution variability of the most important biogeochemical variables driving $\mathrm{Hg}$ concentrations. Finally, model results for the Augusta Bay suggest a permanent and relevant long-term (at century scale) mercury benthic flux associated with negative effects for the biota of the investigated marine ecosystem and with significant health risks.

Code and data availability

The experimental data used in this study are available and properly referenced in the paper or collected in the tables in the Supplement. The software code files are available on GitHub (10.5281/zenodo.3384784, ) and can be also downloaded from http://biomatlab.iasi.cnr.it/models/hr3dhgv1.zip (last access: 20 March 2020).

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-2073-2020-supplement.

Author contributions

GD devised the HR3DHG model. GD, AB, and ADG designed the software used to numerically solve the equations of the model. GD and MS jointly wrote the paper. DV and BS supported the HR3DHG model development. AB, DV, BS, and ADG managed the simulations. DSM and MB performed the Hg data collection. DSM developed the sampling strategy. DSM and MB performed the study of Hg biogeochemistry. MS investigated the Hg biogeochemical dynamics. AC investigated the hydrodynamics in Augusta Bay. AC performed the ocean modeling and generated the code of SHYFEM. EQ performed the data statistics and mapping. All authors contributed to reviewing the paper.

Competing interests

The authors declare that they have no conflict of interest.

Financial support

This research has been supported by the Ministry of University, Research and Education of the Italian government under project “Centro Internazionale di Studi Avanzati su Ambiente, ecosistema e Salute umana – CISAS” (grant no. CUP: B62F15001070005, CIPE n. 105/2015).

Review statement

This paper was edited by Carlos Sierra and reviewed by Ginevra Rosati and Dusan Zagar.