Assessing the spatial and temporal variability of methylmercury biogeochemistry and bioaccumulation in the Mediterranean Sea with a coupled 3D model

Human exposure to mercury (Hg) is a cause of concern, due to the biomagnification of the neurotoxic species monomethylmercury (MMHg) in marine ecosystems. Previous research revealed that commercial fish species in the Mediterranean Sea ecosystems are particularly enriched in Hg, due to a combination of physical and ecological factors. Since the fate of Hg depends on the interactions among several biogeochemical and physical drivers, biogeochemical modeling is crucial to support the integration and interpretation of field data. Here, we develop and apply a coupled transport–biogeochemical–metal bioaccumulation numerical model (OGSTM–BFM–Hg) to simulate the biogeochemical cycling of the main Hg species (Hg^II, Hg⁰, MMHg, and DMHg) in seawater, organic detritus, and through the planktonic food web. The model is applied to a 3D domain of the Mediterranean Sea to investigate the spatial and temporal variability of methylmercury (MeHg) distribution and bioaccumulation and major uncertainties in Hg cycling. Model results reproduce the strong vertical and zonal gradients of MeHg concentrations related to primary production consistently with the observations and highlight the role of winter deep convection and summer water stratification in shaping the MeHg vertical distribution, including subsurface MeHg maximum. The modeled bioaccumulation dynamics in plankton food webs are characterized by a high spatial and temporal variability that is driven by plankton phenology and is consistent with available field data of Hg concentrations in plankton, as well as with other indicators, such as bioconcentration factors (BCFs) and trophic magnification factors (TMFs). Model results pointed out that the increment in water temperature linked to a decline of deep convection can cause an increase in water MeHg concentrations with cascading effects on plankton exposure and bioaccumulation.

1 Introduction

The anthropogenic alteration of the mercury (Hg) biogeochemical cycle has led to global enrichment of Hg concentrations in all the environmental compartments (Amos et al., 2015; UNEP, 2019; Zhang et al., 2014a, b) and concerns over human exposure to neurotoxic monomethylmercury (MMHg), which is produced in the ocean and biomagnifies in marine food webs (UNEP, 2019). The ocean has absorbed 50 % of anthropogenic Hg historical emissions, 35 % of which are currently stored in the water and the rest in the sediments (Zhang et al., 2014b): relative to pre-anthropogenic levels, the Hg enrichment is 230 % in surface ocean waters, 25 % in intermediate waters, and 12 % in deep waters (UNEP, 2019).

Models and observations (Lamborg et al., 2014; Zhang et al., 2014b) suggest that the oceanic water column contains 280–370 Mmol of anthropogenic Hg (equivalent to about 70 000 t), a large part of which is recycled within surface and intermediate water due to the biological carbon pump and the microbial loop (e.g., Zhang et al., 2018). Hg is scavenged from surface water by organic particles through adsorption and uptake, transported downward via particles sinking, and released in the dissolved phase following particle degradation. The locally increased microbial activity and availability of inorganic Hg at the depth of particle remineralization are also thought to promote Hg methylation (Cossa et al., 2009; Heimbürger et al., 2010; Sunderland et al., 2009). Concurrent biological and photochemical transformations occur in the water column, driving the interconversions among the main mercury species in seawater (inorganic oxidized Hg, Hg^II; monomethylmercury, MMHg; elemental Hg, Hg⁰; and dimethylmercury, DMHg). Other important processes driving a fast cycling of Hg species in surface water are the atmospheric deposition of oxidized Hg and the gaseous exchange of Hg⁰ (Jiskra et al., 2021). The resulting vertical distribution of Hg species mirrors the local equilibrium among a number of biogeochemical and physical drivers. However, recent comprehensive reviews of the current understanding of Hg fate in the ocean (Bowman et al., 2020; Mason et al., 2012; Sonke et al., 2013) and in the Mediterranean Sea (Cossa et al., 2022) highlighted still high uncertainty in describing and quantifying most of these processes.

Methylmercury (MeHg; defined as the sum of MMHg and DMHg) is more frequently measured than the individual methylated species. In oceanic and Mediterranean areas with high primary production, the MeHg vertical distribution is characterized by low concentrations at the surface due to photodegradation and concentration maxima coincident with the maximum in apparent oxygen utilization, a proxy for heterotrophic activity (Cossa et al., 2009; Sunderland et al., 2009). These observations suggest that most oceanic MeHg is produced in situ. This conclusion is also supported by estimates of MeHg fluxes in the ocean (Mason et al., 2012; Zhang et al., 2020) and by studies based on stable Hg isotopes (Blum et al., 2013; Motta et al., 2019). The role of DMHg in controlling MMHg dynamics remains puzzling: recent observations (Cossa et al., 2017) from the Northwestern Mediterranean Sea support the idea that DMHg is produced from Hg^II in highly productive waters and then degraded to MMHg (Conaway et al., 2009). However, only an abiotic pathway has been directly identified for the formation of DMHg (Jonsson et al., 2016), and rates of DMHg formation from Hg^II measured in polar water were negligible compared to methylation rates from Hg^II to MMHg (Lehnherr et al., 2011).

Tunas and other commercial fish species from the Mediterranean Sea are particularly enriched in Hg in comparison to the same species from other geographical areas (Cossa et al., 2012; Cossa and Coquery, 2005; Harmelin-Vivien et al., 2009; Tseng et al., 2021). It has been suggested that, along with other ecological features of the whole food web, this is due to a shallower occurrence of the MeHg concentration maxima compared to the ocean which results in higher phytoplankton exposure and bioaccumulation (Cossa et al., 2012). Available observations on MeHg distribution in the Mediterranean Sea suggest that concentrations are lower in the oligotrophic waters of the Ionian Sea than in the mesotrophic waters of the Northwestern Mediterranean (Cossa et al., 2009, 2022) and Adriatic Sea (Kotnik et al., 2015).

Recent efforts to couple the dynamics of biological carbon pump and microbial loop with Hg dynamics, including bioaccumulation in the lower food web, provided global-scale assessments of MeHg production and bioaccumulation (Wu et al., 2020, 2021; Zhang et al., 2020). The only study modeling the Hg dynamics in the Mediterranean Sea (Žagar et al., 2007) did not include any biological component. Much of the current knowledge on the dynamics of MeHg in the open sea was acquired only afterwards (Bowman et al., 2016, 2015; Cossa et al., 2017, 2009; Heimbürger et al., 2010; Lehnherr et al., 2011; Monperrus et al., 2007; Munson et al., 2015, 2018; Ortiz et al., 2015; Sunderland et al., 2009).

Here, we couple a model for Hg biogeochemistry (Canu et al., 2019; Rosati et al., 2018, 2020; Zhang et al., 2020) with the 3D transport biogeochemical model OGSTM–BFM (Cossarini et al., 2021; Lazzari et al., 2010, 2021; Salon et al., 2019) to investigate spatial and temporal variability of MeHg dynamics in the Mediterranean Sea. The coupled model, which has a $\mathrm{1}/\mathrm{16}{}^{\circ }$ horizontal resolution, describes the biotic and abiotic interconversions among four Hg species (Hg^II, MMHg, Hg⁰, and DMHg), as well as other key processes of the Hg cycle such as partitioning to organic detritus, sinking to the seabed, exchange of gaseous Hg with the atmosphere, and bioaccumulation in four phytoplankton and four zooplankton functional groups (Lazzari et al., 2012; Vichi et al., 2015). The agreement with observations has been improved by tuning the sinking velocity of organic detritus and by exploring the sensitivity to alternative parameterizations of Hg methylation. A first assessment of the uncertainties concerning the fate of riverine Hg loads is also achieved and discussed along with other terms of the Hg budget in the Mediterranean Sea. Model results are compared to available observations discussing the implications of spatial and temporal variability of modeled MeHg distributions in water and plankton.

2 Methods

2.1 Study area

