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

The biogeochemical dynamics of Hg, and specifically of its three species Hg⁰, Hg^II, and 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 ([Hg_D]) 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 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 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 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 (Driscoll et al., 2013; Batrakova et al., 2014). Over the last few years some theoretical studies have offered innovative tools to reproduce the mass balance and the dynamics of [Hg] 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 (Ciffroy, 2015) and the WASP (Water Analysis Simulation Program) model (Melaku Canu et al., 2015; Canu and Rosati, 2017; Rosati et al., 2018). In particular, the River MERLIN-Expo model (Ciffroy, 2015) 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 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 (Melaku Canu et al., 2015; Canu and Rosati, 2017; Rosati et al., 2018). In particular, a WASP model applied to a 1D domain and calibrated by using experimental data for dissolved Hg and MeHg allowed us to explore [Hg] dynamics in the Black Sea (Rosati et al., 2018). Similarly, the box model approach has been adopted in a 2D configuration (Melaku Canu et al., 2015; Canu and Rosati, 2017) to calculate Hg mass balance in the coastal areas of the Marano–Grado lagoon (northern Italy), where heterogeneous spatial distributions of 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 (Pakhomova et al., 2018) the biochemistry of 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) (Yakushev et al., 2017) to reproduce the vertical dynamics of the total dissolved Hg and MeHg in the marine coastal areas of the Étang de Berre lagoon (France) (Pakhomova et al., 2018). 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 Soerensen et al. (2016), who devised a high-resolution 1D model for Hg species present in water and sediments of the Baltic Sea (Soerensen et al., 2016). In both HR models, however, the strong impact of the horizontal velocity field on the spatiotemporal distribution of [Hg] could not be considered since 1D modeling was used.

In general, the appropriate modeling to reproduce the spatial and temporal variability of Hg species in highly heterogeneous marine ecosystems, such as Augusta Harbour, requires the use of a hydrodynamics model integrated with a biogeochemical model (Zagar et al., 2007, 2014). For this aim, Zagar et al. (2007) 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 Hg species. Although the grid used did not guarantee a high spatial resolution, the modified PCFLOW3D model allowed for obtaining, for all the Hg species, theoretical vertical profiles of [Hg] in acceptable agreement with experimental data for most of the Mediterranean Sea (Zagar et al., 2007) and to improve the flux estimations of Hg mass balance for the whole Mediterranean basin (Rajar et al., 2007). By following the same modeling approach of Rajar et al. (2007), over the last decade several authors have used 3D advection–diffusion–reaction models to simulate the spatiotemporal dynamics of [Hg_D] and [Hg_T] in oceans, lakes, and estuaries (Zhu et al., 2018). However, the 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, Zhang et al. (2014) reproduced the [Hg_T] in oceans and calculated Hg mass balance by using a 3D ocean tracer model (OFFTRAC) coupled with a general circulation model (GEOS-Chem) (Zhang et al., 2014). Here, the sinking flux of 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 [Hg] 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 Hg species in seawater (Hg⁰, Hg^II, and MeHg) coupled with a diffusion–reaction model for dissolved 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 (Sprovieri et al., 2011; Bagnato et al., 2013; Sprovieri, 2015; Oliveri et al., 2016; Salvagio Manta et al., 2016): 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 Hg at highly contaminated coastal marine sites. The HR3DHG model has been designed to predict the biogeochemical behavior of 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 (La Barbera and Spagnolo, 2002; Fiasconaro et al., 2004; Valenti et al., 2004, 2008; Dutkiewicz et al., 2009; Morozov et al., 2010; Valenti et al., 2012; Denaro et al., 2013a, b, c; Valenti et al., 2015, 2016a, b, c, 2017; Morozov et al., 2019) and the mechanisms responsible for the uptake of Hg within cells (Pickhardt and Fischer, 2007; Radomyski and Ciffroy, 2015; Lee and Fischer, 2017; Williams et al., 2010) are taken into account as specific contributions to the scavenging process and 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 [Hg] within the 3D domain; (ii) to estimate the 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 MeHg and Hg^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 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 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 30 km² 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 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 (Bellucci et al., 2012). Although discharge activities were definitively stopped in 2005, the Hg contamination from the chlor-alkali plant remains a critical environmental threat, with extremely high [Hg] in the bottom sediments (ICRAM, 2008; Sprovieri et al., 2011; Oliveri et al., 2016), significant Hg evasion fluxes from sediments to seawater (Salvagio Manta et al., 2016) and to the atmosphere (Bagnato et al., 2013; Sprovieri, 2015), and evident and recently documented risks for the ecosystem (Tomasello et al., 2012; Bonsignore et al., 2013) and for human health (Bianchi et al., 2006; Bonsignore et al., 2015, 2016). The geographical position, together with its geological and oceanographic features, assigns this area a key role in the Hg inventory at Mediterranean scale. The estimate of the Hg export from Augusta Bay to the open sea (0.54 kmol yr⁻¹; Salvagio Manta et al., 2016) corresponds to about 4 % of total input from coastal point and diffuse sources to the Mediterranean Sea (12.5 kmol yr⁻¹; Rajar et al., 2007). A very narrow shelf develops down to 100–130 m with a mean gradient of about 1.0^∘ and a next steep slope characterized by a dense net of canyons dropping to the deep Ionian basin (Budillon et al., 2008). The Augusta Harbour covers a surface of 23.5 km² with two main inlets connecting with the open sea: the Scirocco (300 m wide and 13 m deep) and the Levante inlets (400 m wide and 40 m deep). The bottom is mainly flat with an average depth of 15 m, with the exception of a deeper channel about 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 10 and 20 cm, and the winds are generally from the northwest and northeast with an average speed of around 3 m s⁻¹ (De Marchis et al., 2014). 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.

