The 3D biogeochemical marine mercury cycling model MERCY v2.0 – linking atmospheric Hg to methylmercury in fish

Mercury (Hg) is a pollutant of global concern. Due to anthropogenic emissions, the atmospheric and surface ocean Hg burden has increased substantially since preindustrial times. Hg emitted into the atmosphere gets transported on a global scale and ultimately reaches the oceans. There it is transformed into highly toxic methylmercury (MeHg) that effectively accumulates in the food web. The international community has recognized this serious threat to human health and in 2017 regulated Hg use and emissions under the UN Minamata Convention on Mercury. Currently, the first effectiveness evaluation of the Minamata Convention is being prepared, and, in addition to observations, models play a major role in understanding environmental Hg pathways and in predicting the impact of policy decisions and external drivers (e.g., climate, emission, and land-use change) on Hg pollution. Yet, the available model capabilities are mainly limited to atmospheric models covering the Hg cycle from emission to deposition. With the presented model MERCY v2.0 we want to contribute to the currently ongoing effort to improve our understanding of Hg and MeHg transport, transformation, and bioaccumulation in the marine environment with the ultimate goal of linking anthropogenic Hg releases to MeHg in seafood.

Here, we present the equations and parameters implemented in the MERCY model and evaluate the model performance for two European shelf seas, the North and Baltic seas. With the model evaluation, we want to establish a set of general quality criteria that can be used for evaluation of marine Hg models. The evaluation is based on statistical criteria developed for the performance evaluation of atmospheric chemistry transport models. We show that the MERCY model can reproduce observed average concentrations of individual Hg species in water (normalized mean bias: Hg_T 17 %, Hg⁰ 2 %, MeHg −28 %) in the two regions mentioned above. Moreover, it is able to reproduce the observed seasonality and spatial patterns. We find that the model error for Hg_T(aq) is mainly driven by the limitations of the physical model setup in the coastal zone and the availability of data on Hg loads in major rivers. In addition, the model error in calculating vertical mixing and stratification contributes to the total Hg_T model error. For the vertical transport we find that the widely used particle partitioning coefficient for organic matter of log(k_d)=5.4 is too low for the coastal systems. For Hg⁰ the model performance is at a level where further model improvements will be difficult to achieve. For MeHg, our understanding of the processes controlling methylation and demethylation is still quite limited. While the model can reproduce average MeHg concentrations, this lack of understanding hampers our ability to reproduce the observed value range. Finally, we evaluate Hg and MeHg concentrations in biota and show that modeled values are within the range of observed levels of accumulation in phytoplankton, zooplankton, and fish. The model performance demonstrates the feasibility of developing marine Hg models with similar predictive capability to established atmospheric chemistry transport models. Our findings also highlight important knowledge gaps in the dynamics controlling methylation and bioaccumulation that, if closed, could lead to important improvements of the model performance.

1 Background

Mercury (Hg) is a global pollutant and a dangerous neurotoxin (AMAP/EMEP, 2019a). Since preindustrial times, the global Hg cycle has been significantly altered by anthropogenic emissions (Streets et al., 2019) resulting in a 3-fold pre-anthropogenic-to-present-day increase in the atmospheric and a substantial increase in oceanic Hg burden (Lehnherr, 2014; Amos et al., 2013). The major anthropogenic sources of Hg are emissions from coal-fired power plants, small-scale artisanal gold mining, and metal and cement production (Pirrone et al., 2010; AMAP/EMEP, 2013, 2019a, b). In addition, natural emissions and legacy re-emissions from previously deposited Hg (most of it of anthropogenic origin) also contribute significantly to the atmospheric Hg burden (Pirrone et al., 2010; Driscoll et al., 2013; Obrist, 2018). The atmospheric lifetime of Hg is estimated to be in the range of 0.6 to 1.0 years (Slemr et al., 2018), resulting in a global atmospheric distribution of Hg. Atmospheric Hg will eventually be deposited (Cohen et al., 2016; Jiskra et al., 2018). A large fraction is deposited directly into the ocean, but Hg deposited onto land can also be transported to the ocean via rivers and groundwater. In the aqueous phase, inorganic Hg can be methylated, forming the highly bioaccumulative monomethylmercury (MMHg) and/or dimethylmercury (DMHg). These MeHg compounds are readily accumulated in the food web and pose a risk to food safety and human health (Clarkson, 1990; Mason et al., 1996; Chen et al., 2012; Parks et al., 2013; Puty et al., 2019). Because of this, the international community, under the umbrella of the United Nations Environmental Programme (UNEP), signed the Minamata Convention on Mercury, which came into force in 2017. Under this convention, all participating 184 nations have agreed to assess Hg pollution under their jurisdiction, to minimize Hg usage and release of Hg compounds into the environment, and to regularly assess the impact of the reduction measures taken on the environmental Hg burden and distribution. In order to assess the impact of reduction measures, there is an urgent need to understand the Hg pathways from anthropogenic releases to top predators and humans, with specific attention to the marine ecosystem.

In this paper, we (1) introduce a newly developed numerical multi-compartment model for Hg cycling in the marine environment including accumulation in the marine food web (MERCY v2.0) and (2) evaluate the model performance to reproduce observed concentrations of, seasonality of, and variability in Hg species. For the latter, we apply performance criteria used for evaluation of atmospheric chemistry transport models and also for evaluation of marine Hg models. We use these criteria to (2.1) quantify the models' predictive capabilities based on our current understanding of Hg cycling, (2.2) identify the major sources of model error, and (2.3) quantify the constraints on model improvement based on current process understanding and measurement availability and uncertainty. With this study, we present an evaluation of our marine Hg model and a general framework that provide the basis for future intercomparison studies of marine Hg models.

1.1 Research question

The key question concerning Hg pollution is how changing Hg emissions and other external stressors such as climate and land-use change impact MeHg accumulation in seafood, which is an important global protein source for human consumption (Pauly et al., 2002; Obrist et al., 2018). To anticipate the natural Hg cycle and to identify the impact of human actions on the system, it is necessary to develop multi-compartment chemistry transport models (CTMs) including all relevant compartments: atmosphere, soil/vegetation, rivers and oceans, sediments, and the marine ecosystem. The need to incorporate all compartments into a single multi-compartment model arises from the fact that Hg is non-degradable and constantly cycling between environmental compartments, unlike most pollutants, which tend to accumulate in a single compartment and/or degrade over time. For example, atmospheric deposition of oxidized Hg is a major flux of Hg into the ocean, but reduction reactions in the ocean and the high vapor pressure of elemental Hg⁰ also result in a constant release of Hg from ocean to atmosphere (Fitzgerald et al., 1984; Fitzgerald and Kim, 1986; Andersson et al., 2008; Mason et al., 2012). The only real sink for Hg in the environment is burial in the lithosphere, mainly as stable cinnabar (HgS) in anoxic marine sediments. Thus, coupled Earth system models are needed to gain a deeper understanding of the processes and dynamics governing transport of Hg, Hg methylation, and the variability in Hg accumulation in the marine food web. While there are a large number of emission inventories and atmospheric CTMs, there are still only a limited number of CTMs with a focus on marine Hg cycling and food web transfers.

1.2 Development and state of the art in Hg modeling

Atmospheric Hg modeling is well established, and a large variety of global (ECHMERIT: Jung et al., 2009; De Simone et al., 2014; GLEMOS: Travnikov and Ilyin, 2009; Travnikov et al., 2009; GEM-MACH: Durnford et al., 2012; Kos et al., 2013; Dastoor et al., 2015; GEOS-Chem: Holmes et al., 2010; Amos et al., 2012; Song et al., 2015) and regional (CMAQ: Bullock et al., 2008; Bash, 2010; Zhu et al., 2015; DEHM: Christensen et al., 2004; WRF-Chem: Gencarelli et al., 2017) atmospheric CTMs for Hg cycling have been published. Due to this abundance, many model intercomparison and source apportionment studies have improved our understanding of atmospheric Hg transport and source–receptor relationships and have allowed us to predict future atmospheric Hg levels and deposition fluxes (Bergan et al., 1999; Xu et al., 2000; Petersen et al., 2001; Lee et al., 2001; Seigneur et al., 2001; Bullock and Brehme, 2002; Dastoor et al., 2002; Hedgeock et al., 2004; Selin et al., 2007; Travnikov et al., 2009; Bieser et al., 2014; Gencarelli et al., 2017; Dastoor et al., 2015; Song et al., 2015; Cohen et al., 2016; Travnikov et al., 2017; Bieser et al., 2017; Horowitz et al., 2017). These models and studies are a keystone in informing policymakers to support the implementation and effectiveness evaluation of the Minamata Convention (http://www.mercuryconvention.org, last access: 1 April 2023).

Compared to Hg modeling in the atmosphere, marine Hg modeling is still in its infancy, and only a limited number of models exist so far. The development of marine Hg models can be divided into four phases. At first, the ocean was implemented as a boundary for atmospheric CTMs, and nowadays most atmospheric CTMs implement some kind of surface ocean parameterization to explicitly include Hg air–sea exchange. One of the earliest marine Hg model developments were box models (Sunderland and Mason, 2007), followed by the addition of inorganic Hg redox chemistry and transport in a 2D slab ocean model coupled to the GEOS-Chem model (Selin et al., 2008; Strode et al., 2007; Soerensen et al., 2010). The aim of these early models was to improve air–sea exchange by including horizontal transport, redox chemistry, and river loads. Next came the development of the first marine 3D models. These models, still limited to the inorganic Hg cycle, were used to investigate marine Hg dynamics (Zhang et al., 2014a, b; Bieser and Schrum, 2016). In the next stage, several specialized marine Hg models were developed which were not based on 3D hydrodynamic models. Soerensen et al. (2016a) published a coupled physical–biogeochemical multi-box model including organic Hg chemistry to investigate the Hg budgets in the Baltic Sea. Focusing on bioaccumulation, Schartup et al. (2018) implemented Hg accumulation in a complex food web model and Sunderland et al. (2009, 2018) modeled the consumer exposure to MeHg in seafood. Finally, Pakhomova et al. (2018) developed a model with comprehensive Hg chemistry based on a hydrodynamic 1D model. Only in recent years has the development of comprehensive marine Hg models gained traction. So far, four marine Hg models based on numerical hydrodynamic 3D models have been published (Semeniuk and Dastoor, 2017; Zhang et al., 2019; Kawai et al., 2020; Rosati et al., 2022). All of these models include a complete marine Hg chemistry including MeHg. Yet, only Zhang et al. (2020) and Rosati et al. (2022) have also implemented Hg cycling into a biogeochemical model considering uptake to and release from marine biota, making this model the first hydrodynamic 3D Hg model to include the marine ecosystem.

1.3 Our contribution to the presented problem

Here we present our newly developed biogeochemical multi-compartment model for Hg cycling, MERCY v2.0, and evaluate its predictive capabilities and limitations using evaluation criteria applied for performance evaluation of atmospheric CTMs (Derwent et al., 2010; Thunis et al., 2012, 2013; Carnevale et al., 2014). We focus on the implementation of the marine Hg cycle including a comprehensive marine Hg chemistry and partitioning scheme as well as bioconcentration and biomagnification. We improve on the state of the art by introducing an experimental upper trophic layer that simulates Hg and MeHg accumulation in fish. To our knowledge, MERCY v2.0 includes all currently known processes controlling marine Hg cycling. The model is based purely on processes, reactions, and rates published in peer-reviewed literature, and no additional model tuning was performed.

We investigate the model predictive capabilities, something we consider important before using the model to study budgets or global dynamics. This allows us to quantify our model uncertainty, which for other models has only been loosely constrained to be “orders of magnitude” (Kawai et al., 2020), and discuss the processes and parameters driving it. Set up on a high-resolution regional domain covering a wide range of marine regimes in a region with high primary productivity and a relative abundance of observations, we evaluate the ability of the model to reproduce observed concentrations of, seasonality of, and variability in individual marine Hg species. Using common practice from atmospheric Hg modeling, we establish a quantitative benchmark for the capability of the model to reproduce actual observations of marine Hg concentration and speciation. Based on this we discuss the major knowledge gaps and research questions that need to be tackled in order to improve our understanding of marine Hg cycling. Our ultimate goal is to improve capabilities to link changes in external stressors like anthropogenic emissions and climate change to MeHg accumulation in the marine food web by providing an independent model for marine Hg cycling and by fostering collaboration in the form of model intercomparison studies comparable to the efforts in atmospheric Hg modeling (Ryaboshapko et al., 2002; Bullock et al., 2008; Travnikov et al., 2017; Bieser et al., 2017). Finally, we want to identify and communicate the major needs for monitoring of Hg species in the marine environment.

2 Model description

2.1 Model framework

The marine Hg chemistry scheme we develop for MERCY v2.0 is embedded into GCOAST (Geesthacht Coupled cOAstal model SysTem), a modeling framework coupling physical, chemical, and biological numerical models. It is an update and overhaul of MERCY v1.0 (Bieser and Schrum, 2016), which featured only inorganic Hg chemistry and no ecosystem interactions. As input, MERCY uses hourly model output from four types of 3D hydrodynamic model (atmospheric physics, atmospheric chemistry, marine physics, and marine ecosystem) to drive the marine Hg speciation, transport, and bioaccumulation model. While this approach requires a large amount of storage capacity, it reduces the computational requirements and allows the model to be easily run with input from alternative biogeophysical models. The external variables used by MERCY are listed in Table 1. In brief, the models used in this work are as follows:

1. The regional weather and climate model COSMO-CLM (Rockel et al., 2008; Sørland et al., 2021) provides meteorological variables used to calculate air–sea exchange (temperature and wind speed) and photolytic reactions (surface shortwave radiation). COSMO-CLM is nudged to the atmospheric reanalysis dataset ERA-Interim (Berrisford et al., 2011; Dee et al., 2011; Hersbach et al., 2020).

2. The atmospheric chemistry transport model CMAQ-Hg (Byun and Schere, 2006; Zhu et al., 2015; Bieser et al., 2016) is forced by COSMO-CLM meteorology and used to calculate atmospheric transport, chemistry, particle partitioning, and deposition for atmospheric trace gases. MERCY uses atmospheric Hg concentrations and deposition fluxes from CMAQ-Hg.

3. The physical hydrodynamic ice–ocean model HAMSOM (Backhaus, 1983; Schrum and Backhaus, 1999) is directly coupled to the ecosystem model ECOSMO, enabling it to represent the impact of the ecosystem on the hydrodynamics (e.g., light attenuation by biota). In MERCY the physical variables are used to calculate marine mercury transport as well as temperature and salinity dependence of mercury cycling and speciation. The HAMSOM advection scheme is used to transport all Hg state variables. The model setup is based on GEBCO bathymetry data (GEBCO Bathymetric Compilation Group).

4. The marine end-to-end NPZD (nutrient, phytoplankton, zooplankton, detritus) ecosystem model ECOSMO (Schrum et al., 2006; Daewel and Schrum, 2013; Daewel et al., 2019) is a 3D-resolved food web model directly coupled with HAMSOM. It includes nutrients (nitrogen, phosphorus, and silica) and a food web based on a functional group approach with three phytoplankton species (diatoms, flagellates, and cyanobacteria), two zooplankton species (herbivore and omnivore), a macrobenthos group, and a pelagic fish group representing higher trophic levels. Additionally, oxygen, biogenic opal, detritus, and dissolved organic matter are considered, and the model includes a two-layer sediment compartment to simulate sedimentation and resuspension. In MERCY detritus and dissolved organic matter determine the partitioning of Hg and MeHg, and factors such as light attenuation and oxygen concentration influence Hg speciation. Moreover, concentrations of the various species of the model food web are used to calculate bioconcentration and biomagnification of Hg and MeHg.

All employed models and data are freely available (see “Code availability” and “Data availability” sections at the end of the paper).

Table 1 MERCY input variables and source models.

^∗ Rates for macrobenthos are in mg C m⁻² s⁻¹.

2.2 General equations

MERCY v2.0 implements all processes we identified as relevant to marine (pelagic and benthic) Hg cycling into a 3D ocean-ecosystem model. MERCY is based on basic principles describing Hg transport, transformation, and bioaccumulation. It is set up on the same grid and domain as the coupled ocean-ecosystem model HAMSOM-ECOSMO. Based on archived hourly HAMSOM-ECOSMO output, it is effectively offline-coupled to the marine hydrodynamic and ecosystem models. The HAMSOM-ECOSMO model has been shown to accurately reproduce ecosystem dynamics in the coupled North Sea–Baltic Sea system. The model equations and a model validation on the basis of nutrients are presented in detail by Daewel and Schrum (2013), who showed that the model can reasonably simulate ecosystem productivity in the North Sea and the Baltic Sea on seasonal to decadal timescales. Using the same numerical approximations as described in Daewel et al. (2019), the rate of change in the concentration of Hg state variables over time $\frac{\mathit{\delta }C}{\mathit{\delta }t}$ is estimated by the prognostic equation (Eq. 1). This rate of change is subsequently integrated over the internal time step and applied to the corresponding state variables.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathit{\delta }C}{\mathit{\delta }t}}=V\mathrm{\nabla }C+w_{\mathrm{d}}{\displaystyle \frac{\mathit{\delta }C}{\mathit{\delta }z}}+{\displaystyle \frac{\mathrm{d}z}{\mathit{\delta }z}}\left(A_v{\displaystyle \frac{\mathit{\delta }C}{\mathit{\delta }z}}\right)+R_{(C,B)}\end{array}$$