The Mediterranean Sea (Fig. 1) is a semi-enclosed basin characterized by decreasing west–east gradient of primary production resulting from the superposition of biological pump and inverse estuarine circulation, as well as from the limited impact of riverine loads on the open-ocean dynamics (Crise et al., 1999, Crispi et al., 2001). Due to the formation of intermediate and deep water generating thermohaline circulation, the Mediterranean is considered a miniature of the ocean (Pinardi et al., 2019). The large-scale basin circulation (Pinardi et al., 2015, 2019; Pinardi and Masetti, 2000, and references therein) is composed of three thermohaline cells (Fig. 1). A zonal vertical circulation belt is associated with the inflow of shallow Atlantic Water (AW) at the Strait of Gibraltar that becomes progressively saltier due to evaporation moving eastward (MAW, Modified Atlantic Water) and eventually sinks at intermediate depths forming the Levantine Intermediate Water (LIW), outflowing at Gibraltar. The other cells are driven by dense water formation due to winter cooling, salting, and sinking of surface water. Western Mediterranean Deep Water (WMDW) originates in the Northwestern Mediterranean Sea, while Eastern Mediterranean Deep Water (EMDW) is formed in the Adriatic Sea but also in the Aegean and the Levantine Sea. The mixing between these deep water masses is inhibited by the shallowness of the Strait of Sicily (about 500 m). Long-term observations of ocean color (Bosc et al., 2004; Bricaud et al., 2002; D'Ortenzio and D'Alcalà, 2009) and coupled transport–biogeochemical models (Di Biagio et al., 2019, 2021, Lazzari et al., 2012, 2014, 2021) supported the characterization of different biogeoprovinces based on the seasonal, inter-annual, and high-frequency variability of biogeochemical processes in different areas (subbasins). Biogeochemical modeling also revealed the importance of subsurface plankton blooms, which are not captured by ocean color satellite observations, highlighting the existence of a deep chlorophyll maximum that becomes deeper with increasing water oligotrophy. The most productive areas are the mesotrophic western subbasins of the Alboran Sea (Alb), Northwestern and Southwestern Mediterranean (Nwm and Swm). The eastern Mediterranean is ultra-oligotrophic, apart from the Northern Adriatic Sea (Nad), which is a coastal sea, and other “intermittently blooming” areas (i.e., with high inter-annual variability) such as the Southern Adriatic Sea (Sad), the Aegean Sea (Aeg), and the Rhodes Gyre in the Levantine Sea (Lev).

Figure 1 Bathymetric map of the Mediterranean Sea (from Ocean Data View, Schlitzer, 2014). The white contours indicate the subdivision of the model domain in subbasins: Alboran Sea (Alb), Southwestern Mediterranean (Swm1 and Swm2), Northwestern Mediterranean (Nwm), Tyrrhenian Sea (Tyr1 and Try2), Northern Adriatic Sea (Nad), Southern Adriatic Sea (Sad), Ionian Sea (Ion1, Ion2, and Ion3), Aegean Sea (Aeg), and Levantine Sea (Lev1, Lev2, Lev3, and Lev4). The light blue line shows the zonal cell driven by the surface inflow of Atlantic Water (AW) that moves eastward, forming Modified Atlantic Water (MAW) and Levantine Intermediate Water (LIW), outflowing at intermediate depths at the Strait of Gibraltar. The dark blue lines show the path of meridional cells related to winter convection in the Nwm, Nad, Sad, Aeg, and Lev originating the Western and Eastern Mediterranean Deep Water (WMDW and EMDW).

2.2 The coupled model

The Biogeochemical Flux Model (BFM; Vichi et al., 2015) is a model that simulates the cycling of carbon, nitrogen phosphorus, and silica through dissolved, living, and non-living particulate phases of the marine environment. The description of the planktonic food web is based on the plankton functional type (PFT) approach: the pool of species having similar traits (nutrients affinities, light affinity, prey–predator relationships) are described by a single variable. In the BFM, PFTs have variable stoichiometry and simulate the lower trophic levels of marine food webs up to carnivorous zooplankton (Sect. 2.2.1). The OGSTM–BFM model has been coupled offline to the ocean general circulation model (NEMO) to investigate different issues related to the Mediterranean biogeochemistry (Canu et al., 2015; Cossarini et al., 2015, 2021; Lazzari et al., 2021, 2012, 2014, 2016), and it is now routinely used as the backbone of the biogeochemical component of the Copernicus CMEMS Mediterranean system (https://resources.marine.copernicus.eu/products, last access: 25 July 2022, Salon et al., 2019; Terzić et al., 2019). The integration of a Hg cycling module in the OGSTM–BFM model adds more flexibility to this system and opens the way to novel operational oceanography applications. Here the OGSTM–BFM–Hg model is implemented (Fig. 2) by coupling the OGSTM–BFM model with the equations describing mercury (Hg) species partitioning, transformation, transport (Sect. 2.2.2), bioconcentration in phytoplankton (Sect. 2.2.3), and trophic transfer to zooplankton (Sect. 2.2.4). The coupled model is used to investigate the Hg cycling and bioaccumulation dynamics (Sect. 2.2.5) at the base of the Mediterranean Sea food webs. The model has a $\mathrm{1}/\mathrm{16}{}^{\circ }$ horizontal resolution (approximately 6 km) and 72 vertical levels and runs with a 300 s time step. A model spin-up of 13 years forced with climatological data (Sect. 2.2.6) is performed, and the results of the last year of simulation (2017) are analyzed.

Figure 2 Overview of the coupling in the OGSTM–BFM–Hg model, which integrates the biogeochemistry and bioaccumulation of Hg species in the OGSTM–BFM model, previously coupled to the transport model NEMO and to a multispectral biotical model. The thick grey arrows with text highlight variables and fluxes read by the Hg biogeochemical model (Fig. 3).

2.2.1 BFM model dynamics

The BFM model dynamically simulates nine plankton functional types (PFTs) representative of four phytoplankton groups (picophytoplankton, nanophytoflagellates, diatoms, and large phytoplankton), four zooplankton groups (heterotrophic nanoflagellates, microzooplankton, omnivorous mesozooplankton, and carnivorous mesozooplankton), and one group of heterotrophic bacteria (Vichi et al. 2015). The model reproduces the fluxes of carbon, nitrogen, phosphorous, and silicates from the inorganic nutrient pool to organisms (phytoplankton, zooplankton, and bacteria) and organic compounds pools (particulate and dissolved organic carbon, POC and DOC) with a variable stoichiometric formulation. Non-living organic matter includes one class of organic detritus (POC) and three classes of DOC (labile, semi-labile, and semi-refractory).

2.2.2 Biogeochemical dynamics in the Hg model

The Hg biogeochemical model (Canu and Rosati, 2017; Canu et al., 2019; Melaku Canu et al., 2015; Rosati et al., 2018, 2020, 2022b, Zhang et al., 2020) simulates the cycling of four mercury species in marine water (Fig. 3): inorganic divalent Hg (Hg^II), monomethylmercury (MMHg), elemental mercury (Hg⁰), and dimethylmercury (DMHg). Equations (1)–(4) describe the time variation of non-conservative tracers representing Hg species, omitting advective and diffusive processes that are resolved by the transport model.

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Hg}}^{\mathrm{II}}}{\mathrm{d}t}}& =-k_{\mathrm{met}}\left\{{\mathrm{Hg}}^{\mathrm{II}}\right\}+k_{\mathrm{demet}}\left\{\mathrm{MMHg}\right\}\\ & -(k_{\mathrm{phred}}+k_{\mathrm{biored}})\left\{{\mathrm{Hg}}_{\mathrm{DOC}}^{\mathrm{II}}+{\mathrm{HgCl}}_n\right\}\\ & +\left(k_{\mathrm{phox}}+k_{\mathrm{biox}}\right)\left\{{\mathrm{Hg}}^{\mathrm{0}}\right\}\\ & +k_{\mathrm{phdm}}(\mathrm{1}-\mathit{\alpha })\left\{{\mathrm{MMHg}}_{\mathrm{DOC}}+\mathrm{MeHgCl}\right\}\\ & -{\mathit{\omega }}_{\mathrm{s}}\left\{{\mathrm{Hg}}_{\mathrm{POC}}^{\mathrm{II}}\right\}\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dMMHg}}{\mathrm{d}t}}& =+k_{\mathrm{met}}\left\{{\mathrm{Hg}}^{\mathrm{II}}\right\}{k}_{\mathrm{demet}}\mathrm{MMHg}\\ & -k_{\mathrm{phdm}}\left\{{\mathrm{MMHg}}_{\mathrm{DOC}}+\mathrm{MMHgCl}\right\}\\ & -k_{\mathrm{met}}\left\{\mathrm{MMHg}\right\}+k_{\mathrm{phdm}\mathrm{2}}\mathrm{DMHg}\\ & -{\mathit{\omega }}_{\mathrm{s}}\left\{{\mathrm{MMHg}}_{\mathrm{POC}}\right\}+{\mathrm{Ex}}_{\mathrm{zoo}}{\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{zoo}}}{{\mathrm{C}}_{\mathrm{zoo}}}}\\ & +,D_{\mathrm{zoo}}{\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{zoo}}}{{\mathrm{C}}_{\mathrm{zoo}}}}{\displaystyle \frac{\mathrm{d}{\mathrm{MMHg}}_{\mathrm{phy}}}{\mathrm{d}t}}\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathrm{dHg}}^{\mathrm{0}}}{\mathrm{d}t}}& =+\left(k_{\mathrm{phred}}+k_{\mathrm{biored}}\right)\left\{{\mathrm{Hg}}_{\mathrm{DOC}}^{\mathrm{II}}+{\mathrm{HgCl}}_n\right\}\\ & \left(k_{\mathrm{phox}}+k_{\mathrm{biox}}\right)\left\{{\mathrm{Hg}}^{\mathrm{0}}\right\}+\mathit{\alpha }{k}_{\mathrm{phdm}}\\ & +\left\{{\mathrm{MMHg}}_{\mathrm{DOC}}+\mathrm{MMHgCl}\right\}\\ & -k_{\mathrm{vol}}\left\{{\mathrm{Hg}}_{\mathrm{w}}^{\mathrm{0}}-{\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{0}}}{K_{{\mathrm{H}}^{\prime }}}}\right\}\end{array}\text{(4)}& {\displaystyle }{\displaystyle \frac{\mathrm{dDMHg}}{\mathrm{d}t}}=+k_{\mathrm{met}}\left\{\mathrm{MMHg}\right\}{k}_{\mathrm{phdm}\mathrm{2}}\left\{\mathrm{DMHg}\right\}\end{array}$$

