Assessing the spatial and temporal variability of MeHg 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 modelling is 10 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 0 , 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 MeHg distribution and bioaccumulation. Model results reproduce the strong vertical and zonal gradients of MeHg concentrations 15 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 sub-surface MeHg maximum. The modelled bioaccumulation dynamics in plankton food webs are characterized by a high spatial and temporal variability that is driven by plankton phenology, and are in agreement with available field data of concentrations in plankton and with other indicators, such as bioconcentration factors (BCFs) and trophic magnification factors (TMFs). Model results pointed out that the 20 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.


Introduction
The anthropogenic alteration of 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., 2014aZhang et al., , 2014b) ) 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 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 (UN Environment 2019).
Here, we couple a model for Hg biogeochemistry (Canu et al., 2019;Rosati et al., 2020Rosati et al., , 2018;;Zhang et al., 2020) with the 3D transport biogeochemical model OGSTM-BFM (Lazzari et al., 2021(Lazzari et al., , 2010;;Salon et al., 2019) to investigate spatial and temporal variability of MeHg dynamics in the Mediterranean Sea.The coupled model (Fig. 1), which has a 1/16 o horizontal resolution, describes the biotic and abiotic interconversions among four Hg species (Hg II , MMHg, Hg 0 , and DMHg), as well as other key processes of Hg cycle such as partitioning to detritus, sinking to the seabed, exchange of gaseous Hg with the atmosphere, and bioaccumulation in 4 phyto-and 4 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.Model results are compared to available observations discussing the implications of spatial and temporal variability of modelled MeHg distributions in water and plankton.

Methods
The Biogeochemical Flux Model BFM (Vichi et al., 2015) is a model that describes 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 the same variable.In the BFM, PFTs have variable stoichiometry, and simulate the lower trophic levels of marine food webs up to carnivorous zooplankton (Sect.2.1).
The OGSTM-BFM model has been coupled off-line to the ocean general circulation model (NEMO) to investigate different issues related to the Mediterranean biogeochemistry (Canu et al., 2015;Cossarini et al., 2021Cossarini et al., , 2015;;Lazzari et al., 2021Lazzari et al., , 2016Lazzari et al., , 2014Lazzari et al., , 2012) ) 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,(Salon et al., 2019;Terzić et al., 2019).The integration of a Hg cycling module in the BFM-OGSTM model adds more flexibility to this system and open the way to novel operational oceanography applications.Here the OGSTM-BFM-Hg model is implemented (Fig. 1), by coupling the OGSTM-BFM model with the equations describing mercury (Hg) species partitioning, transformation, transport (Sect.2.2), bioconcentration in phytoplankton (Sect.2.3), and trophic transfer to zooplankton (Sect.2.4).The coupled model is used to investigate the Hg cycling and bioaccumulation dynamics (Sect.2.5) at the base of the Mediterranean Sea food webs.The model has a 1/16° horizontal resolution (approximately 6 km) and 72 vertical levels and runs with a 300 seconds time step.A model spin-up of 13 years forced with climatological data (Sect.2.6) is performed, and the results of last year of simulation (2017) are analysed.