The physical transport terms for advection V∇C with 3D velocity field $V=(u,v,w)$, vertical transport $w_{\mathrm{d}}\frac{\mathit{\delta }C}{\mathit{\delta }z}$ with sinking velocity w_d, and turbulent mixing $\frac{\mathrm{d}z}{\mathit{\delta }z}\left(A_v\frac{\mathit{\delta }C}{\mathit{\delta }z}\right)$ with diffusion coefficient A_v and velocity V are calculated by the hydrodynamic host model. At the upper and lower boundary of the water column, boundary conditions are presented to account for air–sea exchange (Sect. 2.3.6) and sedimentation and resuspension (Sect. 2.3.5). Each Hg state variable C is subject to additional transformations R_(C,B), which include chemical transformations Rc_(C) (Sect. 2.3.1), partitioning Rp_(C) (Sect. 2.3.2), and biological uptake Rb_(C,B) by ecosystem group B (Sect. 2.3.4) (Eq. 2). Marine biota are implemented in the ecosystem model following a functional group approach further described in Sect. 2.3.4. All transformations R_(C,B) are mass-conserving transfer reactions, which means that besides emission inputs and inflow/outflow at the domain boundaries, no Hg is added or removed from the system. The exact formulation of R_(C,B) differs for each Hg species in the model. In this section, we give a general overview of all possible transformations, and the exact formulae and parameterizations are given in Sect. 2.3. A complete list of all Hg state variables is given in Table 2. All chemical reactions Rc_(C) and their respective reaction rates can be found in Table 3, and further physical and biological parameters are given in Table 4.

$$\begin{array}{}\text{(2)}& R_{\left(C,B\right)}={\mathrm{Rc}}_{\left(C\right)}+{\mathrm{Rp}}_{\left(C\right)}+{\mathrm{Rb}}_{(C,B)}\end{array}$$

Chemical transformations Rc_(C) (Eq. 3) are the sum of all reactions, where species C is a reaction product $\sum _{i=\mathrm{0}}^n{k}_i{C}_i$ of another species C_i with reaction rate k_i minus the sum of all reactions where C is an educt $\sum _{j=\mathrm{0}}^n{k}_{j}C$ with reaction rate k_j. Chemical reactions are implemented as pseudo-first-order reactions $\frac{\mathit{\delta }C}{\mathit{\delta }t}=kC$ using either a fixed reaction rate k₁ or a dynamic reaction rate k₂=k₁C₂ dependent on a second reactant C₂ or an associated environmental variable (e.g., temperature). For photolytic reactions the reaction rate is k=k_pE_λ with the integrated photon flux $E_{\mathit{\lambda }}=\underset{{\mathit{\lambda }}_{\mathrm{0}}}{\int }{}_n^{\mathit{\lambda }}{E}_{\mathit{\lambda }}$ for specific wavelengths λ and the photolysis rate k_p.

$$\begin{array}{}\text{(3)}& {\mathrm{Rc}}_{\left(C\right)}=\sum _{i=\mathrm{0}}^n{k}_i{C}_i-\sum _{j=\mathrm{0}}^n{k}_{j}C,\end{array}$$

where n is the number of Hg species.

Partitioning Rp_(C) (Eq. 4) describes sorption and desorption of dissolved C_aq to particulate organic matter (POM) and dissolved organic matter (DOM), where C_POM is particulate Hg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+}$ and C_DOM is Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ bound to DOM. The equilibrium between these species is described by sorption and desorption rates k_s and k_d.

$$\begin{array}{}\text{(4a)}& {\displaystyle }\begin{array}{rl}{\mathrm{Rp}}_{\left(c\right)}\left(C_{\mathrm{aq}}\right)& =k_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{POM}}+k_{\mathrm{d}\mathrm{2}}{C}_{\mathrm{DOM}}\\ & -k_{\mathrm{s}\mathrm{1}}{C}_{\mathrm{aq}}\text{POM}-k_{\mathrm{s}\mathrm{2}}{C}_{\mathrm{aq}}\text{DOM}\end{array}\text{(4b)}& {\displaystyle }{\mathrm{Rp}}_{\left(c\right)}\left(C_{\mathrm{POM}}\right)=k_{\mathrm{s}\mathrm{1}}{C}_{\mathrm{aq}}\text{POM}-k_{\mathrm{d}\mathrm{1}}{C}_{\mathrm{POM}}\text{(4c)}& {\displaystyle }{\mathrm{Rp}}_{\left(c\right)}\left(C_{\mathrm{DOM}}\right)=k_{\mathrm{s}\mathrm{2}}{C}_{\mathrm{aq}}\text{DOM}-k_{\mathrm{d}\mathrm{2}}{C}_{\mathrm{DOM}}\end{array}$$

Biological uptake of Hg into biota Rb_(C,B) (Eq. 5) includes two distinct processes: (1) bioconcentration, which is defined as the passive uptake of dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ through the cell membrane of a functional ecosystem group B, and (2) biomagnification, which is the sum of active uptake and release through feeding. For higher trophic levels, the Hg in biota from active and passive uptake is stored in separate state variables with different release rates due to the differing accumulation patterns for each uptake process.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{Rb}}_{(C,B)}& =\sum _{i=\mathrm{0}}^n\left(v_i{A}_B{C}_i\right)+\sum _{b=\mathrm{0}}^m(r_{B,b}{\mathit{ϵ}}_C{C}_{\mathrm{b}}-r_{\mathrm{b},B}{C}_B)\\ & -(r_{\mathrm{r}\left(B\right)}+r_{m\left(B\right)})C_B,\end{array}\end{array}$$

where n is the number of Hg species and m is the number of ecosystem groups.

Bioconcentration $\sum _{i=\mathrm{0}}^n\left(v_i{A}_{\mathrm{b}}{C}_i\right)-r_{\mathrm{r}}{C}_{\mathrm{b}}$ is the sum of passive uptake with an uptake rate r_u=v_iA_b depending on the permeation velocity v_i of dissolved Hg species C_i and the average ecosystem group surface area A_b minus an ecosystem group and Hg-species-dependent release rate r_r multiplied with the Hg concentration inside biota C_b.

Biomagnification $\sum _{b=\mathrm{0}}^m{r}_{B,b}{\mathit{ϵ}}_C{C}_{\mathrm{b}}-r_{b,\mathrm{B}}{C}_B$ describes the active transfer of Hg driven by feeding rates r_B,b of an ecosystem group B on other ecosystem groups b and the corresponding feeding pressure r_b,B. The efficiency of Hg transfer upon feeding is determined by a Hg-species-dependent uptake efficiency ϵ_C.

Additional release from the biological matrix C_b is described by a mortality rate r_m. For the release of Hg from detritus into the dissolved Hg pool r_m, we have a temperature-dependent remineralization rate k_rem (see Eq. 9 in Sect. 2.3.1).

Finally, the respective change in dissolved Hg concentrations C_aq due to uptake into and release from marine biota is given by Eq. (6), where $\sum _{b=\mathrm{0}}^m\mathit{\left\{}{r}_{B,b}\right(\mathrm{1}-{\mathit{ϵ}}_C\left)C_{\mathrm{b}}\mathit{\right\}}$ is the Hg fraction directly excreted into the dissolved phase upon feeding of ecosystem group B on another ecosystem group b and (r_r(B)+r_m(B))C_B is the release due to a constant release rate r_r(B) and the mortality rate r_m(B) of Hg species C_b in ecosystem group B.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{Rb}}_{(C,B)}& =\sum _{b=\mathrm{0}}^m\mathit{\left\{}{r}_{B,b}\right(\mathrm{1}-{\mathit{ϵ}}_C\left)C_{\mathrm{b}}\mathit{\right\}}\\ & -\sum _{i=\mathrm{0}}^n\mathit{\left\{}{v}_i{A}_B{C}_i\mathit{\right\}}+(r_{\mathrm{r}\left(B\right)}+r_{\mathrm{m}\left(B\right)})C_B,\end{array}\end{array}$$

where n is the number of Hg species and m is the number of ecosystem groups.

2.3 Implemented processes

MERCY implements Hg using 35 variables (Table 2) representing different Hg species in the atmosphere, ocean, and sediment. For each model time step and each grid cell, the species are redistributed accounting for mass conservation based on physical, chemical, and biological processes. Figure 1 gives a graphical overview of transformations between Hg species in MERCY.

Figure 1 Schematic of the chemical mechanism in MERCY. Solid lines indicate chemical reactions, fine dotted lines photolytic reactions, dash-dotted lines instantaneous partitioning processes, and bold dotted lines bioaccumulation and releases from biota into the dissolved phase. Colors codes are white for elemental mercury, yellow for inorganic oxidized mercury, pink for methylated mercury, and green for Hg in biota. The physical state of each species is indicated by “g” for gaseous, “aq” for dissolved, and “s” for solid. The upper row indicates Hg species in the atmosphere, and the lower row indicates those in the sediment. All species and their reactions are given in Tables 2 and 3. Note Reaction (R20) (reductive methylation, Table 3) MMHg–DOM → Hg⁰ extends to the left edge of the figure.

Table 2 Hg species in MERCY. Species can represent state variables in multiple models.

^∗ Hg species in biota (Hg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+}$ and MMHg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+})$ represent one state variable for each functional group in the ecosystem model ECOSMO (see Sect. 2.3.4), giving a total of 40 mercury species.

2.3.1 Chemistry

In this section, we present all chemical state variables and the transformation processes in the model. A complete overview of all chemical transformations and the respective reaction rates k is given in Table 3. All chemical transformations are calculated using pseudo-first-order reactions following Eq. (7). The chemical mechanism is implemented using a tendency approach, where the relative change for each reaction is calculated and all changes to state variables are applied simultaneously. Equation (3) gives the change to the concentration of a single Hg species due to all reactions depleting and producing it. We run the chemistry module with a time step of 60 s but find that it runs stably and efficiently even with much larger time steps of 600 s.

$$\begin{array}{}\text{(7)}& C_t=C_{\mathrm{0}}{e}^{-tk},\end{array}$$

where C_t is the concentration at time t [ng L⁻¹], C₀ is the concentration at time 0 [ng L⁻¹], t is time [s], and k is a pseudo-first-order reaction constant [s⁻¹].

Table 3 Chemical reactions as implemented in the MERCY model.

Pseudo-first-order reaction rates k₂=k₁C depend on the following variables C: ^a temperature-dependent reaction rate, ^b cyanobacteria-concentration-dependent reaction rate, ^c reaction rate dependent on photolytically active radiation, ^d sulfate-concentration-dependent reaction rate, and ^e oxygen-dependent concentration reaction rate. ^f Reaction starts with a lower rate under hypoxic conditions.

Redox reactions