Field observations suggest that particulate Hg dynamics are controlled by regenerative scavenging mechanisms dependent on adsorption, remineralization, and sinking, with little influence of desorption (Lamborg et al., 2016). In the coupled OGSTM–BFM–Hg model, the dynamics POC and DOC including sinking and degradation are spatially resolved by the BFM model, and their concentrations are read by the Hg model (Figs. 2 and 3) to calculate the partitioning of Hg^II and MMHg between particulate species (associated with organic detritus, i.e., Hg${}_{\mathrm{POC}}^{\mathrm{II}}$ and MMHg_POC) and dissolved species (associated with dissolved organic carbon, i.e., Hg${}_{\mathrm{DOC}}^{\mathrm{II}}$ and MMHg_DOC, and ionic species, i.e., HgCl_n and MMHgCl). By combining POC and DOC concentrations with chemical partition constant K_D (Table S1) (Choe et al., 2003; Choe and Gill, 2003; Lamborg et al., 2016), the fraction of each Hg species with respect to its total ($f_{\text{Hg-POC}},f_{\text{Hg-DOC}},f_{\text{HgCl}},f_{\text{MeHg-POC}},f_{\text{MeHg-DOC}},f_{\text{MegCl}}$) is estimated (e.g., Eqs. 5–7 for Hg^II) assuming repartition at thermodynamic equilibrium. The fractions are used to calculate concentrations of particulate and dissolved Hg species from the total pool of Hg^II and MMHg (e.g., Eq. 8). The pool of dissolved gaseous mercury species (DGM, i.e., Hg⁰ and DMHg) is assumed to be entirely in the dissolved phase.

$$\begin{array}{}\text{(5)}& {\displaystyle }{f}_{{\mathrm{Hg}}_{\mathrm{POC}}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{POC}}^{\mathrm{II}}}{{\mathrm{Hg}}_{\mathrm{TOT}}^{\mathrm{II}}}}={\displaystyle \frac{\mathrm{POC}{K}_{D\text{Hg-POC}}}{\mathrm{1}+\mathrm{DOC}{K}_{D\text{Hg-DOC}}+\mathrm{POC}{K}_{D\text{Hg-POC}}}}\text{(6)}& {\displaystyle }{f}_{{\mathrm{Hg}}_{\mathrm{DOC}}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{DOC}}^{\mathrm{II}}}{{\mathrm{Hg}}_{\mathrm{TOT}}^{\mathrm{II}}}}={\displaystyle \frac{\mathrm{DOC}{K}_{D\text{Hg-POC}}}{\mathrm{1}+\mathrm{DOC}{K}_{D\text{Hg-DOC}}+\mathrm{POC}{K}_{D\text{Hg-POC}}}}\text{(7)}& {\displaystyle }{f}_{\mathrm{HgCl}}={\displaystyle \frac{{\mathrm{HgCl}}_n}{{\mathrm{Hg}}_{\mathrm{TOT}}^{\mathrm{II}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{DOC}{K}_{D\text{Hg-DOC}}+\mathrm{POC}{K}_{D\text{Hg-POC}}}}\text{(8)}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{POC}}^{\mathrm{II}}=f_{{\mathrm{Hg}}_{\mathrm{POC}}}{\mathrm{Hg}}_{\mathrm{TOT}}^{\mathrm{II}}\end{array}$$

Once the partitioning is determined, the model computes the transformation fluxes based on first-order rates (k_x, d⁻¹ Table S2), the bioconcentration (Sect. 2.3) and bioaccumulation fluxes (Sect. 2.4), the sinking of Hg_POC and MMHg_POC, the advective transport within the model domain and at the model boundaries (i.e., Strait of Gibraltar), and the exchange of Hg⁰ with the atmosphere. Microbial Hg methylation (k_met Hg^II), reduction ($k_{\mathrm{biored}}\left\{{\mathrm{Hg}}_{\mathrm{DOC}}^{\mathrm{II}}+{\mathrm{HgCl}}_n\right\}$), and oxidation (k_biox Hg⁰) are parameterized as functions of particulate organic matter remineralization (Zhang et al., 2020, 2014), which is read from the BFM model. The direct conversion of Hg^II to DMHg is assumed to be negligible, due to the low rates observed in the only dataset available for this process (Lehnherr et al., 2011). Although several observations linked Hg methylation to oxygen consumption and organic matter remineralization in the water column, a full mechanistic understanding of methylation and demethylation is still lacking (An et al., 2019; Wang et al., 2020). Here we assumed that both reactions involve the entire Hg^II and MMHg pools (Ortiz et al., 2015; Schaefer and Morel, 2009; Zhang et al., 2019). Photochemical transformations, i.e., photoreduction ($k_{\mathrm{phred}}\left\{{\mathrm{Hg}}_{\mathrm{DOC}}^{\mathrm{II}}+{\mathrm{HgCl}}_n\right\}$), photooxidation (k_phox Hg⁰), and photodegradation of MMHg (k_phdem {MMHg_DOC +MMHgCl}) and DMHg (k_phdm2 DMHg) are assumed to act only on the dissolved fractions and are a function of the attenuated shortwave flux (Zhang et al., 2014). MMHg photodegradation is assumed to yield both Hg⁰ and Hg^II (k_phdm {MMHg_DOC+MeHgCl}) due to the contrasting results in experimental studies (Luo et al., 2020).

Figure 3 Hg dynamics in the coupled model OGSTM–BFM–Hg. The Hg model simulates the cycling of inorganic (Hg^II and Hg⁰) and methylated (MMHg and DMHg) Hg species resolving the partitioning of Hg^II and MMHg between particulate (Hg_POC, MMHg_POC) and dissolved species (Hg_Diss, as Hg_DOC+ HgCln, and MMHg_Diss, as MMHg_DOC + MMHgCl) based on partition coefficients (K_D) and concentrations of POC and DOC, the photochemical and biological transformations (dotted colored arrows) as first-order kinetics depending on the rate constants k_x, the evasion of Hg⁰ to the atmosphere controlled by wind speed through the volatilization rate constant k_vol, and the MMHg bioaccumulation (pink arrows) driven by phytoplankton uptake in four different PFTs and trophic transfer to four zooplankton PFTs.