BFM model dynamics
The BFM model dynamically simulates 9 plankton functional types (PFTs) representative of 4 phytoplankton groups (picophytoplankton, nanophytoflagellates, diatoms, and large phytoplankton), 4 zooplankton groups (eterotrophic nanoflagellates, microzooplankton, omnivorous mesozooplankton, and carnivorous mesozooplankton), and 1 group of heterotrophic bacteria (Vichi et al. 2015).The model reproduces the fluxes of carbon, nitrogen, phosphorous, and silicates from the inorganic nutrients 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).

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., 2020, 2018, in prep., Zhang et al., 2020) simulates the cycling of four mercury species in marine water (Fig. 2): inorganic divalent Hg (Hg II ), monomethylmercury (MMHg), elemental mercury (Hg 0 ), and dimethylmercury (DMHg).Eq. (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.(2) Desorption rates of Hg from particles in the ocean are negligible, and dynamics of Hg between the particulate and the dissolved phase can be synthetized as a balance between adsorption from the dissolved phase and particle remineralization (Lamborg et al., 2016).Accordingly, in the model the partitioning of Hg II and MMHg between particulate species (Hg II POC, MMHgPOC, fraction of each Hg species with respect to its total (fHg-POC, fHg-DOC, fHgCl, fMeHg-POC, fMeHg-DOC, fMegCl) is computed (e.g., Eq. 5-7 for Hg II ) from partition coefficients KD (Choe et al., 2003;Choe and Gill, 2003;Lamborg et al., 2016) that are set as model parameters (Tab.S1), and from DOC and POC concentrations.Then, concentrations of particulate and dissolved species of Hg II and MMHg are calculated from the corresponding fraction and the total concentration of the state variable (e.g., Eq. 8).
The pool of dissolved gaseous mercury species (DGM, i.e., Hg 0 and DMHg) is assumed to be entirely in the dissolved phase.are parameterized as functions of particulate organic matter remineralization (Zhang et al., 2020(Zhang et al., , 2014a)), 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) (Zhang et al., 2014a).MMHg photodegradation is assumed to yield both Hg 0 and Hg II (k phdm {MMHg DOC +MeHgCl}) due to the contrasting results in experimental studies (Luo et al., 2020).
The gas-exchange flux of Hg 0 with the atmosphere (k vol {Hg w 0 -Hg atm 0 /K H ' }) is using the Nightingale's parameterization (Nightingale et al., 2000) for the rate constant k vol , (Eq. ( 9)), which depends on wind speed (u10, m/s) and on the ratio of Hg Current knowledge on dynamics and concentrations of DMHg in the ocean and in the atmosphere does not allow a detailed modelling 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 uptaken from the water pool by

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 (eterotrophic 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 -3 d -1 ) and on the MMHg:C ratio of the preys that are grazed (

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 3 (w) were converted into ng gw.w -1 by computing the quota of MMHg (Q MMHg, PFT , Eq. ( 12)) to plankton wet weight (B w.w., PFT , gw.w./m 3 (w)) and correcting for the molar mass of Hg (200.59 g/mol).
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.Logarithmic scale is often used to handle BCF values (e.g.McGeer et al., 2003), and thus here the logarithm is included in Eq. ( 13).
BCF =log (10 6 Q MMHg, PFT MMHg w 200.59 ) The trophic magnification factor (TMF) is used to explore biomagnification across ecosystems, by expressing the increase of Hg along trophic levels (Alava et al., 2018;Harding et al., 2018).TMF is usually calculated in field studies of the wholefood web assessing both concentrations of pollutants and trophic levels though nitrogen isotopes.Here, in agreement with other modelling studies (Wu et al., 2021(Wu et al., , 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 by a lower trophic level "herbivorous" part and an upper "omnivorous" part.

Initial conditions, boundary conditions and forcings
Initial conditions, boundary conditions and physical forcings for the BFM implementation have monthly resolution, and are set in agreement with Lazzari et al., (2021).Initial conditions (Tab.S5) for Hg species in Mediterranean subbasins (Fig. 3) 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., 2015Kotnik et al., , 2007) ) and boundary conditions are set in agreement with the observations from the GEOTRACES GA03 meridional transect (Bowman et al., 2015).River load is estimated assuming a Hg concentration of 3 pM and 3.5% of MMHg (Cossa et al., 2017).A higher Hg concentration was chosen 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; https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.Rimondi et al., 2019), assuming MMHg was the 1% (Hines et al., 2000).The loading estimate is 6.57Mg/y, which would be about 70 kg lower (6.49Mg/y) if we used a uniform Hg concentration of 3 pM.The atmospheric Hg 0 concentrations is set to a constant value Hg atm 0 =1.6 ng/m 3 (Andersson et al., 2007;Fantozzi et al., 2013;Gårdfeldt et al., 2003), and the atmospheric deposition of HgT is set to 67.5 Mg/y, the 2% of which is MMHg (Cossa et al., 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 modelled fluxes (Rosati et al., 2018).

Model sensitivity to POC sinking velocity
The fit against experimental observation has been improved by exploring the importance of the sinking velocity parameter.
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.
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 current version of 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.The base model simulation was run by adopting the default POC sinking velocity of 3 m/d used in previous other works (Lazzari et al., 2021(Lazzari et al., , 2012) ) to model the biogeochemistry of the Mediterranean Sea.Two additional sensitivity simulations were run, by setting the POC sinking velocity to 10 and 20 m/d, assuming community-averaged cell diameters of 35 and 50 µm, respectively.

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

MeHg distribution in the water of 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).Modelled 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 North Western Mediterranean (NWM, Fig. 5b) and Ionian (ION, Fig. 5c) subbasins.
The sensitivity analysis showed that the sinking speed is an important control on the vertical distribution of MeHg (i.e., the shape of the vertical profile) in the water column but has little effect on the values of concentrations maxima (Supplemental Sect.S1 and Fig. S1).The sensitivity also showed that MeHg concentrations increase at depths where concentrations maxima have been observed (Cossa et al., 2009) if a threefold increase in the coefficient   is assumed (Supplemental Sect.S1 and Fig. S2), while the inclusion of DOC remineralization flux result in an excessive increase of surficial MeHg concentrations.
Modelled concentrations of HgT (Fig. S3) are within the experimental uncertainty of the data used to initialize the model, thus Hg methylation in the model does not appear to be limited by availability of inorganic Hg, and a calibration of the parameter   appears appropriate, considering that the initial guess for this parameter came 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).In spite of uncertainties remaining to be addressed, the OGSTM-BFM-Hg model is able to reproduce observed spatial gradients, testifying the improvements in capability to simulate Hg dynamics (Žagar et al., 2014, 2007) through the coupling with key biological processes, e.g. the carbon pump and microbial loop (Bowman et al., 2020;Cossa et al., 2017Cossa et al., , 2009;;Heimbürger et al., 2010;Monperrus et al., 2007;Munson et al., 2018;Sunderland et al., 2009).
Seasonality affects the modelled distribution of MeHg (Fig. 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 to 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 (Fig. 4, 5b and 5c) 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 biological active zone enhancing marine organisms exposure.
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 kmet (Fig. S4).However, while kmet are maxima in April, and net methylation fluxes (Fig. 6) are maxima in June (0.45pmol m -3 d -1 ), 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.Modelled 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. S4) 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. S4), due to the highest water temperature, irradiance, and deeper penetration of shortwave radiation due to the low plankton biomass.

Modelled MMHg in phytoplankton and zooplankton
The spatial and temporal distributions (Fig. 7) of modelled MMHg in phytoplankton (MMHgphy) and zooplankton (MMHgzoo) of the Mediterranean Sea are overall comparable to the concentrations of particulate MMHg (0.7 ± 1.3 pmol m -3 ) measured in the Atlantic Ocean (Bowman et al., 2015).The model reproduced a seasonal increase in MMHgphy (Fig. 7a) during late spring and summer that is driven by blooms of picophytoplankton (Fig. S5) combined with an increase in MMHg availability in surface water (Fig. S6).The amount of MMHg bioaccumulated by picophytoplankton (MMHgphy,P3, Fig. S7) is high (up to 0.6 pmol(hg) m(w) -3 ) compared to the bioaccumulation by other phytoplankton PFTs (up to 0.02 pmol(hg) m(w) -3 , Fig. S8-S10) due to the model assumption of increasing bioconcentration capacity with decreasing phytoplankton cell size (Sect.2.3).However, not all the picophytoplankton blooms results in high values of MMHgphy: the strongest blooms of picophytoplankton, occurring in March and April at 25 m depth (Fig. S5), 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. S6).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. S6).The average MMHgphy,P3 is lower in the ALB subbasin than in the SAD subbasin (Tab.1), the latter being characterized by higher availability of MMHg at 35-45 m depth and higher biomassess of picophytoplankton during summer (Fig. S5 and S6).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 dynoflagellates (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 analysed by that study do not overlap well with the functional groups of  (Gosnell & Mason, 2015), and of the Northwest Atlantic Ocean (0.15 ± 0.06 ng gw.w -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 gw.w -1 in seston of diameter 80-200 µm sampled in the Gulf of Lions during spring and autumn 2004-2006(Cossa et al., 2012)).
The modelled zooplankton quotas compares well with the observations of MMHg in mesozooplankton of the Gulf of Lions, which was 0.52 ± 0.39 ng gw.w -1 for the size range 200-500 µm, 0.57± 0.39 ng gw.w -1 for the size range 500-1000 µm, and 1.8 ± 1.3 ng gw.w -1 for zooplankton >200 µm (Cossa, 2012).MMHg measured in the '90s in zooplankton in the Hg hotspot of Gulf of Trieste (NAD, Fig. 3) was 2.1 ± 2.3 ng gw.w -1 (Horvat et al., 2014(Horvat et al., , 1999)), 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 Soca-Isonzo river (Hines et al., 2000) and coastal lagoons (Melaku Canu et al., 2015).

Bioconcentration Factors and Trophic Magnification Factors
Modelled BCFs 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 https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.
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) grazes 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 on 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 QMMHg,PFT ) are predicted for picophytoplankton and carnivorous zooplankton (Tab. 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 (Tab.1).These dynamics are consistent with a global model taking into account a more complex plankton food web (366 PFTs) that found highest MMHg levels in phytoplankton than in their predators, and biomagnification in the highest trophic levels with high food intake (Wu et al., 2021).
Trophic magnification factors (TMFs) were calculated for each subbasin of the Mediterreanean 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 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 modelling studies suggesting that MMHg biomagnification between small zooplankton groups and phytoplankton is unlikely to happen (Wu et al., 2021(Wu et al., , 2020;;Zhang et al., 2020).Biomagnification (TMF>1) in the upper food web (Fig. 9b) occurs only in Mediterranean subbasins with relatively abundant biomasses of carnivorous zooplankton (Fig. S16), 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 -3 ), have a shorter food web and lower TMFs (<1).Although oligotrophy tends to favour 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).Model results highlight that summer stratification of the water column is an important process for the buildup of a sub-surficial maximum of MeHg concentrations, while the occurrence of deep convection in winter result in substantial redistribution of MeHg smoothing the vertical profiles.A decrease of 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 waters 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.ALB 9.1 10 -2 (± 2.9 10 -2 ) 2.6 10 -2 (± 5.1 10 -3 ) 2.1 10 -3 (± 5.4 10 -3 ) 1.9 10 -2 (± 1.6 10 -3 ) SAD 7.2 10 -2 (± 7.9 10 -2 ) 1.6 10 -2 (± 1.7 10 -2 ) 7.1 10 -3 (± 7.1 10 -3 ) 1.0 10 -3 (± 2.7 10 -4 ) to organic detritus) and dissolved species (Hg II DOC, MMHgDOC, associated to dissolved organic carbon, and ionic species HgCln and MMHgCl) is calculated assuming repartition at thermodynamic equilibrium, taking into account dynamic POC remineralization and sinking based on the biogeochemical processes simulated in the BFM model (Fig.s 1 and 2).The https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.
Figure 6 also shows that net MeHg demethylation prevails throughout the year in deep waters (with the exception of a weak positive methylation in SAD and https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.TYR during summer), and for most of the year in surface waters.This suggests that the presence of MeHg in surface and deep water depends on its formation in the intermediate waters and diapycnal mixing.
OGSTM-BFM-Hg model.More field, in vitro and modelling studies are needed to disentangle the dynamics underlying key processes for Hg bioaccumulation.The highest concentrations of MMHgzoo (Fig.7b) are predicted from June to September at 35 m-depth, showing a delay of about one month with respect to the temporal evolution of MMHgphy due to the combined impacts of plankton phenology and https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.trophic interactions.Heterotrophic nanoflagellates bioaccumulate much more MMHg (up to 0.5 pmol/m3, Fig.S11) than other zooplankton PFTs (up to 0.1 pmol m-3, Fig.S12-S14) because they feed on picophytoplankton (Fig.8) which is more enriched in MeHg than other PFTs, and are an important driver for the distribution of MMHgzoo.However, in spite of heterotrophic nanoflagellates being always abundant from March to October (Fig.S15), and most abundant in the eastern subbasins, the highest MMHg bioaccumulation occurs when the increase of 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).The mean MMHg quota (QMMHg,PFT ) in phytoplankton of the Mediterranean Sea is about 3.0 ng gw.w -1 .for picophytoplankton, 0.3 ng gw.w -1 for autotrophic nanoflagellates, 0.03 ng gw.w -1 for diatoms, and 0.04 ng gw.w -1 for large phytoplankton.The spatial variability is highlighted by the comparison of the quotas estimated for the ALB and SAD subbasins (Tab. 1) that are in the lower and upper range of the vales for the Mediterranean Sea.Modelled values are comparable with the few observations available for open waters of the central Pacific Ocean (0.1-4 ng gw.w -1 ) https://doi.org/10.5194/bg-2022-14Preprint.Discussion started: 14 February 2022 c Author(s) 2022.CC BY 4.0 License.4ConclusionsWe developed the OGSTM-BFM-Hg model, a numerical state of the art model tracking the biogeochemistry of the main marine Hg species (Hg II , Hg 0 , MMHg, and DMHg) 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 observations of MeHg was improved by adopting an intermediate sinking speed for organic detritus and by increasing the rate constant for Hg methylation, while the inclusion in the equation for methylation of the remineralization of dissolved carbon (DOC) increases the mismatch between model and observations, suggesting that this process might not be important in shaping MeHg distributions.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Hg dynamics in the coupled model OGSTM-BFM-Hg.The Hg model simulates the cycling of inorganic (Hg II and Hg 0 ) and 630