Hg redox chemistry is implemented with five reactions. Reduction (Hg²⁺ → Hg⁰) is driven by three processes: (R1) a continuously ongoing chemical reduction, often referred to as dark reduction; (R2) photolytic reduction; and (R3) biogenic reduction (Table 3). We use reaction rates reported by Kuss et al. (2015). This leads to each reduction reaction being of roughly similar importance for the total Hg⁰ production, albeit with specific distinct seasonality (note that the biogenic reduction only plays a role in the Baltic Sea due to cyanobacteria). This is in contrast to other published reaction rates where photolysis is the dominant pathway (Qureshi et al., 2009). We do not use an intermediate oxidation product (Hg^∗) as we found the species to be too short-lived for the given model setup. We chose the values from Kuss et al. (2015) as, contrary to other studies, these were measured under in situ conditions. The oxidation is driven mainly by chemical oxidation (R4), while photolytic oxidation (R5) rates are much smaller, leading to a net photolytic reduction. The photolysis rates are parameterized to the photolytically active radiation based on observations. The biogenic reduction reaction rate is scaled by the cyanobacteria biomass and is not triggered by other phytoplankton species (Kuss et al., 2015). For the chemical reduction, we consider a temperature-dependent reaction rate k_rd defined as 100 % at 15 ^∘C (50 % at 0 ^∘C and 125 % at 20 ^∘C) (Kuss et al., 2015) (Eq. 8). Finally, we consider reductive demethylation of MeHg⁺ (Reaction R20), which is only a minor source of Hg⁰ in the model.

$$\begin{array}{}\text{(8)}& k_{\mathrm{rd}}=\mathrm{2.92}\times {\mathrm{10}}^{-\mathrm{7}}{e}^{\left(\mathrm{0.045}{T}_{\mathrm{w}}\right)},\end{array}$$

where T_w is the water temperature [^∘C] and k_rd is the dark reduction rate [s⁻¹] (Reaction R1, Table 3).

Cinnabar formation

Additionally, we implemented Hg sulfur chemistry using oxygen concentrations calculated by ECOSMO, whereas sulfur ions (S²⁻) are represented by negative oxygen concentrations in order to reduce the number of transported state variables (Table 1) (Neumann, 2000). In anoxic waters, cinnabar (HgS) is formed by reaction with sulfide species (H₂S, HS⁻, S²⁻) (Reaction R6, Table 3). This reaction is kinetically fast and scavenges the majority of the inorganic Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ within a few hours. The product of this reaction is considered particulate but without a sinking velocity due to the small size of these particles (Paquette and Helz, 1995; Soerensen et al., 2018). In the slower Reaction (R7), HgS subsequently binds to −SH groups of DOM, a reaction that can lead to the dissolution of 50 % of the HgS within 24 h. After 1 d, the dissolution reaction is in equilibrium with the re-crystallization Reaction (R8). In the presence of oxygen, sulfur is quickly oxidized and HgS is readily transformed back into soluble Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ species (HgS_(s)+ 2O₂→ HgSO_4 (aq)) (Reaction R9). In the model, HgSO₄ is attributed back to the dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ pool and not tracked by an additional state variable.

Organic chemistry

The organic chemistry doubles the number of variables introduced for the inorganic Hg chemistry mechanism (Fig. 1). In the model, we implemented three sources for MMHg⁺: (1) methylation in anoxic waters (Reaction R10), (2) methylation in oxic waters (Reaction R11), and (3) methylation due to biologic activity (Reaction R12). The anoxic methylation is thought to be due to anaerobic bacteria and is in our model the fastest methylation process ($\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{7}}$ s⁻¹). Studies have found that methylation also occurs in oxic waters although at much slower rates (Lehnherr, 2014; Heimbürger et al., 2015; Bowman et al., 2020; Soerensen et al., 2018). We implemented an additional constant methylation reaction ($\mathrm{3.47}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻¹) and a biologically induced methylation in oxic water to reflect the fact that numerous bacteria have been shown to actively methylate Hg (Soerensen et al., 2018; Capo et al., 2020). We use the amount of remineralized organic material as a proxy for anoxic microenvironments in the oxic water column. The remineralization is dependent on temperature (Eq. 9), with DOM being mineralized at a higher rate of $k_{{\mathrm{rem}}_{-\mathrm{DOM}}}=\mathrm{10}{k}_{{\mathrm{rem}}_{-\mathrm{POM}}}$. Following Eq. (7) we calculate the amount of remineralized organic matter and use this to scale the biologic methylation rate (Reaction R12). The reaction rate (Reaction R12) has been chosen such that the effective biological methylation rate mostly lies between Reactions (R10) and (R11), ranging from zero to $\mathrm{1.13}\times {\mathrm{10}}^{-\mathrm{7}}$ s⁻¹.

$$\begin{array}{}\text{(9)}& k_{\mathrm{rem}{-}_{\mathrm{POM}}}=\mathrm{0.006}\left\{\left(\mathrm{1}+{\mathrm{20}}^{\left(\frac{T_{\mathrm{w}}^{\mathrm{2}}}{{\mathrm{13}}^{\mathrm{2}}+T_{\mathrm{w}}^{\mathrm{2}}}\right)}\right)\right\},\end{array}$$

where $k_{\mathrm{rem}{-}_{\mathrm{POM}}}$ is the POC remineralization rate [d⁻¹] and T_w is water temperature [^∘C].

Besides MMHg⁺ we also consider double methylation reactions producing DMHg (Reactions R13, R14). For the degradation DMHg → MMHg${}^{+}\to$ Hg²⁺, we consider constant demethylation reactions (Reactions R15, R16), photolytic degradation (Reactions R17–R19), and reductive demethylation (Reaction R20). Finally, we apply methylation and demethylation only to dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ and MeHg${}_{\left(\mathrm{aq}\right)}^{+}$ species. Thus, high loads of DOM and POM influence the effective net methylation and produce a non-linear behavior in the system (Olsen et al., 2018).

Chemical reactions in the sediment

In the sediments, we consider only two species: Hg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+}$ and MMHg${}_{\left(\mathrm{s}\right)}^{+}$. These undergo methylation and demethylation using the same reactions and rates as in the pelagic zone (Table 3). We consider the sediments to always be at least partially anoxic depending on the oxygen concentration in the adjacent water layer (50 %–100 % anoxic for O₂ between 2 and 0 mL L⁻¹). All abiotic methylation reactions (Reactions R10 and R11, Table 3) thus take place in the model sediment. Additionally, Hg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+}$ is subject to dark reduction and subsequently released from the sediment as Hg⁰ (Capo et al., 2022).

2.3.2 Partitioning

The speciation of Hg²⁺ and MMHg⁺ plays a major role in transport, chemical reactions, and bio-availability. In the partitioning scheme we distinguish between three phases: (1) dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ and MeHg${}_{\left(\mathrm{aq}\right)}^{+}$, which are stored in two advected state variables and are further resolved into Hg(OH)_2(aq), HgOHCl_(aq), HgCl_2(aq), MeHgOH_(aq), and MeHgCl_(aq), which in turn are diagnostic variables dependent on salinity; (2) Hg bound to dissolved organic material Hg²⁺–DOM_(aq) and MeHg⁺–DOM_(aq); and (3) the particulate Hg species Hg²⁺–POC_(s) and MeHg⁺–POC_(s).

Three-way partitioning is calculated as a function of Hg concentration, particle load, and dissolved organic matter concentration (Eqs. 10–12). As we could not obtain sorption and desorption rates and because our carbon representation does not capture the number of O- and S-binding sites available for Hg, we implemented partitioning based on partitioning coefficients instead of a dynamic sorption/desorption process as described in Eq. (4). We use a value of log(k_d)=66 for Hg²⁺ associated with DOC based on the work of Tesán Onrubia et al. (2020). This k_d is higher than what is used in other models (Lyon et al., 1997; Zhang et al., 2019; Kawai et al., 2020). Moreover, we use distinct partitioning coefficients for binding to POC k_d and DOC k_l for inorganic Hg²⁺ (log(k_d) = 5.4 L kg⁻¹ and log(k_l)=5.6) and organic MMHg⁺ (log(k_d)=4.9 and log(k_l)=5.0) (Allison and Allison, 2005; Batrakova et al., 2014) (Table 4).

$$\begin{array}{}\text{(10a)}& {\displaystyle }{k}_{\mathrm{d}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{POC}}}{{\mathrm{POCHg}}_{\mathrm{aq}}^{\mathrm{2}+}}}\text{(10b)}& {\displaystyle }{k}_{\mathrm{l}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{DOC}}}{{\mathrm{DOCHg}}_{\mathrm{aq}}^{\mathrm{2}+}}}\end{array}$$

$$\begin{array}{}\text{(11a)}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{aq}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{T}}}{\mathrm{1}+k_{\mathrm{d}}+k_{\mathrm{l}}}}\text{(11b)}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{DOC}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{T}}{k}_{\mathrm{l}}}{\mathrm{1}+k_{\mathrm{d}}+k_{\mathrm{l}}}}\text{(11c)}& {\displaystyle }{\mathrm{Hg}}_{\mathrm{POC}}={\displaystyle \frac{{\mathrm{Hg}}_{\mathrm{T}}{k}_{\mathrm{d}}}{\mathrm{1}+k_{\mathrm{d}}+k_{\mathrm{l}}}}\end{array}$$

Here, k_d is the Hg–POM_(s) $/$ Hg_(aq) partitioning coefficient [1], k_l is the Hg–DOM_(aq) $/$ Hg_(aq) partitioning coefficient [1], POM denotes particulate organic matter [1], SPM denotes suspended particles [1], and DOM denotes dissolved organic matter [1].

The model assumes instantaneous equilibrium and redistributes Hg²⁺ and MeHg⁺ between the three states on each time step. This approach is supported by lab studies that indicate the partitioning equilibrium is reached within an hour (Mason et al., 1995). Finally, mass conservation is ensured by Eq. (12).

$$\begin{array}{}\text{(12)}& {\mathrm{Hg}}_{\mathrm{T}}={\mathrm{Hg}}_{\mathrm{aq}}+{\mathrm{Hg}}_{\mathrm{POC}}+{\mathrm{Hg}}_{\mathrm{DOC}}\end{array}$$

Table 4 Physical and biological constants used in MERCY v2.0.

2.3.3 Radiation

The radiation available for photolytic reactions is determined from hourly input fields using shortwave radiation reaching the surface as modeled by the meteorological model COSMO-CLM (Table 1). As the reaction rates for Hg photolysis are usually reported in relation to photolytically active radiation (PAR), we convert the modeled shortwave radiation using an average factor of 0.5211, not taking into account diurnal variations (Jacovides et al., 2004). We then calculate the cumulative light extinction E_tot (Eq. 13) by water (Eq. 14), phytoplankton (Eq. 15), dissolved organic matter (Eq. 16), and suspended particles (SPM) (Eq. 17), whereby we estimate the total particulate matter concentration for light attenuation using a constant ratio of 0.1 times the particulate organic carbon (POC) concentration (Sharif et al., 2014) (Eq. 18). Finally, the remaining radiation R_z at half the depth of each layer is calculated following the Lambert–Beer law (Eq. 19). All parameters used to calculate light extinction are given in Table 4.

$$\begin{array}{}\text{(13)}& {\displaystyle }{E}_{\mathrm{tot}}=E_{\mathrm{phy}}+E_{\mathrm{POC}}+E_{\mathrm{DOC}}+E_{\mathrm{water}}\text{(14)}& {\displaystyle }{E}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\sum _{z=\mathrm{0}}^n\mathrm{0.05}(h_{z+\mathrm{1}}-h_z)\text{(15)}& {\displaystyle }\begin{array}{rl}{E}_{\mathrm{phy}}& =\sum _{z=\mathrm{0}}^n\mathrm{0.000377}(C_{\mathrm{FLA}}+C_{\mathrm{DIA}}+C_{\mathrm{CYA}})\\ & \cdot (h_{z+\mathrm{1}}-h_z)\end{array}\text{(16)}& {\displaystyle }{E}_{\mathrm{DOC}}=\sum _{z=\mathrm{0}}^n\mathrm{0.00029}{C}_{\mathrm{DOC}}(h_{z+\mathrm{1}}-h_z)\text{(17)}& {\displaystyle }{E}_{\mathrm{POC}}=\sum _{z=\mathrm{0}}^n{\displaystyle \frac{\mathrm{0.0002}{C}_{\mathrm{POC}}}{\mathrm{PSR}}}(h_{z+\mathrm{1}}-h_z)\text{(18)}& {\displaystyle }{C}_{P_{\mathrm{total}}}=C_{\mathrm{POC}}/\mathrm{PSR}\text{(19)}& {\displaystyle }{R}_{\mathrm{z}+\mathrm{1}}=R_z{\exp }^{E_{{\mathrm{tot}}_z}}\end{array}$$

Here, C_FLA is the flagellate concentration [mg C m⁻³], C_DIA is the diatom concentration [mg C m⁻³], C_CYA is the cyanobacteria concentration [mg C m⁻³], C_DOC denotes dissolved organic carbon [mg C m⁻³], C_POC denotes particulate organic carbon [mg C m⁻³], C_P total denotes the total particle load [mg C m⁻³], PSR denotes the POC fraction of total particles [1] = 0.1 (Sharif et al., 2014), E_phy denotes extinction by phytoplankton [1], E_DOC denotes extinction by DOC [1], E_POC denotes extinction by POC [1], $E_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ denotes extinction by water [1], E_tot denotes total light extinction [1], z is the number of vertical layers [1], n is the number of layers [1], h is the height of grid cell z [m], and R is radiation at layer z+1 [W m⁻²].

2.3.4 Biological uptake

Hg bioaccumulation has been implemented directly into the HAMSOM-ECOSMO framework (Daewel and Schrum, 2013; Daewel et al., 2019). ECOSMO is based on a functional group approach lumping species based on properties like nutrient requirements (NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, PO${}_{\mathrm{4}}^{\mathrm{3}+}$, SiO₂) and feeding habits (herbivorous, omnivorous, carnivorous). ECOSMO includes three phytoplankton species (flagellates, diatoms, and cyanobacteria) and two zooplankton species (micro- and mesozooplankton), as well as a macrobenthos and a fish group with the latter representing mass fluxes to higher trophic levels (Fig. 2).