The gas-exchange flux of Hg⁰ with the atmosphere ($k_{\mathrm{vol}}\left\{{\mathrm{Hg}}_{\mathrm{w}}^{\mathrm{0}}-{\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{0}}/K_{{\mathrm{H}}^{\prime }}\right\}$) uses Nightingale's parameterization (Nightingale et al., 2000) for the rate constant k_vol (Eq. 9), which depends on wind speed (u₁₀, m s⁻¹) and on the ratio of Hg Schmidt number to CO₂ Schmidt number ($\mathit{\text{Sc}}{\mathrm{Hg}}_{\mathrm{w}}/\mathit{\text{Sc}}{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{w}}$).

$$\begin{array}{}\text{(9)}& k_{\mathrm{vol}}=\left(\mathrm{0.222}{u}_{\mathrm{10}}^{\mathrm{2}}+\mathrm{0.333}{u}_{\mathrm{10}}\right){\left({\displaystyle \frac{\mathit{\text{Sc}}{\mathrm{Hg}}_{\mathrm{w}}}{\mathit{\text{Sc}}{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{w}}}}\right)}^{-\mathrm{0.5}}\end{array}$$

Current knowledge on dynamics and concentrations of DMHg in the ocean and in the atmosphere does not allow a detailed modeling representation (Baya et al., 2015; Coale et al., 2018; De Simone et al., 2014; Melaku Canu et al., 2015; Nerentorp Mastromonaco et al., 2017; Rosati et al., 2018; Zhang et al., 2020), and simulated variations of DMHg (Eq. 4) depend on photodegradation, production upon MMHg methylation, and advective transport. MMHg is taken up from the water pool by phytoplankton (MMHg_phy) (Sect. 2.3), and it is released from the biotic pool back to the water column by zooplankton excretion (${\mathrm{Ex}}_{\mathrm{zoo}}\frac{{\mathrm{MMHg}}_{\mathrm{zoo}}}{{\mathrm{C}}_{\mathrm{zoo}}}$) and death ($D_{\mathrm{zoo}}\frac{\mathrm{zooMMHg}}{{\mathrm{C}}_{\mathrm{zoo}}}$) (Sect. 2.4). Loadings and initial and boundary conditions are presented in Sect. 2.6.

2.2.3 MMHg uptake by phytoplankton in the Hg model

The bioaccumulation model simulates the MMHg uptake by the phytoplankton PFTs of the BFM model assuming that the uptake flux (Eq. 10, nmol_(Hg) m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$ d⁻¹) decreases at high DOC concentrations and increases with phytoplankton surface area : volume ratios (R_{SV_PFT}), following the empirical relation estimated by Schartup et al. (2018) from multi-species phytoplankton uptake data (Lee et al., 2017; Lee and Fisher, 2016). The bioaccumulation model also considers the phytoplankton density (ρ = 10⁻¹² g${}_{(\mathrm{w}.\mathrm{w}.)}$ µm${}_{(\mathrm{w}.\mathrm{w}.)}^{-\mathrm{3}}$), wet weight biomass ($B_{\mathrm{w}.\mathrm{w}.,\mathrm{PFT}}$, in g${}_{(\mathrm{w}.\mathrm{w}.)}$ m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$), and the water MMHg concentrations ({MMHg}_w, pM). Parameter values are given in Table S3.

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{MMHg}}_{\mathrm{phy},\mathrm{PFT}}}{\mathrm{d}t}}=& \mathrm{0.118}\cdot R_{\mathrm{SV}\mathrm{\_}\mathrm{PFT}}\cdot {\exp }^{(-\mathrm{0.08}\cdot \mathrm{DOC})}\\ & {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{B}_{\mathrm{w}.\mathrm{w}.,\mathrm{PFT}}{\displaystyle \frac{{\left\{\mathrm{MMHg}\right\}}_{\mathrm{w}}}{{\mathrm{10}}^{\mathrm{3}}}}\end{array}\end{array}$$

$B_{\mathrm{w}.\mathrm{w}.,\mathrm{PFT}}$ values were derived from the plankton biomasses expressed on a intracellular carbon content basis (C_PFT, in mg_(C) m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$) by adopting carbon to wet weight conversion factors (Table S4) that we estimated from a dataset of experimental data for various phytoplankton and zooplankton species (Jørgensen et al., 1979; Mahlmann et al., 2008).

2.2.4 Zooplankton bioaccumulation in the Hg model

According to previous studies (Schartup et al., 2018; Zhang et al., 2020), MMHg bioaccumulation in zooplankton is mostly driven by grazing (> 80 %); therefore, the trophic transfer of MMHg from phytoplankton to zooplankton PFTs (heterotrophic nanoflagellates, microzooplankton, omnivorous mesozooplankton, and carnivorous mesozooplankton) in the model is set proportionally to the carbon (C) transfer, neglecting direct uptake from seawater. Thus, the bioaccumulation of MMHg in each zooplankton PFT (Eq. 11) depends on carbon grazing fluxes (G_zoo, PFT, mg_(C) m⁻³ d⁻¹) and the MMHg : C ratio of the preys that are grazed ($\frac{{\mathrm{MMHg}}_{\mathrm{phy},\mathrm{PFT}}}{{\mathrm{C}}_{\mathrm{phy},\mathrm{PFT}}}$, nmol_(Hg) mg${}_{\left(\mathrm{C}\right)}^{-\mathrm{1}}$). MMHg release to the water pool occurs following carbon excretion (Ex_zoo, PFT) and zooplankton death (D_zoo, PFT, mg_(C) m⁻³) and depends on the variable MMHg : C quota ($\frac{{\mathrm{MMHg}}_{\mathrm{zoo},\mathrm{PFT}}}{{\mathrm{C}}_{\mathrm{zoo},\mathrm{PFT}}}$, nmol_(Hg) mg${}_{\left(\mathrm{C}\right)}^{-\mathrm{1}}$) in each zooplankton group.

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{MMHg}}_{\mathrm{zoo},\mathrm{PFT}}}{\mathrm{d}t}}& =G_{\mathrm{zoo},\mathrm{PFT}}{\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{phy},\mathrm{PFT}}}{{\mathrm{C}}_{\mathrm{phy},\mathrm{PFT}}}}\\ & -{\mathrm{Ex}}_{\mathrm{zoo},\mathrm{PFT}}{\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{zoo},\mathrm{PFT}}}{{\mathrm{C}}_{\mathrm{zoo},\mathrm{PFT}}}}\\ & -D_{\mathrm{zoo},\mathrm{PFT}}{\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{zoo},\mathrm{PFT}}}{{\mathrm{C}}_{\mathrm{zoo},\mathrm{PFT}}}}\end{array}\end{array}$$

2.2.5 Indicators for MMHg bioaccumulation and biomagnification

To compare the model content of MMHg in phytoplankton and zooplankton against experimental data, the output in nmol_(Hg)m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$ was converted into ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ by computing the quota of MMHg (Q_MMHg, PFT, Eq. 12) to plankton wet weight ($B_{\mathrm{w}.\mathrm{w}.,\mathrm{PFT}}$, g${}_{\mathrm{w}.\mathrm{w}.}$ m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$) and correcting for the molar mass of Hg (200.59 g mol⁻¹).

$$\begin{array}{}\text{(12)}& Q_{\mathrm{MMHg},\mathrm{PFT}}={\displaystyle \frac{{\mathrm{MMHg}}_{\mathrm{phy},\mathrm{PFT}}}{B_{\mathrm{w}.\mathrm{w}.,\mathrm{PFT}}}}\mathrm{200.59}\end{array}$$

The bioconcentration factor (BCF, L kg⁻¹) is used to synthesize the enrichment of Hg in plankton with respect to water concentrations (Harding et al., 2018; McGeer et al., 2003), and due to the normalization for MMHg water levels, it allows better intercomparability across ecosystems than Hg concentrations. We calculated BCFs (Eq. 13) for each plankton PFT as the ratio between the quota Q_MMHg, PFT and water MMHg concentrations, using model output at monthly resolution spatially averaged. The logarithmic scale is often used to handle BCF values (e.g., McGeer et al., 2003).

$$\begin{array}{}\text{(13)}& \mathrm{BCF}={\mathrm{10}}^{\mathrm{6}}{\displaystyle \frac{Q_{\mathrm{MMHg},\mathrm{PFT}}}{{\mathrm{MMHg}}^{\mathrm{w}}\mathrm{200.59}}}\end{array}$$