Figure 5 :
Figure 5: Vertical profiles of monthly averaged MeHg (MMHg+DMHg) concentrations in the water column reproduced by 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. 3): (a) Southern Adriatic Sea (SAD), (b) North Western Mediterranean (NWM), and (c) Ionian Sea (ION2).The right panels show the entire water column and the left panels show a detail of modelled concentrations in the top 600 m of the same subbasins.

Figure 6 :
Figure 6: Monthly evolution of the net methylation fluxes (pmol m -3 d -1 ) 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 modelled Hg methylation fluxes and demethylation 655

Figure 8 .
Figure 8. Conceptual model of MMHg bioaccumulation in the plankton food web of the SAD subbasin of the OGSTM-BFM-Hg model, subdivided into a lower "herbivorous" (red arrows) and upper "omnivorous" (blue arrows) food webs.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 kmet, converting Hg II into MMHg.The squares with bold numbers show estimated values of BCF (Bioconcentration factors, log(kg l -1 )) for each PFT, calculated from spatially averaged model output at monthly resolution.The green circles indicate the four phytoplankton PFTs (P1, P2, P3, P4, representative respectively of diatoms, autotrophic nanoflagellates, picophytoplankton and large phytoplankton), the light blue circles indicate the four zooplankton PFTs (Z3, Z4, Z5, Z6, representative respectively of carnivorous mesozooplankton, omnivorous mesozooplankton, microzooplankton, and heterotrophic nanoflagellates), the purple circle indicate 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 (DOCsr).

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