In MERCY we consider bioaccumulation of inorganic Hg²⁺ and organic MeHg⁺ for each of the seven functional groups. Moreover, we distinguish between passive uptake directly from the water column (bioconcentration) and active uptake due to the consumption of contaminated food (biomagnification). The first is accumulated as Hg attached to the organism (zooplankton carapace, fish gills) and the second incorporated internally. Figure 3 depicts a schematic overview of the rate constants used to describe bioaccumulation in MERCY with phytoplankton, which only undergoes passive uptake, on the left and with higher trophic species, which also actively feed on other species, in the middle and on the right. All bioaccumulation processes are calculated separately for inorganic Hg²⁺ and organic MeHg⁺, and the accumulated Hg is transported consistently with the movement of the associated biota. In total, this leads to 22 bioaccumulation state variables (6 phytoplankton, 8 zooplankton, 4 macrobenthos, and 4 fish), which roughly doubles the number of chemical state variables (20) in the model (Table 2). All parameters used for bioaccumulation modeling are given in Table 5.

Figure 2 Overview of the ECOSMO marine ecosystem nutrient and functional group model (Daewel et al., 2019).

Figure 3 Schematic overview of Hg²⁺ and MMHg⁺ bioaccumulation for phytoplankton (left), microzooplankton (middle), and mesozooplankton (right). Dashed lines indicate passive uptake and release rates (Eq. 19); solid lines indicate active uptake due to feeding, with a fraction being instantly released back into the water column (Eq. 22); and dotted lines show Hg loss due to mortality (Eq. 23).

Bioconcentration

In MERCY dissolved Hg${}_{\mathrm{aq}}^{\mathrm{2}+}$ and MMHg${}_{\mathrm{aq}}^{+}$ are accumulated via passive uptake U_P (Eq. 20) through the cell membrane of the phytoplankton functional groups (diatoms, flagellates, cyanobacteria). For zooplankton, macrobenthos, and fish, the passive uptake is thought to lead to Hg accumulation on the surface or areas that are exposed to water like the mouth or gills in the case of fish (Fig. 3: r_u). The uptake rate is r_u calculated based on the surface area A_(B) dependent on an ecosystem functional group B and a Hg-species-dependent permeation velocity v (Eq. 21). We estimate average volume and surface areas for phytoplankton species based on observations of size and geometric shape (Table S1) (Olenina et al., 2003). The cell volume is used to estimate the organic carbon content, which is then used to estimate the ratio of organic carbon to the cell surface (Menden-Deuer and Lessard, 2000). This ratio allows us to model the total phytoplankton cell surface per functional group based on the organic carbon content as modeled in ECOSMO. The estimated surface area is used to calculate the Hg-species-dependent uptake rate based on Mason et al. (1996). Diffusive uptake by zooplankton is implemented based on experimental uptake studies but is less important compared to phytoplankton due to the comparably small surface areas of these species (Tsui and Wang, 2004).

$$\begin{array}{}\text{(20)}& U_{\mathrm{P}\left(B\right)}=r_{\mathrm{u}}{\mathrm{Hg}}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+},\end{array}$$

where U_P (B) is the passive uptake of ecosystem group B [ng s⁻¹], r_u is the passive uptake rate [s⁻¹], and Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ is dissolved Hg [ng m⁻³].

$$\begin{array}{}\text{(21)}& r_{\mathrm{u}}=v_C{A}_B{C}_B,\end{array}$$

where v_C is the permeation velocity for Hg species i [m s⁻¹], A_B is the average surface area of ecosystem group B [m² mg⁻¹ C], and C_B is the concentration of ecosystem group B [mg C m⁻³].

Biomagnification

For all non-phytoplankton species, we consider the active uptake U_A due to feeding rates r_B,b and r_b,B which lead to a fraction ϵ_(C) of the Hg in prey to be incorporated into the predator (Fig. 3: r_B,b r_b,B). Through this process, Hg²⁺ and MMHg⁺ are magnified along the food web (Eq. 22). Zooplankton feeds on detritus, phytoplankton, and other zooplankton, while fish feed on mesozooplankton and macrobenthos following Daewel et al. (2019) (Fig. 4). Moreover, there is macrobenthos that exists only in the marine bottom layer and feeds on these species. We base our uptake on studies that show that only a fraction of Hg²⁺ (ϵ_(Hg)=0.45) and a fraction of MMHg⁺ (ϵ_(MeHg)=0.97) are incorporated into the predator, while the rest is excreted directly back into the water column (Mason et al., 1996; Wang and Wong, 2003; Tsui and Wang, 2004; Pickhardt et al., 2006).

$$\begin{array}{}\text{(22)}& U_{\mathrm{A}\left(B\right)}=\sum _{b=\mathrm{0}}^m{r}_{(B,b)}{\mathit{ϵ}}_{\left(C\right)}{\mathrm{Hg}}_{\left(b\right)}-r_{(b,B)}{\mathrm{Hg}}_{\left(B\right)},\end{array}$$

where U_A (B) is the active uptake rate in ecosystem group B [ng s⁻¹], Hg_B is the Hg concentration in ecosystem group B [ng m⁻³], Hg_b is the Hg concentration in ecosystem group b [ng m⁻³], m is the number of ecosystem groups [1], r_B,b is the feeding rate [m³ s⁻¹], r_B,b is the predation rate [m³ s⁻¹], and ϵ_C is the feeding efficiency [dimensionless between 0–1].

Figure 4 Flowchart of Hg bioaccumulation due to feeding following the ECOSMO end-to-end functional group approximation (Daewel et al., 2019). Rates for all depicted flows are given in Table 5.

Release

Mercury accumulated by active U_A and passive uptake U_P can also be released back into the water column (Eq. 23). There are three distinct processes in the bioaccumulation model that release Hg accumulated in the food web back into the water column. Firstly, there are species-dependent fixed release rates for Hg adsorbed to (r_r1) and adsorbed in (r_r2) the biological species (Eq. 25). Secondly, upon feeding described by feeding rates r_B,b and r_b,B, a fraction 1−ϵ_(C) of the Hg accumulated in prey is not incorporated into the predator, and this is directly released back into the water column (Eq. 24). Finally, based on the ECOSMO mortality and respiration rates r_m for each ecosystem group, Hg is released (Eq. 26). Feeding, mortality, and respiration rates are directly taken from ECOSMO (Table 1), and the relevant equations are described in detail in Daewel et al. (2019). For detritus, the mortality rate is a temperature-dependent remineralization rate (Eq. 9).

$$\begin{array}{}\text{(23)}& R_{\left(B\right)}=R_{\mathrm{R}\left(B\right)}+R_{\mathrm{F}\left(B\right)}+R_{\mathrm{M}\left(B\right)},\end{array}$$

where R_(B) is the release rate from ecosystem group B [ng s⁻¹], R_R(B) (Eq. 24) is the constant release rate [ng s⁻¹], R_F(B) (Eq. 25) is the feeding-related release rate [ng s⁻¹], and R_M(B) (Eq. 26) is the mortality-related release rate [ng s⁻¹].

$$\begin{array}{}\text{(24)}& R_{\mathrm{R}\left(B\right)}=\sum _{B=\mathrm{0}}^m{r}_{\mathrm{r}\mathrm{1}\left(B\right)}{\mathrm{Hg}}_{\left(B\right)\mathrm{ext}}+r_{\mathrm{r}\mathrm{2}\left(B\right)}{\mathrm{Hg}}_{\left(B\right)\mathrm{int}},\end{array}$$

where Hg_(b) ext is Hg on ecosystem group B [ng m⁻³], Hg_(b) int is Hg inside ecosystem group B [ng m⁻³], r_r1 is the release rate of external Hg [m³ s⁻¹], and r_r2 is the release rate of internal Hg [m³ s⁻¹].

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{R}_{\mathrm{F}\left(B\right)}& =\sum _{B=\mathrm{0}}^m\mathit{\left\{}\sum _{b=\mathrm{0}}^m{r}_{(B,b)}\right[(\mathrm{1}-{\mathit{ϵ}}_{\left(C\right)\mathrm{ext}}){\mathrm{Hg}}_{\left(b\right)\mathrm{ext}}\\ & +(\mathrm{1}-{\mathit{ϵ}}_{\left(C\right)\mathrm{int}}){\mathrm{Hg}}_{\left(b\right)\mathrm{int}}\left]\mathit{\right\}},\end{array}\end{array}$$

where r_B,b is the feeding rate of group B on group b [s⁻¹], ϵ_C int is the feeding efficiency for external Hg species C [dimensionless between 0–1], ϵ_C ext is the feeding efficiency for external Hg species C [dimensionless between 0–1], Hg_(b) ext is Hg on ecosystem group b [ng m⁻³], and Hg_(b) int is Hg inside ecosystem group b [ng m⁻³].

$$\begin{array}{}\text{(26)}& R_{\mathrm{M}\left(B\right)}=\sum _{B=\mathrm{0}}^m{r}_{\mathrm{m}\left(B\right)}({\mathrm{Hg}}_{\left(B\right)\mathrm{ext}}+{\mathrm{Hg}}_{\left(B\right)\mathrm{int}}),\end{array}$$

where r_m is the mortality rate of ecosystem group B [m³ s⁻¹].

Table 5 Overview of bioaccumulation parameters. External variables taken from the ecosystem model ECOSMO such as mortality (r_m) and feeding rates (r_f) are given in Table 1. The abbreviated phytoplankton references are (A) Mason et al. (1996), Menden-Deuer and Lessard (2000), and Olenina et al. (2003) and (B) Pickhardt and Fisher (2007) and Nfon et al. (2009).

2.3.5 Benthic–pelagic coupling

Following the sediment modeling concept by Daewel et al. (2019), we implemented a simple two-layer sediment system, where the first layer interacts with the lowest water column grid cell and the second layer represents a permanent sink.

Sedimentation

Sedimentation occurs due to the settling of Hg bound to particles and detritus. The sedimentation flux F_s is calculated using a sinking velocity w_d of 5 m d⁻¹ for Hg bound to particles (POC) (Daewel and Schrum, 2013) (Eq. 27).

$$\begin{array}{}\text{(27)}& F_{\mathrm{s}}=w_{\mathrm{d}}{\mathrm{Hg}}_{\mathrm{POC}}^{\mathrm{2}+},\end{array}$$

where F_S is the sedimentation flux [ng s⁻¹ m⁻²], Hg${}_{\mathrm{POC}}^{\mathrm{2}+}$ is the particulate mercury concentration in water [ng m⁻³], and w_d is the sinking velocity [m s⁻¹].

Resuspension

Resuspension F_r is triggered by a critical ocean current velocity of U 0.01 m s⁻¹. In the case that a critical current velocity is reached, no sedimentation takes place and a resuspension rate r_res of 25 [d⁻¹] is used to release Hg${}_{\left(\mathrm{s}\right)}^{\mathrm{2}+}$ from the first sediment layer into the lowest water grid cell (Eq. 28). Depending on the depth (<1 m) of the lowest grid cell and current velocity, resuspension can also directly affect the second-lowest water grid cell.

$$\begin{array}{}\text{(28)}& F_{\mathrm{r}}=r_{\mathrm{res}}{\mathrm{Hg}}_{\mathrm{S}}^{\mathrm{2}+},\end{array}$$

where F_r is the resuspension flux [ng s⁻¹ m⁻²], r_res is the resuspension rate [s⁻¹], and Hg${}_{\mathrm{S}}^{\mathrm{2}+}$ is the mercury concentration in sediment [ng m⁻²].

Burial

Hg²⁺ and MMHg⁺ in the first layer are constantly transported to the second layer, which represents a permanent sink in the model. The burial flux F_b is based on a constant burial rate of k_bur= $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ [d⁻¹] (Eq. 29, Table 4).

$$\begin{array}{}\text{(29)}& F_{\mathrm{b}}=r_{\mathrm{bur}}{\mathrm{Hg}}_{\mathrm{S}}^{\mathrm{2}+},\end{array}$$

where F_b is burial flux [ng s⁻¹ m⁻²], r_bur is the burial rate [s⁻¹], and Hg${}_{\mathrm{s}}^{\mathrm{2}+}$ is the mercury concentration in sediment [ng m⁻²].

2.3.6 Air–sea exchange

Air–sea exchange of elemental Hg⁰ is one of the most important processes in the global Hg cycle. Here, we use the approach of Kuss et al. (2009), Kuss (2014)), which is based on Henry's law constant H by Andersson et al. (2008) to determine the equilibrium between Hg⁰ in water ${\mathrm{Hg}}_{\left(\mathrm{aq}\right)}^{\mathrm{0}}$ and air ${\mathrm{Hg}}_{\left(\mathrm{air}\right)}^{\mathrm{0}}$ (Eq. 30). Next, the transfer velocity for CO₂k₆₀₀ is approximated using a quadratic parameterization depending on the 10 m wind speed U₁₀ (Eq. 31). We then calculate the transfer velocity k_w for Hg⁰ by scaling k₆₀₀ using the temperature (T)-dependent and salinity (S)-dependent diffusivity of Hg⁰ in water (Eqs. 32 to 35) (Kuss, 2014). The actual inter-compartmental Hg⁰ flux F_Hg is then calculated based on surface concentrations in the adjacent compartments (Eq. 36). The air–sea exchange is also applied for DMHg. However, the CMAQ-Hg model does not consider DMHg yet. Hence, the atmosphere is only a sink for DMHg, which is instantaneously transformed into Hg²⁺ (Niki et al., 1983), and its fate is currently not explicitly resolved.

$$\begin{array}{}\text{(30)}& H_{\mathrm{Hg}}=e^{\left(\frac{-\mathrm{2404.3}}{T}+\mathrm{6.915}\right)}\end{array}$$

(Anderssen et al., 2008)

$$\begin{array}{}\text{(31)}& k_{\mathrm{600}}=\mathrm{0.222}{U}_{\mathrm{10}}^{\mathrm{2}}\cdot \mathrm{333}{U}_{\mathrm{10}}\end{array}$$

(Nightingale et al., 2000)

$$\begin{array}{}\text{(32)}& {\mathit{Sc}}_{\mathrm{35}}=-\mathrm{0.0398}{T}^{\mathrm{3}}+\mathrm{3.3910}{T}^{\mathrm{2}}-\mathrm{118.02}T+\mathrm{1948.2}\end{array}$$

(Kuss, 2014)

$$\begin{array}{}\text{(33)}& {\mathit{Sc}}_{\mathrm{0}}=-\mathrm{0.0304}{T}^{\mathrm{3}}+\mathrm{2.7457}{T}^{\mathrm{2}}-\mathrm{118.13}T+\mathrm{2226.2}\end{array}$$