The trophic magnification factor (TMF) is used to explore biomagnification across ecosystems by expressing the increase in Hg along with trophic levels (Alava et al., 2018; Harding et al., 2018). TMF is usually calculated in field studies of the whole food web assessing both concentrations of pollutants and trophic levels through nitrogen isotopes. Here, in agreement with other modeling studies (Wu et al., 2021, 2020; Zhang et al., 2020), TMFs were estimated from the ratio of the MMHg quota in zooplankton PFTs to that of grazed phytoplankton PFTs (Eq. 14). The analysis was carried out considering the model plankton food web as composed of a lower trophic level “herbivorous” part and an upper “omnivorous” part.

$$\begin{array}{}\text{(14)}& \mathrm{TMF}={\displaystyle \frac{Q_{\text{MMHg,\hspace{0.17em}zoo-PFTs}}}{Q_{\text{MMHg,\hspace{0.17em} phy-PFTs}}}}\end{array}$$

2.2.6 Initial conditions, boundary conditions, and forcings

Model initial conditions, boundary conditions, and physical forcings at monthly resolution are set according to Lazzari et al. (2021). Initial conditions for Hg species (Table S5) in Mediterranean subbasins (Fig. 1) are extrapolated from published data from different investigations in the Mediterranean Sea (Cossa et al., 2017; Cossa and Coquery, 2005; Ferrara et al., 2003; Horvat et al., 2003; Kotnik et al., 2007, 2015), and boundary conditions are set in agreement with the observations from the GEOTRACES GA03 meridional transect (Bowman et al., 2015). River Hg load is estimated based on dissolved Hg concentration, assuming 3 pM and 3.5 % of MMHg (Cossa et al., 2017), using a higher concentration for three rivers that are known to be highly impacted by mining (Isonzo-Soča river, 25 pM) (Hines et al., 2000), industrial activities (Po river, 6 pM) (Vignati et al., 2008), or both of them (Tiber river, 6 pM) (Lanzillotta et al., 2002; Rimondi et al., 2019), assuming MMHg is the 1 % (Hines et al., 2000). Further, a sensitivity simulation accounting for the load of particulate Hg was carried out (Sect. 2.3.3) to investigate model uncertainty regarding this source. The atmospheric Hg⁰ concentration is set to a constant value ${\mathrm{Hg}}_{\mathrm{atm}}^{\mathrm{0}}=\mathrm{1.6}$ ng m⁻³ (Andersson et al., 2007; Fantozzi et al., 2013; Gårdfeldt et al., 2003), and the atmospheric deposition of Hg_T is set to 67.5 Mg yr⁻¹, 2 % of which is MMHg (Cossa et al., 2022, 2017). The Hg concentrations in the Black Sea outflow at the Dardanelles Strait are set at 0.2 pM in agreement with the observations from the GEOTRACES GA04 cruise and modeled fluxes (Rosati et al., 2018).

2.3 Model sensitivity to key parameters and processes

A set of sensitivity simulations was implemented to investigate model uncertainty and improve the fit against experimental observations, initially showing substantial underestimation of MeHg concentration maxima. The sensitivity simulations aimed at assessing the impacts of variations of the values of the sinking velocity of organic detritus and associated Hg species (Sect. 2.3.1) and of the coefficient for Hg methylation rates (Sect. 2.3.2), as well as the loading of Hg from rivers (Sect. 2.3.3), which is a process poorly characterized by field observations (Liu et al., 2021). All the analyses of model results focus on the same simulation, which adopts a sinking velocity of 10 m d⁻¹, a methylation rate constant higher than the one used in global modeling (Zhang et al., 2020), and accounts only for dissolved input of Hg from rivers.

2.3.1 Sensitivity simulations for POC sinking velocity

A sensitivity analysis of POC sinking velocity was carried out, hypothesizing that a change in this parameter would impact the vertical distribution of all Hg species due to a change in the distribution and remineralization of POC and in the transport of particulate Hg and MeHg. We also speculate that a deeper occurrence of POC remineralization, leading to deeper MeHg maximum, could result in a decreased amount of MeHg photochemically degraded and thus in higher MeHg concentrations. The current version of the BFM model includes only one class of non-living organic detritus (POC) that is a sink for all the plankton state variables (Vichi et al., 2015), resulting in a homogenous POC sinking velocity representative of an ecosystem-averaged dynamic. In the field, organic particles and aggregates that originate from a continuous size spectrum of planktonic cells of various composition show a wide range of sinking speeds (Cael et al., 2021). The reference model simulation was run by adopting the default POC sinking velocity of 3 m d⁻¹ previously used to model the biogeochemistry of the Mediterranean Sea (Lazzari et al., 2021, 2012). Here, two additional sensitivity simulations were run by setting the POC sinking velocity at 10 and 20 m d⁻¹, assuming community-averaged cell diameters of 35 and 50 µm, respectively.

2.3.2 Sensitivity simulations for Hg methylation

Hg methylation was parameterized in the model by scaling the flux of POC remineralization for a constant x_met (Table S2), adopting the same value used in a global ocean model (Zhang et al., 2020), due to the unavailability of more specific data. We explored the sensitivity of modeled MeHg concentrations to a 3-fold increase in x_met, as well as to an alternative parameterization that included both POC and DOC remineralization flux (Eq. S10).

2.3.3 Sensitivity simulation for river Hg load

The largest Mediterranean rivers are characterized by proximal-accumulation-dominated dispersal system, i.e., deltaic systems with fast and substantial (> 50 %–90 %) sediment accumulation in the proximity of their mouth (Walsh and Nittrouer, 2009). Consistently with this picture, a recent budget for the Mediterranean Sea estimated that the magnitude of total Hg input from rivers (6 Mg yr⁻¹) is comparable to the Hg flux to the shelf sediments (6.8 Mg yr⁻¹), suggesting that most of the Hg associated with riverine particles settles before reaching the open ocean (Cossa et al., 2022). We, therefore, included only the dissolved load of Hg and MMHg in the reference simulation and explored the uncertainty related to this choice with a sensitivity simulation accounting for the total load Hg species, i.e., including both the particulate and dissolved inputs of Hg^II and MMHg. Limited data are available to characterize the spatial and temporal variability of these sources. The most robust observations are from a time series (2008–2010) for the Rhône River (Cossa et al., 2017) indicating 0.85 ± 0.45 nmol g⁻¹ of Hg_T and 0.017 ± 0.008 nmol g⁻¹ of MMHg. Assessments in the Ebro River, impacted by an industrially polluted water reservoir, and in the Po River, which drain the largest Italian industrial and agricultural areas, found respectively 4.9 and 1.15 nmol g⁻¹ of particulate Hg_T at the river mouth, also highlighting the role of episodic stormy events in transporting large amounts of particles and associated pollutants (Palanques et al., 2020; Vignati et al., 2003). Further, at the mouth of the Soča/Isonzo River, impacted by upstream Hg mining, concentrations of unfiltered Hg higher than 200 pM have been reported (Hines et al., 2000). Hg loads for the Mediterranean tributary were estimated based on this information as synthesized in Table S6, also qualitatively referring to the spatial distribution of Hg stock and erosion fluxes in European topsoils (Panagos et al., 2021) and other information on anthropogenic activities (Ahmed et al., 2018). In the reference simulation, including only the dissolved input, the loading estimate is 0.28 Mg yr⁻¹ for total dissolved Hg (and 0.008 Mg yr⁻¹ for dissolved MMHg), while in the simulation with total Hg inputs the loading is about 4.6 Mg yr⁻¹, including 0.27 Mg yr⁻¹ of MMHg. Despite the relatively high mean concentrations (Table S6) assumed, the estimated Hg_T load is slightly lower than the estimate (6 Mg yr⁻¹) from Cossa et al. (2022), most likely because the water discharge used in that work (> 350 km³ yr⁻¹) is higher than the total discharge used in this model setup (> 291 km³ yr⁻¹). On the other hand, the load to the Mediterranean Sea from solely European rivers (2.9 Mg yr⁻¹) estimated from Hg concentrations in topsoils and water erosion rates from river catchments (Panagos et al., 2021) is slightly lower than the load from European rivers in our model setup (3.25 Mg yr⁻¹). Finally, all these estimates are lower than the 18 Mg yr⁻¹ computed by Liu et al. (2021) for the Mediterranean Sea through a global river model explicitly considering water discharge, suspended solids, and Hg loadings. Based on these calculations, Liu et al. (2021) proposed a global river load (1000 Mg yr⁻¹) much higher than the estimate from the 2018 global mercury assessment (300 Mg yr⁻¹) (UNEP, 2019). They also estimated that 70 % of the load is buried in coastal and shelf sediment.