3 Model description

The HR3DHG model has been designed and implemented to reconstruct, at high spatiotemporal resolution, the behavior of [Hg_T] and [Hg_D]. 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 Zagar et al. (2007), 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 (Burchard and Petersen, 1999; Umgiesser et al., 2004; Umgiesser, 2009; Umgiesser et al., 2014; Ferrarin et al., 2014; Cucco et al., 2016a, b, 2019), 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 (Denaro et al., 2013a, c; Valenti et al., 2015, 2016a, b, 2017) (see Sect. S4 of the Supplement for details). By using the curve of the mean vertical profile obtained by Brunet et al. (2007) 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 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 (Ciffroy, 2015; Radomyski and Ciffroy, 2015) (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 Hg released by dead picoeukaryote cells in the seawater compartment (see Fig. 2).

Figure 2 Basic structure of the HR3DHG model.

Figure 3 Basic scheme used for implementing the HR3DHG model. The scheme describes the transformation processes (k₁ – photooxidation, k₂ – photoreduction, k₃ – biological oxidation, k₄ – biological reduction, k_Ph-de – photodemethylation, k_deme – demethylation, k_me – methylation, K_demeth – demethylation, K_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 Hg species in seawater (${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{II}}$, MeHg_D, ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{0}}$, ${\mathrm{Hg}}_{\mathrm{P}}^{\mathrm{II}}$, MeHg_P) and sediments (${\mathrm{Hg}}_{\mathrm{pw}}^{\mathrm{II}}$, MeHg_pw, ${\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}$).

3.1 The advection–diffusion–reaction model for the Hg species in seawater

The dynamics of the [Hg_D] 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 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 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 Me₂Hg concentration is very low in the Augusta Harbour (Sprovieri, 2015), the behavior of this 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 (Sprovieri et al., 2011). Specifically, z represents the depth of the barycenter of each sub-domain, localized between the surface (z=0) and the bottom (z=z_b), while x and y indicate the distance in meters measured from a reference point (lat. 37^∘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 10 and 18 elements regularly spaced at 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 (Tveito and Winther, 1998). The stability analysis, performed according to previously published methods (Roache, 1998; Tveito and Winther, 1998; Thi et al., 2005), indicates that the convergence of our algorithm is guaranteed.

The dynamics of the Hg species in seawater are represented through six processes (Zhang et al., 2014; Melaku Canu et al., 2015): (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 Hg⁰ and Hg^II have been described as reaction terms with a first-order kinetic (Batrakova et al., 2014; Zhang et al., 2014; Melaku Canu et al., 2015). 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) (Han et al., 2007; Zhang et al., 2014). 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 (Zhang et al., 2014).

The model includes three reaction terms regulated by first-order kinetics, which describe the photodemethylation of MeHg, the methylation of Hg^II, and the biotic demethylation of MeHg (Batrakova et al., 2014; Melaku Canu et al., 2015). The first is the amount of Hg⁰ produced by MeHg through photochemical reactions. The second is the amount of MeHg obtained by Hg^II through biotic and abiotic pathways in seawater. The third is the amount of Hg^II produced by 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 (Monperrus et al., 2007b, a; Lehnherr et al., 2011; Batrakova et al., 2014; Melaku Canu et al., 2015).

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 Hg_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 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 (Umgiesser et al., 2004; Umgiesser, 2009; Cucco et al., 2016a, 2019) (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 Hg^II and 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 Hg_D species is regulated by the sinking flux of particle-bound mercury along the water column (Zhang et al., 2014), which depends on variables calculated by using the NP model. The amount of Hg released by particulate organic matter is primarily estimated through variables defined in the NP model and Phytoplankton MERLIN-Expo model (Valenti et al., 2012; Denaro et al., 2013a, b, c; Valenti et al., 2015, 2016a, b, c, 2017; Radomyski and Ciffroy, 2015) (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 (Brunet et al., 2007; Baines et al., 1994) (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 Hg within the suspended particulate matter. The Phytoplankton MERLIN-Expo model gives the spatiotemporal dynamics of the Hg^II and MeHg contents within the picoeukaryote cells (Radomyski and Ciffroy, 2015). These two variables are then used, together with the picoeukaryote abundance, to obtain the amount of Hg^II and 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 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₁, k₂, k₃, and k₄ are the rate constants for the photooxidation of Hg⁰, the photoreduction of Hg^II, the biological oxidation of Hg⁰, and the biological reduction of Hg^II, respectively (h⁻¹); k_Ph-de is the rate constant for the photodemethylation of MeHg (h⁻¹); k_deme and k_me are the rate constants for the biotic demethylation of MeHg and the methylation of Hg^II, respectively (h⁻¹); $S_{\mathrm{L}}^{\mathrm{0}}$, $S_{\mathrm{L}}^{\mathrm{II}}$, and $S_{\mathrm{L}}^{\mathrm{MM}}$ are the direct loads for Hg⁰, Hg^II, and 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 MeHg_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 Hg^II and MeHg, respectively $\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{h}}^{-\mathrm{1}}\right)$.

The photochemical rate constants (k₁ and k₂) are directly proportional to the shortwave radiation flux (RAD) at the water–atmosphere interface (Zhang et al., 2014; Soerensen et al., 2010; Qureshi et al., 2010; Batrakova et al., 2014), while the biological rate constants (k₃ and k₄) are calculated by the organic carbon remineralization rate (OCRR) of the microbial reactions (Zhang et al., 2014) (see Sect. S1 of the Supplement). The k_me and k_deme values are fixed according to Lehnherr et al. (2011), while k_Ph-de is set according to Melaku Canu et al. (2015).

The two sinking fluxes ($S_{\mathrm{SPM}}^{\mathrm{II}}$ and $S_{\mathrm{SPM}}^{\mathrm{MM}}$) are obtained according to previous works (Zhang et al., 2014; Rosati et al., 2018), 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 Hg for Hg^II and 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 Hg for Hg^II and 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⁻² h⁻¹); the pe ratio is the ratio of particulate organic carbon (POC) export to NPP out of the euphotic zone (dimensionless); z₀ is the depth of the euphotic zone (m); k_D is the seawater–SPM partition coefficient for Hg_D (L kg⁻¹); f_oc is the fraction of suspended particulate matter as organic carbon (dimensionless); v_silt is the silt settling velocity (m h⁻¹); $k_{\text{D silt}}^{\mathrm{II}}$ is the partition coefficient of Hg^II to silt (L kg⁻¹); $k_{\text{D silt}}^{\mathrm{MM}}$ is the partition coefficient of MeHg to silt (L kg⁻¹); SPIM is the suspended particulate inorganic matter (kg L⁻¹). The NPP is obtained by Baines et al. (1994) using the conversion equation for the chl a concentration (Baines et al., 1994). This is calculated by the picoeukaryote abundance using the conversion curve of Brunet et al. (2007). The pe ratio is calculated by the surface atmospheric temperature measured from remote sensing and the chl a concentration obtained by the NP model (Zhang et al., 2014). The k_D value was measured within the Augusta Harbour during the last oceanographic survey, while the spatial distributions of f_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_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 Rosati et al. (2018).

The loads of ${\mathrm{Hg}}_{\mathrm{D}}^{\mathrm{II}}$ and MeHg_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 PHg^II and PMeHg are the Hg^II content and MeHg content, respectively, in each cell of picoeukaryotes (µg per cell); b is the picoeukaryote abundance (cells per cubic meter); m is the mortality of picoeukaryotes (h⁻¹); λ is the Hg recycling coefficient for picoeukaryotes (dimensionless). The spatiotemporal dynamics of PHg^II and PMeHg are obtained by solving the ODEs of the Phytoplankton MERLIN-Expo model for Hg^II and MeHg, respectively (Radomyski and Ciffroy, 2015) (see Sect. S5 of the Supplement). The spatiotemporal distribution of b is reproduced by using the NP model (Valenti et al., 2017) (see Sect. S4 of the Supplement). The parameter m is set according to Valenti et al. (2017), while λ is fixed as equal to the nutrient recycling coefficient for picoeukaryotes (Valenti et al., 2015, 2017).

The concentrations [Hg_D] and [Hg_T] 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_D is the seawater–SPM partition coefficient for Hg_D (only Hg^II and MeHg) (L kg⁻¹), and SPM is the suspended particulate matter concentration (kg L⁻¹). The partition coefficient k_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 Hg⁰ at the water–atmosphere interface (Bagnato et al., 2013; Zagar et al., 2007); (ii) the lack of Hg⁰ diffusion at the water–sediment interface (Ogrinc et al., 2007); (iii) the wet and dry deposition of Hg^II at the water–atmosphere interface (Rajar et al., 2007; Zagar et al., 2007); (iv) the wet and dry deposition of MeHg at the water–atmosphere interface (Mason et al., 2012); (v) the diffusion of Hg^II and MeHg at the water–sediment interface; (vi) the constant fixed value of [Hg_D] out of Augusta Bay (Ionian Sea) (Horvat et al., 2003); and (vii) the exchange of elemental mercury, Hg^II, and MeHg between the Augusta basin and the Ionian Sea through the two inlets (Salvagio Manta et al., 2016). 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 (Salvagio Manta et al., 2016). 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, Hg_gas-atm is the gaseous elemental mercury (GEM) concentration in the atmosphere (ng m⁻³); Pr is the amount of precipitation (mm); Δt is the exposition time to precipitation (h); Drydep${}_{{\mathrm{Hg}}^{\mathrm{II}}}$ is the atmospheric dry deposition of Hg^II (ng m⁻² h⁻¹); MTC_water-atm is the gas-phase overall mass transfer coefficient (m h⁻¹); H is the Henry's law constant (dimensionless); ${\mathrm{Hg}}^{\mathrm{0}}{|}_{z=\mathrm{0}}$ is the [Hg⁰] at the sea surface (µg m⁻³); ${\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{II}}$ is the [Hg^II] in the atmosphere (ng m⁻³); ${\mathrm{MTC}}_{\text{sed-water}}^{\mathrm{II}}$ is the mass transfer coefficient for Hg^II at the water–sediment interface (m h⁻¹); ${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$ is the [Hg^II] in the pore water of the surface layer (upper 10 cm) of the sediments (µg m⁻³); ${\mathrm{Hg}}^{\mathrm{II}}{|}_{z=z_{\mathrm{b}}}$ is the dissolved [Hg^II] at the deepest layer of the water column (z=z_b) (µg m⁻³); MTC${}_{\text{sed-water}}^{\mathrm{MM}}$ is the mass transfer coefficient for MeHg at the water–sediment interface (m h⁻¹); MeHg_pore water is the [MeHg] in the pore water in the surface layer (upper 10 cm) of the sediments (µg m⁻³); MeHg${|}_{z=z_{\mathrm{b}}}$ is the dissolved [MeHg] in the deepest layer of the water column (z=z_b) (µg m⁻³); ${\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{0}}$, ${\mathrm{Hg}}_{\mathrm{ext}}^{\mathrm{II}}$, and MeHg_ext are the average [Hg⁰], [Hg^II], and [MeHg], respectively, reported from the Ionian Sea (µg m⁻³). The dynamics of the GEM and Hg^II concentrations in the atmosphere (Hg_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 (Bagnato et al., 2013), 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_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 (Covelli et al., 1999; Sørensen et al., 2001; Schulz and Zabel, 2006; Ogrinc et al., 2007; Bryant et al., 2010; Ciffroy, 2015) (see Sect. S1.2.2 and S1.3.2 of the Supplement).

3.1.1 Mass balance of 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. (10)–(16), according to the following equation (Sprovieri et al., 2011; Salvagio Manta et al., 2016):

$$\begin{array}{}\text{(18)}& A+\mathrm{AD}+B=O+D+V,\end{array}$$

where A is the input of Hg_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. (10)–(16), 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 (Sprovieri et al., 2011; Salvagio Manta et al., 2016).

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.

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 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., Hg^II (${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$) and MeHg (MeHg_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^′, y^′, z^′) 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^′ 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 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=0 by interpolating the experimental data collected by ICRAM in 2005 (ICRAM, 2008) (see Sect. S1.2.4 of the Supplement). We then calculated both [Hg^II] and [MeHg] in pore water using the following equations:

$$\begin{array}{}\text{(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}\end{array}$$

where ${{\mathrm{Hg}}_{\mathrm{T}}}^{\mathrm{sed}}\left(\mathrm{0}\right)$ represents the spatial distribution of [Hg_T] in the sediments at an initial time (mg kg⁻¹), f_MeHg is the fraction of MeHg in the sediments (dimensionless), $K_{\mathrm{d}}^{\mathrm{II}}$ is the sediment–pore water distribution coefficient for Hg^II (L kg⁻¹), and $K_{\mathrm{d}}^{\mathrm{MM}}$ is the sediment–pore water distribution coefficient for MeHg (L kg⁻¹).

In pore water, the dynamics of [Hg^II] and [MeHg] are modeled by considering three chemical–physical processes (Schulz and Zabel, 2006; Melaku Canu et al., 2015; Oliveri et al., 2016): (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 Hg^II and 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 (Hines et al., 2012).

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 Hg^II ($D_{\mathrm{sed}}^{\mathrm{in}}(x^{\prime },y^{\prime },z^{\prime })$) and 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 (Schulz and Zabel, 2006; Melaku Canu et al., 2015).

The desorption term estimates the increase in ${{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}$ and MeHg_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_demeth is the rate constant for the demethylation of MeHg (h⁻¹); K_meth is the rate constant for the methylation of Hg^II (h⁻¹); α is the desorption rate for the $\left[{\mathrm{Hg}}_{\mathrm{T}}^{\mathrm{sed}}\right]$ bound to the sediment particles (h⁻¹).

As boundary conditions, we assume a null value of mercury flux at the bottom of the sediment column (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 (20)–(26) represent the three-dimensional diffusion–reaction model used to describe and reproduce the spatiotemporal dynamics of [Hg^II] 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. (13)–(14), (20)–(21), and (19), 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 Hg_D and Hg_T, set to 1.9 ng L⁻¹, corresponding to the experimental detection limit. Specifically, the initial concentration of each Hg species was fixed according to the percentage observed in seawater samples from the Ionian Sea (Horvat et al., 2003) in such a way as to respect the detection limit for total [Hg_D]. The numerical results were not affected by the chosen initial conditions; indeed, the same spatial distribution of [Hg] at nearly steady state was obtained when initial 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 (Horvat et al., 2003; Schulz and Zabel, 2006; Monperrus et al., 2007b, a; Strode et al., 2010; Lehnherr et al., 2011; Melaku Canu et al., 2015; Sprovieri et al., 2011; Salvagio Manta et al., 2016; Sunderland et al., 2006; Zhang et al., 2014; Batrakova et al., 2014), while other parameters, including those most sensitive for the model, were calibrated to correctly reproduce the experimental data collected during the six oceanographic surveys (Sprovieri et al., 2011; ICRAM, 2008; Bagnato et al., 2013; Salvagio Manta et al., 2016; Oliveri et al., 2016). 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 2005 and 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 Bagnato et al. (2013), Salvagio Manta et al. (2016), and Oliveri et al. (2016) for a detailed description of the measured parameters, the related dynamics, and the analytical methods used (ICRAM, 2008; Bagnato et al., 2013; Salvagio Manta et al., 2016; Oliveri et al., 2016). 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. (20)–(21) the sediment–pore water distribution coefficients were calibrated to guarantee the best theoretical [Hg] in pore water in agreement with the value ranges experimentally observed in a previous work (Oliveri et al., 2016), 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 2017 (see Table S1). In Eq. (22), the desorption rate α was calibrated to obtain the best fit between the theoretical results and experimental observations for [Hg] 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 Sørensen et al. (2001), while the spatial distribution of the sediment porosity within Augusta Harbour was reproduced according to previous works (Covelli et al., 1999; Ogrinc et al., 2007) by exploiting the measurements on the sediment samples performed by ICRAM in 2005. Also, the molecular diffusion coefficient was that reported by Schulz and Zabel (2006).

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 (Denman and Gargett, 1983) 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 Massel (1999). 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 (Massel, 1999) 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 (Melaku Canu et al., 2015; Oliveri et al., 2016; Liu et al., 2012; Cossa and Coquery, 2005; Ciffroy, 2015).

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: [MeHg_D], [Hg_D], and [Hg_T] measured in seawater and all annual 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. (1)–(17), together with the equation system (20)–(26), 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 [Hg_D] are obtained by the model for the sites closest to stations 1, 10, 12, and 20 (sampling May 2011), station 15 (sampling June 2012), and station 26 (sampling February 2012).

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 (Zhang et al., 2014; Melaku Canu et al., 2015). 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 Hg^II and 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 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 [Hg] 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 χ² test, indicates good agreement between the model results and experimental findings for [MeHg] at stations A3 (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0005}$) and A7 (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0005}$), while differences can be observed at stations A9 (${\stackrel{\mathrm{̃}}{\mathit{\chi }}}^{\mathrm{2}}=\mathrm{0.0955}$) and A11 (${\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 MeHg benthic fluxes at these two stations.

Figure 5 Distribution of MeHg and Hg_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).

In our analysis, the spatiotemporal behavior of [Hg_D] is obtained as the sum of the three dissolved mercury species. On the other hand, the dynamics of the spatial distribution of [Hg_T] are estimated according to Eq. (9), assuming a linear correlation between the modeled [Hg_D] and the experimental SPM concentrations. The spatial distributions of [Hg_D] and [Hg_T] are reported for May 2011 in Figs. S2–S5 of the Supplement.

The numerical results for [Hg_D] 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 [Hg_D] is less than the experimental error (σ=3.2 ng L⁻¹) in 59 % of sampling points, while it exceeds 2σ in only 17 % of sampling sites. As a conclusion, the comparison between experimental data and theoretical results for [Hg_D] 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 [Hg_T] 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 [Hg_T] is less than σ in 44 % of cases, while it exceeds 2σ at 32 % of sampling sites. The differences (larger than σ=3.2 ng L⁻¹) 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. 9), 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 Hg⁰ 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).

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. 5. Here, very high benthic Hg^II and 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 Hg_D fluxes are mainly generated by the diffusion process at the seawater–sediment interface and that the amount of Hg_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 Hg_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 3 (November 2011) and 5 (June 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 (Bagnato et al., 2013), while a reduced evasion flux is observed at the northern sector of the basin (see Fig. 6).

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

In this work, we calculate the annual mass balance of the Augusta Bay to study the fate of Hg coming from sediments and to estimate the Hg outflows at the inlets of basin. In Fig. 7, 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 MeHg flux appears to be 1 to 2 orders of magnitude lower. The model results are compared with experimental information reported by Salvagio Manta et al. (2016) for three different sampling sites and in two different periods (September 2011 and June 2012). The modeled Hg_D benthic fluxes (2.65 kmol yr⁻¹ for the year 2011 and 2.61 kmol yr⁻¹ 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⁻¹ in September 2011 and 1.4 kmol yr⁻¹ in June 2012) (Salvagio Manta et al., 2016). 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⁻¹ for the year 2011 and $\mathrm{1.90}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr⁻¹ for the year 2012) are in good agreement with those estimated by Salvagio Manta et al. (2016) for each year ($\mathrm{1.70}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr⁻¹). 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⁻¹) and that estimated in the experimental work ($\mathrm{0.42}\times {\mathrm{10}}^{-\mathrm{2}}$ kmol yr⁻¹) (Salvagio Manta et al., 2016). 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 Bagnato et al. (2013) the AD is calculated by averaging the experimental data acquired during a time-limited sampling period (from 29 August 2011 to 23 April 2012), namely without considering the period in which the amount of precipitation is very low. In this way, the AD obtained by Bagnato et al. (2013) 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 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 Hg_T contribution. In Fig. 7c, we show the annual Hg_T outflow from the Augusta Bay towards the open sea. This has been estimated to be 0.13 kmol yr⁻¹ for the year 2012 and appears significantly lower than the 0.51 kmol yr⁻¹ calculated by Salvagio Manta et al. (2016) for the same year. Our hypothesis to explain this discrepancy is that the previous study does not consider the dynamics of [Hg_T] at the inlets (the Hg_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. (18), as was done in previous works on the Augusta Bay (Sprovieri et al., 2011; Salvagio Manta et al., 2016). The model results for the recycled mercury flux are shown in Fig. 7d. Here, the values calculated by our model (2.50 kmol yr⁻¹ for the year 2011 and 2.46 kmol yr⁻¹ for the year 2012) are larger and probably more realistic than those estimated in Salvagio Manta et al. (2016) (0.84 kmol yr⁻¹).   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 (Salvagio Manta et al., 2016).

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⁻¹ for the year 2011) and the sinking flux (0.07 kmol yr⁻¹ for the year 2011) is observed, probably due to the underestimation of the amount of mercury captured by POM (see Eqs. 4–5). More specifically, this behavior could be caused by the underestimation of NPP, which is calculated by using a conversion equation calibrated for oceans (Baines et al., 1994) rather than for coastal zones.

In contrast, very high values of the annual Hg_T accumulation rate in the surface sediment layer (12.07 kmol yr⁻¹ for the year 2011), with respect to those of the annual recycled flux (2.50 kmol yr⁻¹ for the year 2011), are obtained by our model. This result is caused by the high sedimentation rate (11.7 mm yr⁻¹) estimated experimentally (Sprovieri, 2015; ICRAM, 2008) and used in our calculations for the annual Hg_T accumulation rate (Covelli et al., 1999). 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 (Hg^II and MeHg) dissolved in pore water have been obtained by solving Eqs. (20)–(26). 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 [Hg^II] and [MeHg] in pore waters (b, c, e, f) at stations 8 and 16 (sampling May 2011).

In Fig. 8, 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 [MeHg_pore water] always reach their maximum value within the surface layer of the sediments (<0.5 m of depth). However, the shape of the vertical profiles for $\left[{{\mathrm{Hg}}^{\mathrm{II}}}_{\text{pore water}}\right]$ and [MeHg_pore water] 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 [MeHg_pore water] 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 [MeHg_pore water] 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 (Melaku Canu et al., 2015; Yakushev et al., 2017; Pakhomova et al., 2018; Valenti et al., 2017; Dutkiewicz et al., 2009) 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 (Pickhardt and Fischer, 2007; Radomyski and Ciffroy, 2015) 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 (Dutkiewicz et al., 2009; Morozov et al., 2010; Valenti et al., 2012; Denaro et al., 2013a, b, c; Valenti et al., 2015, 2016a, b, c, 2017; Morozov et al., 2019). 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 [Hg] in seawater and sediment, as well as 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 (Soerensen et al., 2016; Pakhomova et al., 2018), 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 (Rajar et al., 2007; Zagar et al., 2007; Canu and Rosati, 2017) 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 (Soerensen et al., 2016; Pakhomova et al., 2018) provided neither the NPP from the Nutrient–Phytoplankton model (Baines et al., 1994; Brunet et al., 2007; Zhang et al., 2014) nor the load of POM-released Hg_D obtained using the Phytoplankton MERLIN-Expo model (Pickhardt and Fischer, 2007; Radomyski and Ciffroy, 2015) (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 [Hg_D] primarily depend on the velocity field of the marine currents obtained from the hydrodynamic model (Umgiesser, 2009; Umgiesser et al., 2014; Cucco et al., 2016a, b), even if the role played by the vertical and horizontal diffusivities (Pacanowski and Philander, 1981; Massel, 1999; Denman and Gargett, 1983) cannot be neglected. Specifically, the spatiotemporal behavior of [Hg_D] changed significantly when alternative velocity fields for the Augusta Bay were used in the biogeochemical module, confirming a feature already observed in previous models (Zagar et al., 2007). Conversely, limited changes in the spatial distribution of [Hg_D] 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 47 % of the Hg_T on average). Because of the exponential decay of [Hg_T] 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 79.0 % of Hg^II, 18.8 % of elemental mercury, and 2.2 % of MeHg, namely values very similar to those observed experimentally in other contaminated sites (Zhang et al., 2014; Melaku Canu et al., 2015). 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 Hg^II and 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 (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⁻¹ 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 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 (Williams et al., 2010; Tomasello et al., 2012; Lee and Fischer, 2017).

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. (4)–(5). This behavior is probably due to the chl a concentration conversion equation of Baines et al. (1994), 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 Hg dynamics, would be necessary from a theoretical point of view to obtain a better estimation of the Hg_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 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 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 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.