(Kuss, 2014)

$$\begin{array}{}\text{(34)}& \mathit{Sc}={\displaystyle \frac{{\mathit{Sc}}_{\mathrm{35}}S+{\mathit{Sc}}_{\mathrm{0}}(\mathrm{35}-S)}{\mathrm{35}}}\end{array}$$

(Kuss, 2014)

$$\begin{array}{}\text{(35)}& k_{\mathrm{w}}=k_{\mathrm{600}}\sqrt{{\displaystyle \frac{\mathit{Sc}}{\mathrm{600}}}}/\mathrm{360}\mathrm{000}\end{array}$$

(Kuss, 2014)

$$\begin{array}{}\text{(36)}& F_{\mathrm{Hg}}={\displaystyle \frac{{\mathrm{Hg}}_{\left(\mathrm{air}\right)}^{\mathrm{0}}-H_{\mathrm{Hg}}\cdot {\mathrm{Hg}}_{\left(\mathrm{aq}\right)}^{\mathrm{0}}}{k_{\mathrm{w}}}}\end{array}$$

(Schwarzenbach et al., 2003)

In the above equations (Eqs. 30–36), H_Hg is Henry's law constant [1], T is water temperature [^∘C], S is salinity [PSU], Sc₃₅ is the Schmidt number for salt water [1], Sc₀ is the Schmidt number for fresh water [1], k₆₀₀ is the transfer velocity of CO₂ [cm h⁻¹], k_w is the transfer velocity of Hg [cm h⁻¹], Hg${}_{\left(\mathrm{air}\right)}^{\mathrm{0}}$ is the Hg⁰ concentration in air [ng m⁻³], Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{0}}$ is the Hg⁰ concentration in water [ng m⁻³], and F_Hg is the net Hg⁰ flux from atmosphere to water [ng m⁻² h⁻¹].

2.3.7 Technical implementation

As a basis for the presented model development, we build upon the setup used for the earlier-published inorganic marine Hg model MERCY (Bieser and Schrum, 2016). All processes are implemented as stand-alone routines which are called from a main driver function containing several time loops (Fig. 5). Data for the wet cells (pelagic) are stored in vector form to reduce overhead, and data for sediments (benthic) and the lowest atmospheric layer are stored in 2D fields. Input data (Table 1) are read directly during run time from binary ECOSMO output as hourly mean values. This approach was chosen because there is no feedback from the Hg chemistry on the physical and biological models and because it allows us to reduce the computational costs of running the marine Hg model. All output files are created with daily mean values and saved in netCDF format using the IO-API interface (Byun and Schere, 2006; IO-API). The model is set up in a way that it runs for a single year using the last output time step of the previous year as an initial condition. For this initial model evaluation, we run MERCY for 17 years from 2000 to 2016.

Figure 5 Schematic overview of the MERCY model routines and main time loop.

3 Model evaluation

We determine the model performance in reproducing observed concentrations and dynamics (e.g., variability and seasonality) of individual Hg species. Based on this analysis, we identify the processes and parameters responsible for the model error. The model is not specifically calibrated to the area of application, the North and Baltic seas. It is built on the current understanding of mercury cycling in the ocean and should be generally applicable. Major factors that need to be considered before applying the MERCY model to other regions are (1) partitioning coefficients for organic material (OM) as the type of OM varies regionally, (2) the parameterization for biogenic reduction as the values presented here are based on cyanobacteria in the Baltic Sea, (3) the uptake and release rates for bioaccumulation which might not be representative of other regions, and (4) the ecosystem model used that is needed to drive MERCY.

3.1 Statistics

Because there are no established quality criteria for marine models, we use criteria commonly used for evaluation of atmospheric CTMs (Derwent et al., 2010; Thunis et al., 2012, 2013; Carnevale et al., 2014). We start by comparing the observed and predicted means (Eq. 37) using daily model averages in the 10×10 km² grid cell corresponding to the observation. As statistical metrics, we use bias (Eq. 38), error (Eq. 39), standard deviation (Eqs. 40, 41), and the correlation coefficient (Eq. 42) to evaluate systematic error, random error, amplitude error, and phase error. However, for most Hg species, the observations lack the temporal coverage to determine the phase error. Moreover, we use the centered root-mean-square error (CRMSE) because it allows us to distinguish between systematic error (bias) and random error (CRMSE) (Eq. 43) (Carnevale et al., 2014). For our analysis, we normalize the statistical metrics to get concentration-independent values and allow for better comparability between different Hg species.

Equation (37) gives the mean.

$$\begin{array}{}\text{(37)}& O={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{0}}^N{O}_{i}P={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{0}}^N{P}_i,\end{array}$$

where P_i is the predicted value from the model, O_i denotes observed values from measurement, N is the sample size/number of observations, and i is the index.

Equation (38) gives the normalized mean bias.

$$\begin{array}{}\text{(38)}& \mathrm{NMB}={\displaystyle \frac{P-O}{O}}\end{array}$$

Equation (39) gives the normalized centered root-mean-square error.

$$\begin{array}{}\text{(39)}& \mathrm{NCRMSE}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{N}\sum _{i=\mathrm{0}}^N\mathit{\left\{}\right(O_i-O)-(P_i-P){\mathit{\}}}^{\mathrm{2}}}}{O}}\end{array}$$

Equation (40) gives the standard deviation.

$$\begin{array}{}\text{(40)}& {\mathit{\sigma }}_o=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}(O_i-O{)}^{\mathrm{2}}}{\mathit{\sigma }}_p=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}(P_i-P{)}^{\mathrm{2}}}\end{array}$$

Equation (41) gives the normalized mean standard deviation.

$$\begin{array}{}\text{(41)}& \mathrm{NMSD}={\displaystyle \frac{{\mathit{\sigma }}_P-{\mathit{\sigma }}_O}{{\mathit{\sigma }}_O}}\end{array}$$

Equation (42) gives the correlation coefficient (r).

$$\begin{array}{}\text{(42)}& r={\displaystyle \frac{\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^N(O_i-O)(P_i-P)}{{\mathit{\sigma }}_O{\mathit{\sigma }}_P}}\end{array}$$

Equation (43) gives the root-mean-square error.

$$\begin{array}{}\text{(43)}& {\mathrm{RMSE}}^{\mathrm{2}}={\text{CRMSE}}^{\mathrm{2}}+{\text{bias}}^{\mathrm{2}}\end{array}$$

Finally, we use additional quality criteria to determine model performance: firstly, the percentage of model values within a factor of 2 (FAC2), which gives an easy-to-understand estimate of the model quality (Eq. 44). We argue that model values within a factor of 2 are within the combined uncertainty. The uncertainty consists of the measurement uncertainty, the sampling uncertainty when comparing observations with time-integrated (24 h) and space-integrated (100 km²) model grid cells (Schutgens et al., 2016), and the error propagation in the biogeochemical modeling framework. We estimate the measurement error U to be in the range from 20 % for Hg⁰ and Hg_T to 50 % for MeHg.

Equation (44) gives the factor of 2.

$$\begin{array}{}\text{(44)}& \mathrm{FAC}\mathrm{2}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{n}_i,\end{array}$$

with n_i=1 for $\mathrm{5}<\frac{P_i}{O_i}<\mathrm{2}$ and otherwise n_i=0.

Secondly, we use the more technical model quality objective (MQO) as defined by Carnevale et al. (2014). The MQO (Eq. 45) relates the root-mean-square error (Eq. 46) to the root-mean-square uncertainty (Eq. 47). The MQO can be interpreted as follows: for MQO<0.5 on average the model values lie within the measurement uncertainty and thus the model cannot be improved upon unless more precise observations become available. For MQO>1 the model error is on average larger than the measurement uncertainty but the model may be closer to the “true” environmental value than the observations. Thus, the aim is to achieve an MQO<1. Moreover, we determine model performance criteria for NMB, NMSD, and RMSE as proposed by Carnevale et al. (2014) (Eqs. 49–51).

Equation (45) gives the model quality objective.

$$\begin{array}{}\text{(45)}& \mathrm{MQO}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{RMSE}}{{\mathrm{rms}}_U}}\end{array}$$

Equation (46) gives the root-mean-square error.

$$\begin{array}{}\text{(46)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(P_i-O_i{)}^{\mathrm{2}}}\end{array}$$

Equation (47) gives the root-mean-square uncertainty.

$$\begin{array}{}\text{(47)}& {\mathrm{rms}}_U=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{\mathrm{i}=\mathrm{0}}^N{U}^{\mathrm{2}}},\end{array}$$

where U denotes the measurement uncertainty.

Equation (48) gives the model performance criterion for NMB.

$$\begin{array}{}\text{(48)}& {\text{MPC}}_{\mathrm{NMB}}={\displaystyle \frac{\mathrm{2}U}{O}}\end{array}$$

Equation (49) gives the model performance criterion for NMSD.

$$\begin{array}{}\text{(49)}& {\text{MPC}}_{\mathrm{NMSD}}={\displaystyle \frac{\mathrm{1}U}{{\mathit{\sigma }}_O}}\end{array}$$

Equation (50) gives the model performance criterion for RMSE.

$$\begin{array}{}\text{(50)}& {\text{MPC}}_{\mathrm{RMSE}}\le \mathrm{1.0}\end{array}$$

These quality criteria have been developed for atmospheric pollutants like ozone, nitrogen oxides, and fine particles, which have been studied and modeled for decades. For modeling of Hg in the marine environment, the observations are still very limited compared to those of pollutants in the atmosphere. This affects the ability to use these criteria, and we therefore do not expect the MERCY model to meet the criteria at this point. However, we define these as our future goal for marine Hg modeling.

3.2 Model domain (North Sea and Baltic Sea)

Here, we evaluate the model for the North and Baltic seas in northern Europe (model domain shown in Fig. 6). This area was chosen for model evaluation as it covers a large range of different physical and biological conditions: the Baltic Sea (Fig. 6; marine regions 8–15) is an enclosed shelf sea with a surface area of 377 000 km². It is connected to the North Sea (marine regions 1–5) via the shallow Kattegat and Skagerrak (marine regions 6–7) in the southwest. It is a brackish waterbody strongly influenced by freshwater inflow, and it covers a salinity range from <2 PSU in the north that increases towards the southwest to up to 35 PSU in the transition zone between the North and Baltic seas (Fischer and Matthäus, 1996; Lehmann and Post, 2015; Mohrholz et al., 2015). The central Baltic Sea has several deep basins reaching a depth of 460 m in the Landsort Deep (Fig. 6b). It exhibits strong stable stratification with more saline, in parts anoxic, deep water, resulting from an estuarine circulation system with upper-layer outflow of fresh water and lower-layer saline inflow. Every few years, large quantities of oxic and saline waters are transported from the North Sea to the Baltic Sea during so-called major Baltic inflows (MBIs). During the simulation period 2000 to 2015, three MBIs occurred; one of these was an especially strong event during the winter of 2014/15 (Fischer and Matthäus, 1996; Lehmann and Post, 2015). In the northern part of the Baltic, the Bothnian Sea and the Bothnian Bay are seasonally covered by ice, possibly leading to accumulation of Hg from rivers during winter due to the suppression of Hg⁰ evasion. Finally, the Baltic Sea is a system with cyanobacteria, which makes it an interesting study area as these cyanobacteria have been shown to actively reduce Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ (Kuss et al., 2015). Moreover, cyanobacteria can lead to pronounced early-spring–late-summer biomass blooms that affect bioaccumulation (Soerensen et al., 2016a).

The North Sea has a surface area of 575 000 km² and is connected to the Atlantic Ocean at its northern border and via the English Channel to the south. It is a shallow shelf ocean that is well mixed during autumn and winter, and it experiences frequent resuspension events during autumn and winter storms. The southern North Sea is characterized by strong tidal mixing, and thus water masses are well mixed and sediments are resuspended regularly within the tidal cycle. It is an area of high primary productivity and an important fishing ground. Thus, it is an important study area for Hg methylation and bioaccumulation.

Due to their close vicinity to the coast and national monitoring programs, there are a comparably large number of Hg observations available for both the North Sea and the Baltic Sea. However, the data on Hg are still sparse in some areas, especially regarding Hg speciation, which is a major obstacle for model evaluation.

Figure 6 (a) COSMO-CMAQ-Hg atmospheric model domain with the North and Baltic seas highlighted. (b) MERCY marine model domain and topography. (c) Marine regions: (1) English Channel, (S) Scheldt Estuary, (2) southern North Sea, (3) northern North Sea, (4) Norwegian Trench, (5) German Bight, (6) Kattegat, (7) Skagerrak, (8) Belt Sea, (6–8) Swedish west coast, (9) Arkona Sea, (10) Bornholm Sea, (11) Gotland Sea or Baltic Proper, (12) Bay of Riga, (13) Neva Bay, (14) Bothnian Sea, and (15) Bothnian Bay.

3.2.1 Forcing data

To generate the necessary forcing data (Table 1) to run the MERCY model, we used the four models described in Sect. 2.1. For the atmosphere, COSMO-CLM was run on a regional domain for Europe driven by ERA-Interim re-analysis data (Berrisford et al., 2011; Dee et al., 2011; Hersbach et al., 2020). The atmospheric model domain covers the entire European landmass, including northern Africa and western Russia, with a resolution of 24×24 km and 35 vertical layers (Fig. 6a). The calculated meteorology is then used as forcing for the atmospheric CTM CMAQ-Hg, which is set up for the same domain and resolution (Byun and Schere, 2006; Zhu et al., 2015; Bieser et al., 2016). CMAQ-Hg uses boundary concentrations for Hg by an ensemble of the global Hg models GEOS-Chem, GLEMOS, ECHMERIT, and GEM-MACH-Hg (Travnikov et al., 2017) and all other relevant trace gases from the global CTM MOZART (Horowitz et al., 2001). Emissions for the year 2010 were created with the SMOKE-EU emission model (Bieser et al., 2011). Hg emissions are based on the AMAP emission inventory (AMAP/EMEP, 2019b). This is a similar setup to that used in previous studies (Bieser and Schrum, 2016). For computational reasons, we calculated only 1 year (2010) of atmospheric Hg concentration and deposition fields. These were used as boundary conditions for the marine Hg model for all years of the simulation. The ocean-ecosystem model HAMSOM-ECOSMO was run on a model domain covering the entire Baltic Sea and the North Sea with open boundaries in the English Channel and at 63^∘ N, where the North Sea is connected to the Atlantic (Fig. 6b). The resolution of the model is about 10×10 km² (spherical grid) with 20 layers and a maximum water depth of 630 m. The vertical resolution is 5 m for the four uppermost layers with a bottom layer depth of 250 m.