3 Results and discussion

3.1 Model sensitivity and uncertainty

The sensitivity analysis showed that the sinking rate of organic detritus (POC) is an important control on the vertical distribution of MeHg (i.e., the shape of the vertical profile) along the water column but has little effect on the values of MeHg concentration maxima and the distribution of inorganic Hg species (Supplement Sect. S1.1 and Figs. S1 and S2). On the other hand, the increase in sinking velocity from 3 to 10 m d⁻¹ improved the agreement of modeled seasonal fluxes of POC at 200 m depth with available observations (Ramondenc et al., 2016) (Supplement Sect. S1.1 and Fig. S3). The sensitivity analysis also showed that an increase in the coefficient for Hg methylation (x_met) results in a higher agreement with observed MeHg concentrations (Cossa et al., 2022, 2009), while the inclusion of DOC remineralization flux results in an overestimation of surficial MeHg concentrations (Supplement Sect. S1.2 and Fig. S3).

The inclusion of higher riverine Hg^II, MMHg, and POC inputs does not improve the agreement between modeled and observed concentrations of MeHg (Supplement Sect. S1.3 and Figs. S5–S10), suggesting that the underestimation of MeHg maxima in subsurface waters is not due to an underestimation of watershed sources. These results also corroborate the idea that a substantial fraction of particulate Hg and MeHg carried by rivers settles in coastal areas along with POC (Fig. S5), as estimated in a recent budget based on field observations (Cossa et al, 2022). Nonetheless, large uncertainties in the spatial and temporal variability of riverine loadings and coastal processes remain to be addressed both at the basin and global scale (Cossarini et al., 2021; Liu et al., 2021).

The underestimation of MeHg in the model does not appear to be related to inorganic Hg availability, since the concentrations of Hg_T are within the experimental uncertainty of the data (Fig. S11).

Modeled fluxes of Hg_T and Hg⁰ compare well with other budget estimates (Cossa et al., 2022; Žagar et al., 2007). The net exchange of gaseous Hg⁰ with the atmosphere indicates 44.8 Mg yr⁻¹, which is close to a previous estimate (50 Mg yr⁻¹) based on aquatic biogeochemical modeling (Žagar et al., 2007) and slightly lower than the flux adopted in the budget from Cossa and coauthors (67.5 Mg yr⁻¹). Considering both Hg⁰ exchange and atmospheric deposition of Hg^II (Sect. 2.2.6 and Fig. 3), the net exchange of Hg species (∼ 7 Mg yr⁻¹) is here significantly lower than the estimate (29.8 Mg yr⁻¹) from Cossa et al. (2022). However, while the latter estimate is based on an atmospheric model assuming constant ocean surface Hg⁰ concentrations (Gencarelli et al., 2014), here, the flux is estimated accounting for 3D dynamics in the marine system, including sequestration of Hg species from surface waters due to the biological carbon pump. Moreover, recent evidence from isotopes studies pointed out that the ocean might be a net sink for atmospheric Hg rather than a net source (Jiskra et al., 2021).

The modeled export of Hg_T from the Alboran Sea to the Atlantic Ocean is 9.3 Mg yr⁻¹, close to the value (9.7 Mg yr⁻¹) estimated by Cossa et al. (2022). On the other hand, our estimate (2.5 Mg yr⁻¹) of MeHg outflow is significantly lower than the 5.6 Mg yr⁻¹ of the budget by Cossa et al. (2022), even assuming an increased value for x_met and, consequently, of subsurface MeHg.

All in all, a calibration of the parameter x_met appears appropriate based on the current understanding, also considering that the initial guess for this parameter comes from a global ocean application (Zhang et al., 2020). Indeed, such a global model study underestimated the observed MeHg concentrations in the Mediterranean Sea, while reproducing with good agreement, or overestimating, observations from various cruises in the North Atlantic and Pacific Ocean, pointing out the need to better resolve spatial and temporal variability of biological processes and linked Hg dynamics (Zhang et al., 2020). A comprehensive validation of the OGSTM–BFM biogeochemical variables (Cossarini et al., 2021) highlighted good model accuracy in surface water but higher mismatches with observed oxygen and nutrient concentrations in the mesopelagic zone, attributed to uncertainties in remineralization and deep ventilation processes, which also affect methylation dynamics. In spite of uncertainties remaining to be addressed, the OGSTM–BFM–Hg model is able to reproduce the spatial gradients observed in the Mediterranean Sea, showing the improvements in capability to simulate Hg dynamics (Žagar et al., 2007, 2014) through the coupling with key biological processes, e.g., the carbon pump and microbial loop (Bowman et al., 2020; Cossa et al., 2017, 2009; Heimbürger et al., 2010; Monperrus et al., 2007; Munson et al., 2018; Sunderland et al., 2009).

3.2 MeHg distribution in the water of the Mediterranean Sea

The zonal and vertical gradients of MeHg concentrations reproduced by the OGSTM–BFM–Hg model show the highest concentrations in subsurface waters of the most productive subbasins of the Mediterranean Sea (Fig. 4), consistently with the available observations (Cossa et al., 2009). Modeled vertical profiles of MeHg concentrations (Fig. 5) are in good agreement with the observations in the Southern Adriatic Sea (Sad, Fig. 5a) and are in the lower range of the observations for the Northwestern Mediterranean (Nwm, Fig. 5b) and Ionian (Ion, Fig. 5c) subbasins. Seasonality affects the modeled distribution of MeHg (Figs. 4 and 5) with more marked effects in the Sad subbasin (Fig. 5a). The strong winter convection occurring in the Sad in 2017 caused remixing in the water column (Mihanović et al., 2021), leading to the disruption of the subsurface MeHg maximum from January to March and an increase in MeHg surface water concentrations, in spite of net demethylation prevailing at all depths during winter months (Fig. 6). A slighter increase in surface water MeHg concentrations during winter months is visible also for the other subbasins (Figs. 4, 5b, and c), but, in the absence of strong winter convection the variation is small, and the model predicts that the subsurface MeHg maximum is a permanent feature throughout the year. The Nwm is known to be an area of winter convection; however, other authors reported a recent decline of deep convection for this area that was weak in 2014 and never occurred in the period 2015–2017 (Margirier et al., 2020). Based on these model results, we speculate that in the long run the reduction of winter mixing can cause MeHg accumulation in the biologically active zone enhancing the exposure of marine organisms.

Figure 4 Spatial distribution of monthly averaged MeHg concentrations (pM) in the Mediterranean Sea at different depths in the water column (25, 90, 175, 270, and 710 m depth) simulated with the OGSTM–BFM–Hg model.

In the Sad subbasin, a progressive buildup of MeHg subsurface maxima (up to 0.18 pM) is predicted from April to September, driven by an increase in primary production triggering higher POC remineralization and in turn higher Hg methylation rate constants k_met (Fig. S4). However, while k_met values are maxima in April and net methylation fluxes (Fig. 6) are maxima in June (0.45 pmol m⁻³ d⁻¹), MeHg concentrations are maxima in September and remain stable until October (Fig. 5), highlighting the importance of water stratification for the formation of the subsurface maxima. Figure 6 also shows that net demethylation prevails throughout the year in deep waters (except for weak positive methylation in Sad and Tyr during summer) and for most of the year in surface waters, suggesting that the presence of MeHg in surface and deep water depends on its formation in the intermediate waters and diapycnal mixing.

Modeled vertical profiles of MeHg in the oligotrophic Ion subbasin (Fig. 5c) have lower maxima (up to 0.16 pM) and a smoother concentration gradient compared to the more productive waters of the Nwm (up to 0.24 pM). The west–east gradient is sustained by the highest methylation rates (Fig. S12) in the western Mediterranean, driven by higher primary productivity and remineralization of organic detritus, and it is further reinforced by high rates of both photochemical and biological MeHg degradation in the eastern subbasins (Fig. S12), due to the highest water temperature, irradiance, and deeper penetration of shortwave radiation caused by the low plankton biomass.