3.2.2 Initial conditions

As initial conditions, we interpolated observations in water, biota, and sediment using a traditional kriging methodology to produce realistic initial starting conditions (mostly the pronounced vertical gradient) and minimize the spin-up time required (Cressie, 1990). The observational Hg data were retrieved from the database of the German Federal Maritime and Hydrographic Agency (MARENET, 2020). We ran the model using initial conditions multiplied by factors of 0.5 and 2.0 and tested the time necessary for the two runs to converge. For our model domain, which is a relatively small and in parts enclosed shelf sea area, the model runs started to converge after only a few years in the water column but took several years for Hg in sediments and biota (especially at higher trophic levels). For this study, we used a spin-up time of 30 years to reach realistic initial conditions for the production runs.

3.2.3 Boundary conditions

The chosen domain, including only the North Sea and Baltic Sea, has only a very small open boundary: the English Channel in the southwest, which forms a narrow connection to the Atlantic Ocean, and the wider opening in the northern channel. The North Sea in the north of the domain, which receives most of the Atlantic inflow, is connected to the open Atlantic Ocean at the shelf break. This region is characterized by an outflow in the eastern part and inflow in the western part. At the open boundaries, we prescribe constant Hg concentrations using 1.0 pM Hg_T for the North Atlantic and 3.0 pM Hg_T in the English Channel (Cossa et al., 2018; Leermakers et al., 2001).

3.2.4 River loads

River loads are taken from OSPAR and HELCOM reports and the Norwegian Tilførsel program (Green et al., 2011; HELCOM, 2007, 2011). We implemented rivers using monthly load data in the North Sea and annual values for the Baltic Sea as described in Bieser and Schrum (2016). The annual inflow of Hg through rivers is 1100 kg a⁻¹ for the Baltic and 2800 kg a⁻¹ for the North Sea. In the North Sea, the largest fluxes are from the Elbe (1160 kg a⁻¹) and Rhine (800 kg a⁻¹) rivers.

3.2.5 Deposition of Hg²⁺ and atmospheric Hg⁰

Dry and wet Hg deposition is read in as hourly totals from CMAQ netCDF output files. The deposited Hg${}_{\left(\mathrm{g}\right)}^{\mathrm{2}+}$ and Hg${}_{\left(\mathrm{p}\right)}^{\mathrm{2}+}$ are added to the dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$ species assuming instant dissolution of atmospheric particles. In CMAQ, the exchange of Hg⁰ is set to zero for all grid cells with the land-use category water to avoid a doubling of the air–sea exchange calculation in the atmospheric model.

3.3 Observational data

For the model performance, we start by evaluating total Hg (Hg_T) concentrations in the water column. We then look at the individual species, elemental Hg⁰ and organic MeHg. Next, we evaluate the model skill in reproducing Hg concentrations in biota. For this, we compare Hg and MeHg loads in phytoplankton and zooplankton and finally total Hg in fish (Hg_Fish).

Figure 7 Annual averages: (a) Hg_T concentrations in the top 50 m with superimposed observations (Kuss et al., 2017; Soerensen et al., 2018; MARENET, 2020), (b) Hg_T concentrations below 50 m with superimposed observations, (c) average fraction of filterable Hg_Filt. as a proportion of Hg_T, and (d) intra-annual variability in modeled daily average Hg_T concentrations.

3.3.1 Total Hg in water

The available Hg_T observations cover offshore and coastal areas in the North and Baltic seas. Hg_T has been measured as both unfiltered (Soerensen et al., 2018) and the filterable fraction Hg_Filt. (Kuss et al., 2017; MARENET, 2020). In MERCY, Hg_Filt. corresponds to the sum of eight species, namely Hg${}_{\left(\mathrm{g}\right)}^{\mathrm{0}}$, Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$, Hg–DOM_(aq), HgS–DOM_(aq), MMHg${}_{\left(\mathrm{aq}\right)}^{+}$, MMHg–DOM_(aq), DMHg_(g), and HgS_(s) (Table 2). Hg_T additionally includes Hg_P, which is comprised of the two particulate species Hg–POC_(s) and MeHg–POC_(s). In our model Hg_Filt. makes up about 95 % of Hg_T on average (Fig. 7c). Hg_P only plays a significant role (>5 % on annual average) in the southern North Sea and the Wadden Sea. Especially in the Wadden Sea, observed Hg_T concentrations are extremely high with values ranging from 6 to 117 pM. For the model performance evaluation, we removed measurements taken in the Wadden Sea and other areas not resolved in our model setup (e.g., the area between the coastline and barrier islands and lagoons in the Baltic Sea). As depicted in Fig. 7a, virtually no observations are from regions where particles play a major role. Thus, for simplification we will use the term Hg_T to refer to all of these observations, but we compare them to concentrations of the respective model species.

In the North Sea we use 435 measurements of Hg_T sampled between 2007 and 2011 (MARENET, 2020). The samples are taken over the entire year, but BSH (the German Federal Maritime and Hydrographic Agency) sampling focuses mostly on the German exclusive economic zone, although it also includes a few years with data for the greater North Sea. Finally, all measurements are surface samples (0–20 m), which is due to the shallow nature of the North Sea. For the Baltic Sea, there are 111 observations from the MARENET database (MARENET, 2020); 168 observations from three cruises in March 2014, March 2015, and July–August 2015, which cover the southern part of the Baltic Sea from the Belt Sea to the Gotland Deep (Kuss et al., 2017); and 90 observations from three cruises in July and August of 2015 and 2016, which in addition include observations on the Bothnian Sea and Bothnian Bay (Soerensen et al., 2018). Figure 7a and b depict all Hg_T observations used for model evaluation.

3.3.2 Individual marine species: Hg⁰ and MMHg⁺

For the evaluation of Hg⁰, we use 580 measurements from four Baltic Sea cruises in February, April, July, and November 2006 (Kuss, 2014). This dataset allows us to investigate the seasonality of redox reactions. For MMHg⁺, we were able to obtain 310 measurements from six cruises in 2014, 2015, and 2016 covering coastal and offshore areas of the Baltic Sea (Kuss et al., 2017; Soerensen et al., 2017, 2018). For 160 of these, both MeHg and Hg_T were available, which enables a relative evaluation of the methylated fraction M_frac= MeHg $/$ Hg_T. For the North Sea, no Hg⁰ or MMHg⁺ observations are available at all with the exception of nine MeHg measurements from 1999 in the English Channel and the Scheldt Estuary, which we used to set the MMHg⁺ boundary conditions in the English Channel (Leemakers et al., 2001). Thus, we are forced to limit the model evaluation for Hg⁰ and MMHg⁺ to the Baltic Sea.

3.3.3 Hg in biota

Bioaccumulation in the marine biota is evaluated by comparing their total Hg and MeHg content to measured concentrations in biota in the Baltic Sea (Nfon et al., 2009). For evaluation of fish total Hg, we use Hg_T concentration in the muscle tissue of 1166 herring from coastal and offshore locations in the Baltic Sea (Soerensen and Faxneld, 2020). As the biota measurements are in wet weight and our model is in milligrams of carbon dry weight, the ratio of milligrams of carbon per milligram biota of 0.2 for diatoms, 0.33 for flagellates and cyanobacteria, and 0.5 for zooplankton and fish was assumed (Sicko-Goad et al., 1984; Walve and Larsson, 1999). This was combined with a conversion factor of dry weight to wet weight of 0.2 for phytoplankton, 0.16 for zooplankton, and 0.1 for fish (Cushing, 1958; Ricciardi and Bourget, 1998). For phyto- and zooplankton, the model is compared to the observed average, minimum, and maximum concentrations, but due to limited data, no seasonal or regional comparison was possible. For fish we analyze Hg accumulation for five Baltic Sea regions ranging from the western Baltic to the Bothnian Bay.

3.4 Model performance

3.4.1 North Sea (Hg_T)

Figure 8 compares the frequency distribution of 435 Hg_T measurements to the associated model values. It can be seen that the majority of observations lie between 1 and 3 pM, which is captured well by the model. However, the observed high values between 5 and 10 pM cannot be reproduced. We argue that these are due to input from the coastal area (e.g., major rivers, Wadden Sea) not included as input to the model in the current river discharge scenario.

Hg_T concentrations in the North Sea do not exhibit a pronounced seasonality, and the observed variability is driven by a strong land–sea gradient along the European coastline where Hg from rivers is transported northeastwards from the English Channel by the Coriolis force (Fig. 7d). For the analysis, we split the North Sea Hg observations into two groups: (1) the Elbe Estuary (N=366) and (2) the open North Sea including a few observations near the remaining coastline (N=69). Due to the significant Hg inflow from the Elbe (1164 t a⁻¹), Hg concentrations are higher in the Elbe Estuary with a mean concentration of 3.44 pM (Table 6). The model is able to reproduce the observed average (NMB $=-\mathrm{21}$ %) but has a better agreement with the median values (−12 %). In this region, random and amplitude errors are dominant. This is indicative of subgrid dynamics and our inability to resolve the seasonality of Hg from rivers stemming from the use of static river loads for the entire run (OSPAR Commission et al., 2016). However, with 70 % of model values within a factor of 2 of the observations and an MQO of 1.48, the model is still close to our quality goal.

In the less dynamic open North Sea, the model performs better (FAC2 = 84 %, MQO = 1.22) (Table 6). The observed average of 1.92 pM is matched by the model (2.03 pM), and the bias is close to zero (NMB = 6 %). Nevertheless, due to the inhomogeneous distribution of observations, this value is not indicative of the actual background Hg concentrations in the open North Sea. We find that Hg concentrations there are mostly in the range of 1.1 to 1.5 pM.

Figure 8 Frequency distribution of observed (red) and model (blue) Hg_T. Concentrations for the North Sea (N=435).

In summary, for Hg_T, the model is close to our quality criterion (MQO ≤1.0). Improvements to the MQO could likely be achieved by increasing model resolution in the complex coastal regions and including more detailed input from rivers and particle resuspension at the European coastline. Especially for the Wadden Sea, a hydrodynamical model that can model the intertidal zone would be preferential.

Table 6 Regional model performance for Hg_T in the North Sea. The evaluation is based on 435 measurements from the MARENET database (MARENET, 2020). Obs: observations; Mod: model.

3.4.2 Baltic Sea (Hg_T)

In the Baltic Sea, model performance for Hg_T is similar to in the North Sea (FAC2 = 70 %, MQO = 1.28) with a low bias (NMB = −19 %) and a high random error (NCRMSE = 102 %) (Table 7). Unlike for the North Sea, the model predicts a pronounced seasonality with surface Hg_T concentrations around 50 % higher during March (Fig. 9a) compared to August (Fig. 9b), which is in line with observations from Kuss et al. (2017) taken in March and July–August (Fig. 9). The two processes governing this are (1) stratification and particle settling in the central Baltic deep basins after the onset of primary production – this is the biological pump as POC particles here are mainly of biological origin (detritus) – and (2) increased photoreduction and subsequent atmospheric exchange of Hg⁰ (air–sea exchange). Additionally, during winter higher atmospheric Hg⁰ concentrations due to heating-related emissions and a shallow planetary boundary layer reduce and sometimes even reverse the Hg⁰ air–sea gradient. In the open Baltic Sea, Hg concentrations are mostly between 1.0 and 1.5 pM. In stratified areas, Hg_T concentrations can drop down to 0.5 pM during summer. During autumn and winter, mixing and upwelling can occasionally transport Hg from deeper waters upwards, sometimes leading to surface concentrations above 2 pM in some areas.

Figure 9 Hg_T surface concentrations for (a) March and (b) August. Dots indicate observations (Kuss et al., 2017; Soerensen et al., 2018; MARENET, 2020).

Figure 10 Surface transect of the Hg_T concentrations in the Baltic Sea. The x axis gives 365 daily averages for the year 2010; the y axis represents a transect from the German coast in the western Baltic (y=0) to the mouth of the Bothnian Sea (y=87).

Table 7 Spatially aggregated observed and modeled Hg_T concentrations in the Baltic Sea.

For a more detailed analysis, we separate the Baltic Sea into three regions: (1) the western part, which includes the Belt, Arkona, and Bornholm seas; (2) the Gotland Sea in the central Baltic; and (3) the northern part which includes the Bothnian Sea and Bothnian Bay. Moreover, we evaluate the oxic surface/intermediate waters and the deep anoxic waters in the Gotland area separately (Table 7). It is seen that the model is able to reproduce surface concentrations in the western and central areas with a bias close to zero. The model bias is larger in the deep basins, but model performance is still comparable to the North Sea. Here, the low vertical resolution in the model setup below 100 m will play a role. In the northern part, the model strongly overestimates Hg_T concentrations. This overestimation was also seen in the Soerensen et al. (2018) model. Northern Baltic rivers tend to be low in POC but rich in DOC compared to temperate rivers (McClelland et al., 2016; Soerensen et al., 2017), highlighting the importance of DOC flocculation at the point where river water encounters higher-salinity water for the settling and removal of Hg in Bothnian Bay estuaries, something that is currently not included in our model.

Figure 10 depicts the seasonality and Fig. 11 three vertical profiles in the Gotland Sea. It is seen how quickly Hg concentrations can change in this region and, depending on physical drivers, how different the seasonality of vertical mixing can be. At location A (Gotland Deep) Hg concentrations are around 1.5 pM for most of the year with a strong surface depletion (1 pM) during August and September. At location C, located at the opposite side of Gotland, the seasonality is reverse with the highest concentrations (1.2–1.4 pM) during August and September and much lower concentrations (0.9–1.1 pM) throughout the rest of the year.

In summary, our conclusion is similar to that of the North Sea, i.e., that better data on Hg inputs from rivers and a better resolution of the physical processes in the domain seem the most promising options for improving model performance. Especially in the Bothnian Bay, Hg cycling seems to be strongly influenced by terrestrial organic matter. In the central Baltic, we found that typically used k_d values around log(k_d) = 5.5 are not sufficient to reproduce the observed Hg depletion in the surface waters. Here, as described in Sect. 2.3.2, we use a k_d based on Tesán Onrubia et al. (2020) which is an order of magnitude higher and leads to improved correlation to observations (Table 4). In addition, a higher vertical resolution is advised as vertical transport has proven to be an issue with the current model setup. Due to the low model resolution below 150 m, numerical diffusion leads to an overestimation of mixing in the deep basins. Finally, for further model evaluation, it would be useful to increase the seasonal coverage of observations in the area.

Figure 11 Vertical Hg profiles in the central Baltic Sea observations (red) (Soerensen et al., 2018) and model values (blue) for the three central Baltic deep basins given in Fig. 9.

Finally, as the deep basins of the Baltic Sea are anoxic, in this area sulfur chemistry becomes relevant (Reactions R6–R9, Table 3). The effect of adding HgS and HgS–DOM to the chemistry scheme leads to particulate Hg–POM transforming into dissolved HgS species. The effect of this is two-fold: (1) firstly, Hg that is scavenged from the stratified surface layer by detritus (biological pump) accumulates directly at the boundary between oxic and anoxic waters. (2) Secondly, as eventually all inorganic Hg is transformed into HgS species, particle settling stops being a sink and Hg persists in the water column, whereas Hg is effectively transported to the sediment in model runs without sulfur chemistry. This leads to Hg concentrations being constant in the anoxic layer with higher values found only directly at the seafloor. Comparing to observations, we find that the model with sulfur chemistry is better able to capture the observed Hg distribution (Soerensen et al., 2018).

3.4.3 Elemental mercury (Hg0)

In the marine environment, elemental Hg⁰ makes up only a few percent of the total Hg_T. However, it is the species that determines the air–sea exchange and thus is a major driver for atmospheric long-range transport. With the oceans being the largest Hg emitter into the atmosphere (roughly twice as large as current anthropogenic emissions), marine Hg⁰ determines global transport patterns. Moreover, errors in modeled Hg⁰ concentrations propagate to all other Hg species and lead to wrong estimates for the compartmental budgets. Thus, it is of utmost importance to correctly reproduce Hg⁰ concentrations in surface waters. A detailed model study on Hg air–sea exchange in the North and Baltic seas has been published using a previous model version (Kuss, 2014; Kuss et al., 2018; Bieser and Schrum, 2016). The four main drivers of Hg⁰ concentrations are as follows:

1. the reducible fraction of Hg²⁺, which is typically estimated to be 40 % of the dissolved Hg${}_{\left(\mathrm{aq}\right)}^{\mathrm{2}+}$;

2. the parameterization of biologically induced reduction processes;

3. the modeled photon flux and wavelength-dependent extinction in water impacting photolytic reduction;

4. air–sea exchange parameterizations, especially during high wind speeds.

Due to the fast exchange between atmosphere and water, Hg⁰ concentrations converge towards the equilibrium as described by Henry's law constant (Andersson et al., 2008). Therefore, in shelf seas a change in the redox chemistry directly affects the total Hg^T in the system. Due to the mixing in the coastal ocean, this impacts almost the complete waterbody. Moreover, the different reduction pathways produce a distinct seasonal pattern with Hg⁰ concentrations ranging from as low as 5 pg L⁻¹ during winter to up to peaks >60 pg L⁻¹ during cyanobacteria blooms. Thus because of the high intra-annual variability, the model needs to be evaluated against Hg⁰ observations throughout the year, as good agreement for a single cruise does not imply good model performance throughout the year.

Figure 12 Comparison of modeled and observed Hg⁰ concentrations in surface waters for four cruises in the Baltic Sea (x axis: observations; y axis: model values) (Kuss, 2014; Kuss et al., 2018).

Table 8 Comparison of modeled and observed Hg⁰ concentrations for four cruises in the Baltic Sea in 2006 (Kuss, 2014; Kuss et al., 2018).

The observed annual average Hg⁰ concentration for 580 measurements is 14.6 pg L⁻¹; the modeled value is 14.9 pg L⁻¹ with a systematic error of 2 % and a random error of 35 %. The random error is largest in summer (NCRMSE = 42 %) and is due to the biogenic reduction which depends on cyanobacteria biomass. The model shows good correlation with observations (R=0.60) and is able to reproduce the observed variability (NMSD = 38 %) (Table 8). All statistical metrics for Hg⁰ are well inside model performance criteria, and 97 % of the model values are within a factor of 2 of the observations. The model quality objective is below 1.0 (MQO = 0.84). Thus, the model performance is, at least for the given model resolution, in the range where further improvements are hardly feasible (Carnevale et al., 2014; Schutgens et al., 2016).

We acknowledge that the redox chemistry used is based on measurements in the Baltic Sea (Kuss et al., 2015). Thus, it needs to be investigated whether it shows equally good performance for other marine regions. We find that the model performs similarly well throughout the year with the largest bias during summer, when the dynamics driving biological and photolytic reduction lead to a higher variability in Hg⁰ concentrations (Table 9).

Table 9 Seasonal breakdown of Hg⁰ model performance.

Figure 13 depicts the seasonality of a mean Hg⁰ for the Baltic Sea. Moreover, the contributions of the four reduction reactions (1) chemical reduction, (2) photolytic reduction, (3) biogenic reduction, and (4) reductive methylation (Table 3) are given. We find that the dark reduction is the dominant process, producing 55 % of the total Hg⁰ in the Baltic Sea and 70 % in the North Sea. Photolytic reduction contributes 34 % and biogenic reduction contributes 12 % annually. However, from July to mid-August the photolytic reduction becomes dominant (>50 %). When the cyanobacteria bloom starts, light penetration reduces significantly due to the increased marine particle load, and until the end of November the biogenic reduction becomes the dominant process (Fig. 13a). In contrast, as there are no cyanobacteria in the North Sea, photolytic reduction is dominant throughout the summer (Fig. 13b). The reductive methylation reaction plays a negligible role for Hg⁰ surface concentrations but can be a source for Hg⁰ in deeper waters with a high MeHg fraction. It can be seen that there is a background Hg⁰ concentration of about 5 to 15 pg L⁻¹ due to the chemical (dark) reduction process. During model development, we recognized a systematic error in the seasonality (overestimation during summer and underestimation during winter) that could be resolved by introducing a temperature dependency of the chemical reduction reaction, a process which was detectable in the observations by Kuss et al. (2015) (Eq. 8, Sect. 2.3.1). For the photolytic reaction, we found that it is important to validate the radiation fields. Testing the model using different radiation fields resulted in a change in the annual net Hg⁰ production of >10 %. The main driver here is cloud coverage, which is a particularly uncertain state variable in meteorological models. Moreover, we want to note that photolysis rates from observations and incubation experiments are solely reported based on the photolytically active radiation. Due to the highly wavelength-selective light extinction it would be favorable to parameterize photolysis using the actual wavelengths absorbed by Hg. Finally, the biogenic reduction term in the model is driven only by the concentration of cyanobacteria. This creates the observed late-summer–early-autumn Hg⁰ peak. Allowing other phytoplankton species in the model to reduce Hg²⁺, leading to unrealistically high concentrations during the spring flagellate bloom.

Figure 13 Annual profile of mean Hg0 concentration in the Baltic Sea (a) and North Sea (b). The colored areas indicate the contribution of individual reduction pathways (Reactions R1, R3, R20; Table 3).

3.4.4 Methylmercury (MMHg⁺ and DMHg)

Due to the complexity of the analytical methods and the extremely low environmental levels of observed concentrations, MeHg observations in the marine environment are rare. Additionally, they are the most uncertain observations. Here, to calculate the MQO, we assume an uncertainty of 50 %. We evaluate the model predictive capabilities in reproducing (1) MeHg concentrations and (2) the methylated Hg fraction M_frac= MeHg $/$ Hg_T. The latter allows us to evaluate the modeled net methylation independently of the Hg_T model error (Table 10). In the Baltic Sea, observed MeHg concentrations are in the range of 191 (20–603) fM, while the modeled range is 213 (57–350) fM. For M_frac the observed range is 11.4 % (1.3 %–30.4 %) (10th and 90th percentiles), while the model predicts 9.9 % (3.6 %–20.2 %). The frequency distribution for observed and modeled M_frac is given in Fig. 15. The model is in very good agreement with the observations on average but cannot reproduce the observed extreme values. In total there are 17 (6.5 %) samples with a M_frac between 33 % and 73 %, all of which were measured in the intermediate layer between 70 and 160 m.

Evaluating the relative M_frac metric instead of absolute MeHg concentrations reduces the systematic error from −28 % to 5 % and the amplitude error from −74 % to −55 %. This shows that the Hg_T bias accounts for roughly 80 % of the MeHg systematic error and 50 % of the amplitude error. Yet, using M_frac has no significant effect on the random error, indicating a non-linear relationship between the methylated fraction and the absolute amount of MeHg. While systematic error and amplitude error are comparable to the other Hg species, the random error is much larger (NCRMSE = 1.9). This shows that the methylation–demethylation dynamics in the model are too simplified, pointing to missing processes in the model. Figure 15 depicts MeHg and M_frac vertical profiles for the central Baltic Sea deep basins in different years and seasons together with oxygen concentrations. Again, it can be seen that the model is able to reproduce the average vertical profiles but is incapable of capturing the high and low values. Observations indicate an MeHg hotspot near the oxycline. Here, M_frac can become as large as 100 %, meaning that almost all mercury there is present as MeHg. The highest MeHg observations coincide with anoxic conditions, indicating that the availability of dissolved HgS drives methylation in these regions (Soerensen et al., 2018). In Fig. 15, anoxic regions are indicated by negative oxygen concentrations. These are based on measurements of H₂S, and the net oxygen is calculated based on the reaction H₂S + 2O₂→ H₂SO₄.

Figure 14 Observed (Kuss et al., 2017; Soerensen et al., 2018) and modeled frequency distribution of the methylated Hg fraction M_frac.

Table 10 Evaluation of seasonally and vertically clustered M_frac observations against model values.

Figure 15 Vertical MeHg profiles for Baltic deep basins. Negative oxygen concentrations indicate sulfide concentrations. (a) MeHg profiles 2014 (Soerensen et al., 2018); (b) MeHg profiles 2015 (Soerensen et al., 2018); (c) M_frac March (Kuss et al., 2017); (d) M_frac July–August (Kuss et al., 2017).

The model can reproduce the seasonality and vertical gradient of the methylated fraction. On the one hand photolytic demethylation leads to lower MeHg concentrations in the surface ocean during summer. On the other hand, biological activity leads to increased MeHg formation in spring and summer. We find that a biologically induced methylation parameterized with biomass or phytoplankton concentration leads to spring becoming the season with the most effective net methylation. By linking biological methylation to the remineralization of organic carbon, we introduce a temperature dependency that shifts this towards summer (Fig. 16) (Eq. 9, Sect. 2.3.1). Yet, the model still overestimates methylation in spring and underestimates methylation in summer. For a more detailed analysis, we look at surface layer MeHg concentrations on four specific dates. Figure 17 depicts MeHg measurements for 21 March and 1 August of the years 2014 and 2015. In March MeHg concentrations are between 40 and 300 fM and in August between 10 and 200 fM with pronounced spatial gradients. This “spottiness” of the MeHg concentrations partially explains the large random error in the model. Moreover, while the general patterns are similar, methylation shows a significant interannual variability (Fig. 17).

Figure 16 Seasonality of the biologically induced methylation reaction using different parameterizations (Reaction R12, Table 2).

Figure 17 Methylmercury concentrations in the surface ocean on (a) 21 March 2014, (b) 1 August 2014, (c) 21 March 2015, and (d) 1 August 2015. Superimposed are all observations at depths of 0–50 m (Table 9).

Overall, the model reproduces 53 % of MeHg values within a factor of 2. We find that the model performance (MQO = 0.98) is still within the quality criterion. This is due to the much higher uncertainty in MeHg measurements, for which we assumed an error of 50 %. This indicates that further model improvement will be difficult unless more frequent and more precise MeHg measurements become available. Moreover, to reach their full potential MeHg observations need to be combined with extensive auxiliary data. This starts with simple parameters like incoming solar radiation to determine the actual intensity of photolysis or better estimates for special partitioning coefficients for MeHg. In our model, for example, we use a lower k_d of MMHg⁺ compared to Hg²⁺, which means that particle settling will increase M_frac with increasing depth (Table 4, Sect. 3.2.3). Moreover, chemical parameters such as O₂ and H₂S concentrations have been shown to impact the availability of inorganic Hg²⁺ species for methylation. And finally, microbiological observations ranging from chlorophyll concentration to RNA show the activity of methylating bacteria could improve variable methylation rates. From our model evaluation, it seems clear that fixed methylation and demethylation rates cannot account for the observed variability in both MeHg concentrations and the methylated Hg fraction M_frac. Here, we need a better understanding of the parameters modulating methylation and demethylation rates.

3.4.5 Hg in biota

Figure 18 depicts annual average Hg loads in the different ecosystem biota species. The North Sea exhibits higher Hg loads in biota, which can be explained by the high Hg load from rivers, especially the Elbe and Scheldt; the lack of permanent sedimentation; and the earlier onset and higher overall primary production, which increases the effectiveness of the active uptake pathway. The average amount of Hg in biota ranges from 1 % to 5 % of the Hg_T with higher values in the highly productive North Sea. During winter only a little Hg is bound in biota due to the low biomass, while in summer the fraction can be up to 10 %. Due to the high transfer efficiency of MMHg⁺ (97 %), on average, between 5 % and 20 % of the total MMHg⁺ is accumulated in biota. In highly productive regions the amount of MMHg⁺ inside biota can even be larger than the MMHg⁺ remaining in the water column. Flagellates (Fig. 18b) are the most abundant phytoplankton species and thus the most important primary accumulator. However, the diatom bloom occurs earlier in the year and removes MeHg from the water column before the flagellate bloom. The higher Hg load in diatoms (Fig. 18a) is due to their lower carbon content. Finally, cyanobacteria (Fig. 18c), which can lead to major blooms in late summer–early autumn, are the dominant species later in the year, and MeHg loads during the bloom exceed those of the diatoms and flagellates. Due to the active Hg uptake, micro- (Fig. 18d) and mesozooplankton (Fig. 18e) have a higher accumulation factor than the phytoplankton species. Finally, Fig. 18f depicts the Hg load in fish.