Figure 5 Vertical profiles of monthly averaged MeHg (MMHg + DMHg) concentrations in the water column reproduced by the OGSTM–BFM–Hg model for 2017 (colored solid lines), compared with observations for 2004–2005 (dotted lines with crosses) available from the literature (Cossa et al., 2009) in three subbasins (Fig. 1): (a) Southern Adriatic Sea (Sad), (b) Northwestern Mediterranean (Nwm), and (c) Ionian Sea (Ion2). The right panels show the entire water column, and the left panels show a detail of modeled concentrations in the top 600 m of the same subbasins.

Figure 6 Monthly evolution of the net methylation fluxes (pmol m⁻³ d⁻¹) for each subbasin (colored markers) of the Mediterranean Sea (black lines) depth-integrated for (a) surface water (0–100 m depth), (b) intermediate water (100–600 m depth), and (c) deep water (> 600 m depth). Net methylation is calculated as the difference between modeled Hg methylation fluxes and demethylation fluxes (both dark and photochemical).

3.3 Modeled MMHg in phytoplankton and zooplankton

The spatial and temporal distributions of modeled MMHg in phytoplankton (MMHg_phy) and zooplankton (MMHg_zoo) of the Mediterranean Sea (Fig. 7) are overall comparable to the concentrations of particulate MMHg (0.7 ± 1.3 pmol m⁻³) measured in the Atlantic Ocean (Bowman et al., 2015). The model reproduced a seasonal increase in MMHg_phy (Fig. 7a) during late spring and summer that is driven by blooms of picophytoplankton (Fig. S13) combined with an increase in MMHg availability in surface water (Fig. S14). The amount of MMHg bioaccumulated by picophytoplankton (MMHg_phy,P3, Fig. S15) is high (up to 0.6 pmol_(Hg) m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$) compared to the bioaccumulation by other phytoplankton PFTs (up to 0.02 pmol_(Hg) m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$, Figs. S16–S18) due to the model assumption of increasing bioconcentration capacity with decreasing phytoplankton cell size (Sect. 2.3). However, not all the picophytoplankton blooms result in high values of MMHg_phy: the strongest blooms of picophytoplankton, occurring in March and April at 25 m depth (Fig. S13), have a marginal impact on bioaccumulation, especially in the easternmost subbasins, due to the very low concentrations of MMHg in the most surficial waters (Fig. S14). On the other hand, much weaker blooms of picophytoplankton occurring from May to August slightly deeper in the water column (35 and 45 m depth) lead to a maximum in phytoplankton MMHg concentrations (Fig. 7a) because of higher water MMHg concentrations at these depths during these months (Fig. S14).

Figure 7 Spatial distribution of monthly averaged (a) MMHg_phy and (b) MMHg_zoo (pmol_(Hg) m${}_{\left(\mathrm{w}\right)}^{-\mathrm{3}}$) in the Mediterranean Sea at different water depths (25, 35, and 45 m) simulated with the OGSTM–BFM–Hg model for 2017. The plots show the sum of all variables of MMHg in phytoplankton (Figs. S7–S10) and zooplankton (Figs. S11–S14).

The average MMHg_phy,P3 is lower in the Alb subbasin than in the Sad subbasin (Table 1), the latter being characterized by higher availability of MMHg at 35–45 m depth and higher biomasses of picophytoplankton during summer (Figs. S13 and S14). The hypothesis of the cell size effect on bioconcentration is the best synthesis of available knowledge (Schartup et al., 2018); however, the experimental study from which this relationship was derived also showed deviations from this pattern and temperature-induced effects, possibly due to phagocytosis of Hg-OM compounds by mixotroph organisms such as dinoflagellates (Lee and Fisher, 2016). One of the few field studies assessing the size effect on Hg bioaccumulation in plankton (Gosnell and Mason, 2015) found an increase in MeHg % but not in MeHg content for smaller plankton size (< 5 µm); however, the size classes analyzed by that study do not overlap well with the functional groups of the OGSTM–BFM–Hg model. More field, in vitro, and modeling studies are needed to disentangle the dynamics underlying key processes for Hg bioaccumulation.

Table 1 Indicators for MMHg bioaccumulation compared for the Alb and Sad subbasins (Fig. 3). The average (± SD) values of the plankton MMHg content (MMHg_PFT, pmol m⁻³), the MMHg quota (Q_MMHg, PFT, ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$), and the bioconcentration factors (log(BCF_PFT)) are shown for each phytoplankton and zooplankton PFT.

The highest concentrations of MMHg_zoo (Fig. 7b) are predicted from June to September at 35 m depth, showing a delay of about 1 month with respect to the temporal evolution of MMHg_phy due to the combined impacts of plankton phenology and trophic interactions. Heterotrophic nanoflagellates bioaccumulate much more MMHg (up to 0.5 pmol m⁻³, Fig. S19) than other zooplankton PFTs (up to 0.1 pmol m⁻³, Figs. S20–S22) because they feed on picophytoplankton (Fig. 8), which is more enriched in MeHg than other PFTs and is an important driver for the distribution of MMHg_zoo. However, in spite of heterotrophic nanoflagellates being always abundant from March to October (Fig. S23) and most abundant in the eastern subbasins, the highest MMHg bioaccumulation occurs when the increase in the nanoflagellates population is coincident with an enrichment of MMHg in phytoplankton (i.e., in western subbasins, as well as in the Adriatic and part of the Ionian Sea during summer).

Figure 8 Conceptual model of MMHg bioaccumulation in the plankton food web of the OGSTM–BFM–Hg model, subdivided into a lower “herbivorous” (red arrows) and upper “omnivorous” (blue arrows) food webs. The squares with bold numbers show estimated values of the bioconcentration factors (log(BCF)) for each PFT, calculated for the Sad subbasin from spatially averaged model output at monthly resolution. The green circles indicate the four phytoplankton PFTs (P1, P2, P3, and P4, representative respectively of diatoms, autotrophic nanoflagellates, picophytoplankton, and large phytoplankton), the light blue circles indicate the four zooplankton PFTs (Z3, Z4, Z5, and Z6, representative respectively of carnivorous mesozooplankton, omnivorous mesozooplankton, microzooplankton, and heterotrophic nanoflagellates), the purple circle indicates bacteria (B1, representative of all heterotrophic bacteria), and the grey circles indicate particulate organic carbon (POC, representative of all organic detritus) and semi-refractory dissolved organic matter (DOC_sr). The solid arrows show the preferential trophic interactions, the dashed arrows are for secondary trophic interactions (diet switch if the favorite preys are not available), and the circle dashed arrows are for cannibalization within the same zooplankton group. The dash and dotted arrows indicate passive uptake of MMHg by phytoplankton. The dotted arrows indicate the processes of transformation related to carbon recycling by bacteria (remPOC) and of Hg methylation, whose kinetics are coupled through the parameter k_met, converting Hg^II into MMHg.

The mean MMHg quota (Q_MMHg,PFT) in phytoplankton of the Mediterranean Sea is about 3.0 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for picophytoplankton, 0.3 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for autotrophic nanoflagellates, 0.03 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for diatoms, and 0.04 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$for large phytoplankton. The spatial variability is highlighted by comparing the quotas estimated for the Alb and Sad subbasins (Table 1) that are in the lower and upper range of the values for the Mediterranean Sea. Modeled values are comparable with the few observations available for open waters of the central Pacific Ocean (0.1–4 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$) (Gosnell and Mason, 2015) and of the northwest Atlantic Ocean (0.15 ± 0.06 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$) (Hammerschmidt et al., 2013), as well as with data from various coastal and shelf ecosystems (Harding et al., 2018). The speciation of Hg in phytoplankton groups of the Mediterranean Sea has never been assessed (Cinnirella et al., 2019). The MMHg quota was about 0.14 ± 0.1 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ in seston of diameter 80–200 µm sampled in the Gulf of Lion during spring and autumn 2004–2006 (Cossa et al., 2012).