Figure 18 Annual average Hg_T concentrations [ng L⁻¹] in biota: (a) diatoms, (b) flagellates, (c) cyanobacteria, (d) microzooplankton, (e) mesozooplankton, and (f) fraction of Hg in fish.

Figure 19 Seasonality of modeled (a) Hg and (b) MeHg loads in phytoplankton and zooplankton. Superimposed is the water concentration of (a) Hg and (b) MeHg. Panel (c) gives the integrated biomass. All values are averages for the Baltic Sea integrated over the top 100 m.

As the last step of the model evaluation, we compare Hg_T and MeHg loads in biota to observations. Field studies investigating the total Hg content of biota are fairly common and can be used to estimate the model bias. However, only few data on MeHg in biota and species diversity are available. On average the observed Hg_T content in phytoplankton lies in the range of 0.002±0.001 µg g⁻¹ wet weight and that for zooplankton in the range of 0.006±0.005 µg g⁻¹ wet weight (Nfon et al., 2009). Here, we use a conversion factor for wet weight (w.w.) to dry weight (d.w.) of 0.2 for phytoplankton and 0.16 for zooplankton (Cushing, 1958; Ricciardi and Bourget, 1998). Moreover, biomass in ECOSMO is defined in milligrams of carbon, whereas observations are reported in milligrams w.w. or milligrams d.w. We estimate the ratio of milligrams of carbon to milligrams d.w. as 0.2 for diatoms, 0.33 for flagellates and cyanobacteria, and 0.5 for zooplankton (Sicko-Goad et al., 1984; Walve and Larsson, 1999). With this, we estimate the expected average Hg_T loads in biota in the Baltic Sea in the range of 30 (15–45) ng g⁻¹ C in phytoplankton and 75 (10–120) ng g⁻¹ C in zooplankton. MeHg loads in phytoplankton are expected to be around 10 (5–15) ng g⁻¹ d.w., while they are larger for cells with a larger surface-to-volume ratio (Pickhardt and Fisher, 2007; Soerensen et al., 2016b). Figure 19 depicts the average Hg_T and MeHg loads in phytoplankton and zooplankton. For phytoplankton, the model lies well within the expected range for Hg_T (25–80 ng g⁻¹ C) (Fig. 19a) and MeHg (5–15 ng g⁻¹ C) (Fig. 19b). During winter when phytoplankton biomass is low (Fig. 19c), Hg loads reach the maximum of the expected bioaccumulation range, and once production starts, growth dilution lowers the modeled Hg_T and MeHg loads and their concentrations in the water column decline by 10 % and 20 % respectively. For zooplankton, values are well within the expected range. They start low at the beginning and rise over the year. At the end of the year, cyanobacteria start to dominate the phytoplankton community, leading to higher phytoplankton MeHg concentrations and another decline in marine Hg_T concentrations. At the same time, the decrease in diatom and flagellate concentration leads to an increase in the fraction of microzooplankton of the mesozooplankton diet, increasing their trophic level and further increasing the zooplankton Hg_T load. Finally, with Hg loads between 70 and 140 ng g⁻¹ C d.w., fish are the highest trophic level in the ecosystem model. Due to the much more efficient active uptake of MMHg⁺ compared to Hg²⁺ in fish, 60 % to 80 % of the accumulated Hg is in the form of MMHg⁺. Looking at the two uptake pathways of bioconcentration and biomagnification, we find that biomagnification is responsible for 80 % to 90 % of the total Hg uptake for non-phytoplankton species. A more detailed analysis can be found in Amptmeijer et al. (2023).

Next, we evaluate the model capabilities to reproduce Hg content in fish. For this, we compare the modeled bioaccumulation in the functional ecosystem group representing fish to herring. This pelagic species corresponds best to the fish functional group implemented into ECOSMO (Daewel et al., 2019). The analysis is based on 1166 measurements of Hg in fish muscle tissue. We use the same conversion factors as for zooplankton to convert the model carbon dry weight to wet weight total biomass (1 ng g⁻¹ C d.w. = 3.125 ng g⁻¹ w.w.). For this, the dataset is split into five Baltic Sea regions:

1. the Swedish west coast, a stripe from Gothenburg to Oslo;

2. the southern Baltic Proper, which includes the Bornholm Sea and the southern Gotland Sea;

3. the northern Baltic Proper, which includes most of the Gotland Sea;

4. the Bothnian Sea;

5. the Bothnian Bay.

It is not possible to compare the caught fish to an individual model grid cell and time step. Therefore, we compare them to observed average Hg_Fish concentrations. The model reproduces the observed average Hg_Fish of 28 ng g⁻¹ in the Baltic Sea with a systematic error of −9 % (Hg_Fish=25 ng g⁻¹) (Table 11). In order to estimate the model variability in Hg_Fish for each region, we vertically integrate annual average model values for each grid column. The result is a fish dataset in which each member represents a model fish that has spent its life in a single 10×10 km² water column. In reality, herring are not confined to 10×10 km² grid cells, and their Hg accumulation depends on their migration patterns. Yet, we argue that this approach approximates the model spread (Fig. 20). This allows us not only to calculate the bias but also to estimate the model standard deviation. On average, over the whole Baltic Sea, the model captures the observed variability (NMSD = 9 %). The error is driven by the west coast region (NMSD = 309 %), while it varies between (29 %–76 %) in the remaining Baltic Sea. In the west coast region, the observed fish exhibit less than half the variability observed in all other regions. While the model captures the variability in the other regions, it shows the opposite behavior to the observations on the west coast. In this shallow region, we explain the high model concentrations by regular Hg resuspension from sediments, which creates pockets of elevated Hg_T and MeHg concentrations. Thus, the large model spread is an artifact of our methodology based on static fish.

Figure 20 Modeled and observed frequency distribution of Hg in fish in the Baltic Sea regions (Soerensen and Faxneld, 2020).

4 Conclusions

In this paper, we present the regional-scale 3D high-resolution biogeochemical multi-media Hg model MERCY v2.0. The numerical model combines hydrodynamic models for the atmosphere and ocean, including a marine ecosystem model. MERCY includes a comprehensive marine Hg scheme to calculate transport, transformation, and bioaccumulation. The schemes for chemistry, partitioning, and bioaccumulation are based on literature values, and no domain-specific model tuning has been done. We would like to emphasize that MERCY is suitable for any marine region or even for global application. The major factors when applying the MERCY model to other regions are (1) partitioning coefficients to organic material (OM) as the type of OM varies regionally; (2) the parameterization for biogenic reduction as the values presented here are based on cyanobacteria in the Baltic Sea; and (3) the ecosystem model, as trophic dynamics and phytoplankton uptake rates can vary widely between regions. To our knowledge, it is the first model capable of linking atmospheric Hg emissions to MeHg accumulation at higher trophic levels. The intention of this initial model publication is the detailed presentation of the model and first results, focusing on model performance evaluation and the identification of the processes and parameters responsible for the model error. A more comprehensive analysis of the dynamics of and variability in Hg speciation, partitioning, and bioaccumulation is required for future studies. While our model performs more realistically than earlier models for marine Hg cycling, there are still large uncertainties, especially regarding methylation.

Table 11 Model performance evaluation of Hg_T, dissolved elemental Hg⁰, and the methylated Hg fraction. The model error is separated into systematic error (normalized bias), random error (normalized centered root-mean-square error), and amplitude error (normalized mean standard deviation). The model quality objective target value is MQO <1.0 (consult Sect. 3.1 for more information).

n/a: not applicable.

We evaluated model performance for key Hg species based on a simulation for the North and Baltic seas for the years 2000 to 2016. We chose these regions due to the availability of observations. Moreover, the two regions cover a range of regimes, have high primary productivity, and are relevant to fisheries. Unlike atmospheric Hg modeling, there is no precedent or scientific consensus defining the state-of-the-art requirements and limitations of reproducing concentrations of different marine Hg species. Considering the inherent uncertainty in a comparison of model values and observed concentrations (e.g., measurement error, sampling error, error in the hydrodynamic models, the uncertainty in reaction rates, and unknown processes), we define model values within a factor of 2 of the observations as a reasonable agreement. Moreover, we used a statistical model quality objective (MQO <1.0) to assess the model skill (Carnevale et al., 2014) (Sect. 3.1).

A detailed model performance evaluation for the North and Baltic seas demonstrates that the model can reproduce concentrations and seasonality of single Hg species to a degree that validates the model predictive capabilities. For Hg_T, the model is able to reproduce 72 % of the observations within a factor of 2 (Table 12). We find that the model can reproduce background concentrations in the open parts of the shelf seas (1.0–1.5 pM). The model error can mostly be attributed to random and amplitude error. The main source of uncertainty in the model is the transport dynamics of the large Hg influx from rivers and the Wadden Sea. These lead to observed Hg peaks of up to 10 pM. The model resolution of 10×10 km proved insufficient to reproduce the observed temporal and spatial gradients. Because the majority of observations are at the coast near major rivers, the model does not reach the quality objective (MQO = 1.44). Moreover, in the Baltic Sea, the model overestimates vertical mixing from deeper regions with elevated Hg concentrations. This is caused by the coarse vertical resolution below 150 m, which leads to numerical diffusion and an underestimation of stratification. We found that including sulfur chemistry improves model performance in the deep anoxic water layer in the Baltic Sea basins. The mechanism is that Hg transported downwards from the stratified oxic and productive surface layer through the biological pump transforms into dissolved HgS species in anoxic waters. This stops the downward gradient and lessens the role of the sediments in this region as a sink.

We summarize that the improvement of the model performance for Hg_T requires optimizing of the hydrodynamic model. Unless circulation patterns, stratification seasonality, resuspension events, and upwelling regions are correctly represented, hardly any improvement of the model can be achieved. Further, for the coastal ocean, we find that river inflow needs to be better resolved, ideally with daily loads including fluxes of dissolved and particulate carbon. Finally, particle partitioning and subsequent sedimentation comprise a major source of uncertainty. We achieved better results using a log(k_d) of 6.6 (Tesán Onrubia et al., 2020), which is an order of magnitude higher than those used by other models.

The model performed best for elemental Hg⁰. Due to air–sea exchange, Hg⁰ is the key species controlling the exchange between the atmosphere and ocean. Any bias in modeled Hg⁰ fields directly influences the marine total Hg budget and leads to unrealistic results. MERCY is able to reproduce 97 % of Hg⁰ measurements within a factor of 2. We find that the chemical (often referred to as dark) and photolytic reduction processes produce roughly the same amount of Hg⁰ annually although with different seasonality. Moreover, elevated Hg⁰ concentrations in the Baltic Sea between July and October could be reproduced by implementing biological reduction by cyanobacteria. Finally, we find that it is important to consider temperature dependence for the chemical reduction reaction to correctly reproduce the seasonality. With a model skill of MQO = 0.84, we conclude that the model performance for Hg⁰ is in a range where further improvements become marginal. Possible improvements are photolytic reaction rates based on actual wavelengths instead of the photolytic active radiation and a better understanding of biological reducers.

Evaluation of MeHg resulted in the methylated fraction M_frac (MeHg $/$ Hg_T) for which 55 % of model values are within a factor of 2 of observations. The model is able to reproduce the observed mean and seasonality but is unable to capture the observed maxima, resulting in a large random error. Yet, because of the high measurement uncertainty, the model still reaches the quality objective (MQO = 0.98), indicating that the observations are limiting model development. We found that producing realistic MeHg concentrations throughout the year required methylation occurring in oxic waters. Oxic methylation is the primary or sole source (80 %–100 %) of MeHg in large parts of the model domain. The anoxic methylation reaction is dominant in anoxic waters (the deep basins of the Baltic Sea). We found that assumptions made in other models linking methylation to productivity or chlorophyll concentrations pose two problems: firstly, they lead to regions with zero MeHg in seasons with no primary production and very low MeHg concentrations in the deep anoxic basins. And secondly, they produce a phase error in the seasonality due to an overestimation of MeHg during the spring bloom. In MERCY, we parameterize the biogenic methylation with the amount of remineralized organic matter, which adds a temperature dependence to the process that in turn reduces the impact of the spring bloom. Moreover, various sensitivity runs using varying parameters to modulate the biogenic methylation rate to test for possible biological drivers have failed to surpass model formulations including a constant oxic methylation reaction. We summarize that poor model performance for MeHg is the key source of uncertainty in the presented model. In order to improve the model performance, a more detailed understanding of methylation processes is required. Moreover, more high-quality observations, especially on MeHg seasonality, are needed to allow for model-based process studies. The addition of isotopic fractionation to the model might also help to further constrain sources and sinks of MeHg.

Finally, we evaluate the model's ability to reproduce Hg in biota. Our model provides Hg and MeHg loads in phytoplankton, zooplankton, and fish which are inside of the observed range. We find that the modeled phytoplankton concentrations vary within the observed maximum and minimum loads. Zooplankton changes at the trophic level over the course of the year due to changes in diet. As expected, the model predicts the highest MeHg loads in fish, making up 90 % of the total Hg in fish due to its high transfer efficiency. Most parameters used for bioaccumulation are highly uncertain, and there is ample room for improvement in this part of the model. We hypothesize that the ecosystem model which is focused on correctly reproducing carbon fluxes, needs improvements regarding functional traits relevant to bioaccumulation such as size, shape, or feeding behavior.

The presented model allows hypothesis testing within a consistent physical–biological–biogeochemical framework based on basic principles. We are currently working on a model version that allows for seamless coupling with different hydrodynamic ocean and marine ecosystem models to increase the applicability of the model. The model performance is here only cursorily evaluated to limit the length of the paper. For the future, we plan to investigate the sources of model uncertainty and sensitivity in order to identify the insufficient understanding of the processes and find out the imprecise or unknown parameters, especially concerning methylation and biological uptake. Finally, we want to employ and promote the MERCY model as a tool for hypothesis testing and prediction within a consistent physical–biological–biogeochemical framework based on basic principles. This will enable researchers to (1) improve our understanding of the natural variability from seasonal to decadal timescales; (2) investigate forcing dynamics, leading to MeHg accumulation in seafood; and (3) estimate the impact of anthropogenic and natural drivers in support of the Minamata Convention on mercury.