The mean Q_MMHg,PFT values in zooplankton are 0.8 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for heterotrophic nanoflagellates, 0.3 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for microzooplankton, 0.6 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for omnivorous mesozooplankton, and 0.9 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for carnivorous mesozooplankton. The MMHg quota in carnivorous zooplankton (Q_MMHg,Z3) estimated for the Sad subbasin is much higher (3.60 ± 2.2) than the average value due to the high MMHg accumulation in the whole food web (Table 1). However, since biomasses of carnivorous zooplankton in the Sad are very low (Fig. S16), the overall content of MMHg in this group (MMHg_zoo,Z3) is lower (1.0 × 10⁻³ ± 2.7 × 10⁻⁴) than in the Alb subbasin (1.9 × 10⁻² ± 1.6 × 10⁻³).

The modeled zooplankton quotas compare well with the observations of MMHg in mesozooplankton of the Gulf of Lion, which was 0.52 ± 0.39 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for the size range 200–500 µm, 0.57 ± 0.39 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for the size range 500–1000 µm, and 1.8 ± 1.3 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ for zooplankton > 200 µm (Cossa, 2012). MMHg measured in the 1990s in zooplankton in the Gulf of Trieste (Nad, Fig. 1) Hg hotspot was 2.1 ± 2.3 ng g${}_{\mathrm{w}.\mathrm{w}.}^{-\mathrm{1}}$ (Horvat et al., 1999, 2014), consistently with the fact that the Gulf of Trieste is an upper bound for the Mediterranean, being subjected to high loadings of legacy Hg from the Soča/Isonzo river (Hines et al., 2000) and coastal lagoons (Melaku Canu et al., 2015; Rosati et al., 2020).

3.4 Bioconcentration factors and trophic magnification factors

Modeled log(BCFs) values for phytoplankton PFTs of the Mediterranean Sea range 3.7–6.5, which is slightly lower than the range (4.3–6.8) reported for phytoplankton from the central Pacific (Gosnell and Mason, 2015) and is higher than the ranges for Long Island Sound (2.6–5.5) (Gosnell et al., 2017) and other coastal and shelf sites of the world (3.1–4.4) (Harding et al., 2018). Higher BCFs in open waters than in coastal areas are consistent with the enhanced Hg bioavailability in the open sea inferred by other authors (Gosnell and Mason, 2015; Schartup et al., 2015), which in the model is approximated through the inverse relationship between MMHg uptake and DOC (Eq. 10). Modeled zooplankton log(BCFs) values in the Mediterranean Sea span from 4.0–6.4, comparable to the range for mesozooplankton (4.0–6.5) in the Pacific Ocean (Gosnell and Mason, 2015).

The trophic interactions among PFTs simulated by the model can be subdivided into a lower herbivorous part of the food web and an upper omnivorous part of the food web (Fig. 8).

In the lower part of the food web, heterotrophic nanoflagellates (2–20 µm) graze on bacteria (< 1 µm) and picophytoplankton (0.2–2 µm), while microzooplankton (20–200 µm) grazes on heterotrophic and autotrophic nanoflagellates (2–20 µm), diatoms (200–20 µm), and occasionally on other PFTs. In the upper part of the food web, the omnivorous mesozooplankton (200–2000 µm) preys on microzooplankton and all the phytoplankton PFTs > 2 µm, while carnivorous mesozooplankton (200–2000 µm) preys on omnivorous mesozooplankton and large phytoplankton (> 100 µm). The highest BCFs (and Q_MMHg,PFT) are predicted for picophytoplankton and carnivorous zooplankton (Table 1 and Fig. 8). Heterotrophic nanoflagellates display lower bioaccumulation than picophytoplankton because they also feed on bacteria and autotrophic nanoflagellates that dilute the MMHg intake (Fig. 8). Microzooplankton and omnivorous zooplankton, which have an intermediate position in the food web, have bioaccumulation levels comparable to those of heterotrophic nanoflagellates. The effect of biomagnification is visible only in carnivorous zooplankton, having much higher BCFs than their prey (Table 1). These dynamics are consistent with a global model taking into account a more complex plankton food web (366 PFTs) that found higher MMHg levels in phytoplankton than in their predators and biomagnification in the highest trophic levels with high food intake (Wu et al., 2021).

Figure 9 Trophic magnification factor (TMF, unitless) for the (a) omnivorous and (b) carnivorous plankton food webs (Fig. 8) for each subbasin of the Mediterranean Sea (Fig. 1), calculated by integrating over 0–100 m depth the monthly outputs of the OGSTM–BFM–Hg model.

Trophic magnification factors (TMFs) were calculated for each subbasin of the Mediterranean Sea for the lower and upper parts of the food web. Low TMFs (range 0.05–0.38) were estimated for the lower food web, indicating the absence of biomagnification (Fig. 9a), and higher TMFs (range 0.15–2.60) were estimated for the upper food web (Fig. 9b), in agreement with previous modeling studies suggesting that MMHg biomagnification between small zooplankton groups and phytoplankton is unlikely to happen (Wu et al., 2021, 2020; Zhang et al., 2020). Biomagnification (TMF > 1) in the upper food web (Fig. 9b) occurs only in subbasins with relatively abundant biomasses of carnivorous zooplankton (Fig. S24), namely in the Alb and Swm subbasins. The TMF for the Alb subbasin (2.4 ± 0.28) is lower than but comparable to the biomagnification from microseston to zooplankton (TMF = 4) estimated from field data for the northwest Atlantic Ocean (Hammerschmidt et al., 2013). The other Mediterranean subbasins, due to the absence or very low abundance of carnivorous mesozooplankton (< 0.05 mg_(C) m⁻³), have a shorter food web and lower TMFs (< 1). Although oligotrophy tends to favor longer food webs enhancing biomagnification, in extremely oligotrophic ecosystems the very low primary production can limit the presence of carnivorous zooplankton populations, lowering the bioaccumulation potential (Wu et al., 2021). The TMFs < 1 for the upper part of the food web are comparable to the estimates made with a global model for the most oligotrophic ocean regions such as the South Pacific Gyre, the North Atlantic Gyre, and the South Atlantic Gyre (Wu et al., 2021).

4 Conclusions

We developed and released the OGSTM–BFM–Hg model, a numerical state-of-the-art model tracking the biogeochemistry of the main marine Hg species (Hg^II, Hg⁰, MMHg, and DMHg) fully coupled to the biogeochemistry of plankton and nutrients. The model, applied to a high-resolution 3D domain of the Mediterranean Sea, successfully captured the zonal and vertical variability of MeHg concentrations, reproducing a decreasing MeHg gradient from west to east that is supported by the available observations. Sensitivity analysis showed that the agreement with observed MeHg concentrations and POC fluxes was improved by adopting an intermediate sinking speed for organic detritus. Modeled MeHg concentrations also improved by increasing the coefficient constant for the Hg methylation rate, while the inclusion of the dissolved organic carbon remineralization in the equation for Hg methylation increases the mismatch between the model and observations.

Sensitivity simulations on river Hg loadings corroborate the idea that in the Mediterranean Sea most of the inputs of particulate Hg (and POC) from rivers sink in coastal and shelf areas and have a limited impact on open-sea dynamics. However, as pointed out in a recent reassessment of global river inputs, there is a need to better assess the spatial and temporal variability of riverine loadings, especially in relation to extreme events and their increasing occurrence driven by climate changes. The net Hg exchange from Mediterranean waters to the atmosphere is lower than in a previous budget but in line with findings from a recent Hg isotopes study.

Model results highlight that summer stratification of the water column is an important process for the buildup of a subsurficial maximum of MeHg concentrations, while the occurrence of deep convection in winter results in substantial redistribution of MeHg smoothing the vertical profiles. A decrease in winter convection events linked to increasing water temperature, which has already been observed in recent years in the Mediterranean Sea, seems to limit MeHg redistribution in the water column, causing higher water MeHg concentrations in the biologically active zone and, in turn, higher plankton exposure. In fact, spatial and temporal variability of plankton bioaccumulation in the model is controlled by plankton phenology and by the availability of MMHg at the depths at which plankton blooms occur. The biomagnification potential is stronger in productive areas of the Mediterranean Sea characterized by high biodiversity and longer food web length with a significant presence of carnivorous zooplankton, such as the Alb and Swm subbasins.