Modelling the coupled mercury-halogen-ozone cycle in the central Arctic during spring

Each spring halogens released from snow, sea ice, and sea salt aerosols cause the depletion of boundary layer ozone (O₃) and mercury (Hg) in the Arctic (Oltmans, 1981; Barrie, 1986; Barrie et al., 1988; Bottenheim et al., 1990; Schroeder et al., 1998; Lindberg et al., 2001; Skov et al., 2004; Simpson et al., 2007; Abbatt et al., 2012; Wang et al., 2019). Depletion occurs when halogens and other radicals in the atmosphere oxidize elemental mercury, Hg(0) to Hg(II), which is much shorter-lived and is taken up directly onto aerosols, snow, ice and the open ocean (Selin, 2009; Douglas et al., 2012). Oxidized mercury present in particles is also removed to ice, snow, and the ocean via both wet and dry deposition. Once deposited to the Arctic Ocean, mercury can undergo transformation to more toxic forms, including methylmercury, that are harmful to ecosystems and human health (Driscoll et al., 2013). The main form of atmospheric Hg is long-lived, gaseous Hg(0), emitted from both anthropogenic (e.g., coal burning and artisanal gold mines) and natural sources (e.g., volcanoes; Lindberg and Stratton, 1998; Arctic Monitoring and Assessment Programme [AMAP], 2011). Anthropogenic Hg sources in the Arctic are estimated to be less than 1% of global Hg emissions to air, with most Arctic Hg(0) originally from emissions from distant sources (Durnford et al., 2010; Skov et al., 2020; Dastoor et al., 2022b). Understanding the fate of Arctic Hg that involves air-surface exchange processes remains a challenge, including integration of mercury from the atmosphere into the cryosphere and ocean, storage of oxidized mercury in snow and ice, and photochemical re-emissions of mercury back to the atmosphere.

Developing models that describe diurnal variation and cycling of mercury in the lowest part of the Arctic atmosphere has remained a challenge because of the complexity of chemical and physical processes involved. In addition, closely coupled chemical cycles, such as halogens, have recently been improved within Arctic atmospheric models (e.g., Yang et al., 2008; Yang et al., 2010; Toyota et al., 2011; Falk and Sinnhuber, 2018; Fernandez et al., 2019; Yang et al., 2019; Marelle et al., 2021; Herrmann et al., 2022; Swanson et al., 2022). Figure 1 illustrates the main chemical processes of Hg in the Arctic boundary layer. During atmospheric mercury depletion events (AMDEs), Hg(0) is oxidized forming gaseous oxidized mercury, Hg(II), and particulate mercury species, Hg(p), which are deposited to the cryosphere via dry and wet processes (Lindqvist and Rodhe, 1985; Brooks et al., 2006; Skov et al., 2006; Ariya et al., 2015). AMDEs are linked to high concentrations of reactive bromine (e.g., Br₂, Br, BrO, HOBr, and BrONO₂), activated from trace amounts of oceanic bromide on snow and ice surfaces (Lu et al., 2001; Brooks et al., 2006; Sommar et al., 2007; Stephens et al., 2012; Wang et al., 2019). This process occurs simultaneously with the well-studied depletion of surface ozone, known as ozone depletion events (e.g., ODEs; Barrie et al., 1988; Platt and Hönninger, 2003; Simpson et al., 2007; Abbatt et al., 2012; Simpson et al., 2015). Upon deposition of Hg to snow, a fraction can be photoreduced and re-emitted back into the atmosphere as Hg(0) (e.g., Poulain et al., 2004; Brooks et al., 2006; Kirk et al., 2006; Skov et al., 2006; Sommar et al., 2007; Ferrari et al., 2008; Dommergue et al., 2010; Durnford and Dastoor, 2011; Steffen et al., 2013; Douglas and Blum, 2019). Chemical transformation between Hg(0), Hg(II), and Hg(p) can occur rapidly and is driven by several factors including oxidant and aerosol concentrations, air temperature, and solar radiation (Ariya et al., 2015). Currently, many models are limited in capturing the behaviour of Hg in Arctic spring due to: (1) uncertainties in the mercury chemical mechanism; (2) use of offline-calculated oxidant concentrations including reactive bromine; (3) a lack of detailed halogen chemistry/emission processes over snow and sea ice; and/or (4) poorly represented Arctic boundary layer meteorology (Angot et al., 2016; Travnikov et al., 2017).

Bromine is a central species in Arctic O₃ and Hg(0) depletion and is considered to be the major oxidant driving AMDEs (Wang et al., 2019). Elevated Br concentrations in the Arctic are known to occur due to the uptake and emission of reactive bromine on surfaces, which occurs via an autocatalytic process on snow, ice, and aerosols known as the “bromine explosion” (Simpson et al., 2007; Simpson et al., 2015). Two main bromine activation mechanisms in the Arctic have been proposed and tested in 3D chemical transport models. The first method, proposed by Toyota et al. (2011), involves the activation of bromide in the top layer of the snowpack, triggered by ozone deposition or heterogeneous recycling of HOBr and BrONO₂ (Toyota et al., 2011; Falk and Sinnhuber, 2018; Fernandez et al., 2019; Herrmann et al., 2021; Marelle et al., 2021; Swanson et al., 2022). The second method involves the release of bromine from sea salt aerosols (SSA), formed via sublimation of lofted blowing snow (as proposed by Yang et al., 2008; Yang et al., 2010; Huang and Jaeglé, 2017; Rhodes et al., 2017; Huang et al., 2018; Yang et al., 2019; Huang et al., 2020; Yang et al., 2020). As both mechanisms have important impacts on reactive bromine and ozone concentrations in the Arctic, they need to be considered in models for springtime chemistry.

Descriptions of mercury in models often differ in their treatment of physical and chemical processes, partly due to current knowledge gaps (Ariya et al., 2015). Multi-model intercomparisons are therefore useful to assess the impacts of different model parameterizations on the behaviour of atmospheric Hg (Angot et al., 2016; Travnikov et al., 2017; Dastoor et al., 2022a). Arctic-focused model intercomparisons have shown that several models are capable of simulating the seasonality in surface Hg(0) concentrations (spring minima and summer maxima); however, they can underestimate the amplitude of the seasonal variation (Angot et al., 2016; Dastoor et al., 2022a). Several factors are likely to contribute to this disparity, which underscores some of the limitations of current global models. First, most global models that include Hg chemistry use monthly averaged oxidant fields to simulate Hg(0) oxidation (Holmes et al., 2010; Durnford et al., 2012; Song et al., 2015; Horowitz et al., 2017; Shah et al., 2021). For the purpose of studying springtime mercury chemistry, this approach is inadequate as it neglects the diurnal variability of oxidants, including bromine, in the Arctic (Angot et al., 2016). Secondly, the low spatial and temporal resolution of global models may not sufficiently resolve the local chemistry and emissions occurring that contribute to AMDEs during spring (Toyota et al., 2014). Bromine production mechanisms from snow and sea ice, which are crucial for polar atmospheric chemistry, can therefore be inaccurately represented. In comparison to global models, regional models are computationally less expensive at high spatial resolution and can calculate oxidant concentrations online, making them ideal tools for studying AMDEs.

We recently updated our Arctic-specific version of the regional Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) model (Marelle et al., 2021) to include the atmospheric mercury cycle described in Shah et al. (2021). This version of WRF-Chem 4.3.3 currently includes detailed bromine emissions from both Arctic surface snow and blowing snow, previously shown to improve model representation of Arctic ODEs (Marelle et al., 2021). In this dedicated Arctic study, we investigated the impacts of polar bromine activation on springtime Hg and O₃ chemistry and depletion, during spring 2020. To evaluate the model, we compared with data from Arctic stations and from the central Arctic obtained on board the Research Vessel (RV) Polarstern during the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) campaign in spring 2020 (Shupe et al., 2022). The new mercury developments and WRF-Chem model setup are presented in Section 2, with the observations used to evaluate the model described in Section 3. In Section 4, we evaluate the model against meteorological and chemical observations from MOSAiC. In Section 5, we assess the model performance compared to other Arctic observations, as well as discuss the simulated pan-Arctic impacts. Finally, the outcomes and conclusions are presented in Section 6.

We extended the work of Marelle et al. (2021) by including Hg gas-phase chemistry and photolysis (Section 2.1.1), heterogeneous uptake to aerosols and clouds (Section 2.1.2), dry and wet deposition processes (Section 2.1.3), and Hg(0) re-emission from continental snow and snow on sea ice (Section 2.1.4). These updates are based largely on the latest developments in other Hg chemical transport models (e.g., Holmes et al., 2010; Amos et al., 2012; Fisher et al., 2012; Fisher et al., 2013; Gencarelli et al., 2014; Horowitz et al., 2017; Shah et al., 2021), which were implemented in the WRF-Chem model here. The model developments presented in this study are publicly available as Ahmed et al. (2022).

Atmospheric Hg(0) oxidation is understood to proceed via several pathways; however, the relative importance of each oxidant on the global Hg budget remains a subject of debate (Subir et al., 2011; Ariya et al., 2015). Determining the kinetics and speciation of oxidized mercury is an ongoing analytical challenge due to very low atmospheric Hg(II) concentrations (picograms per cubic meter) and limitations in instrument sensitivity (Hynes et al., 2009; Subir et al., 2011). Here, we included gas-phase oxidation of Hg(0) via Br, Cl, and OH to the Statewide Air Pollution Research Center, 1999 version (SAPRC-99; Carter, 2000) chemical mechanism in WRF-Chem, following the chemistry scheme of Shah et al. (2021). We have used the SAPRC-99 chemical mechanism for this study due to the demonstrated performance of this setup for Arctic ozone inland (see Figure 15 from Petäjä et al., 2020) and during bromine-mediated ozone depletion events (see Figures 3 and 4 from Marelle et al., 2021). Chemical concentrations (including oxidants such as OH and Br) are calculated online in the WRF-Chem model. Hg gas-phase reactions are added to a previously developed chemical mechanism (Marelle et al., 2021), which includes chlorine and bromine gas-phase chemistry, using the Kinetic PreProcessor (KPP, Sandu and Sander, 2006). Hg(0) oxidation is treated as a two-step process, first forming a Hg(I) intermediate, which can be reduced to Hg(0) (via photoreduction or thermal dissociation) before undergoing further oxidation to Hg(II) (Goodsite et al., 2004, 2012). All oxidized Hg species included in the model are listed in Table 1 (column 1).

The role of atmospheric Hg(II) reduction in the global Hg cycle is also uncertain, as rate constants are based largely on theoretical estimates yet to be verified experimentally (Goodsite et al., 2004; Balabanov et al., 2005; Dibble et al., 2012; Goodsite et al., 2012; Saiz-Lopez et al., 2018; Saiz-Lopez et al., 2019; Sitkiewicz et al., 2019; Dibble et al., 2020; Francés-Monerris et al., 2020; Khiri et al., 2020; Gómez Martn et al., 2022). Recent modelling studies have shown that atmospheric Hg(II) reduction must be considered in models to match observational estimates of total gaseous mercury lifetime against deposition (Horowitz et al., 2017; Shah et al., 2021; Zhang and Zhang, 2022). We therefore added to WRF-Chem photoreduction reactions of Hg(I) and Hg(II) species, following the implementation of Shah et al. (2021). Photochemical rates were calculated, using the FastJ photolysis scheme in WRF-Chem (Wild et al., 2000), from computationally derived quantum yields (ϕ) and absorption cross sections (σ) of Hg(I) and Hg(II) species. For Hg(I), we used ϕ = 0 for non-dissociative transitions in the wavelength range of 270–460 nm, and ϕ = 1 for transitions in the 460–800 nm range, as in Saiz-Lopez et al. (2019). Values of σ for Hg(I) species were also taken from Saiz-Lopez et al. (2019). For Hg(II), values of ϕ were taken from Francés-Monerris et al. (2020) and σ values were taken from Saiz-Lopez et al. (2018). We did not include the photoreduction of Hg(p) as the rate of this reaction remains highly uncertain.

Surface uptake of gas-phase Hg(II) onto aerosols and clouds is an important process for the formation of Hg(p) in the Arctic, particularly during spring (Steffen et al., 2013; Steffen et al., 2014). Uptake and volatilization of Hg(II) to/from surfaces is a complex process influenced by several variables including temperature, fine particle concentration, and aerosol chemical composition (Subir et al., 2012). We used a parameterized approach to represent Hg heterogeneous chemistry in WRF-Chem, using the 8-bin sectional aerosol scheme MOSAIC (Model for Simulating Aerosol Interactions and Chemistry, not related to the MOSAiC expedition) (Zaveri et al., 2008). We included the heterogeneous uptake of all Hg(II) species onto aerosols and liquid clouds. This process is treated as a kinetic process with the uptake rate, k_het, calculated for each aerosol size bin and summed to obtain a total heterogeneous uptake rate, following Equation 1 (Schwartz, 1986; Jacob, 2000):

where k_het is the total reactive uptake rate (s⁻¹), A_i and r_i are the effective mean aerosol surface area (cm²/$cmair3$⁠) and aerosol radius (cm) for aerosols in size bin i, D_g is the diffusion coefficient (cm² s⁻¹), ν is the mean molecular speed (cm s⁻¹) and α is the mass accommodation coefficient (unitless). We assumed α = 0.1 for uptake of all Hg(II) species onto aerosols and liquid clouds, following Shah et al. (2021). For liquid clouds, we calculated the uptake rate in a similar way to Equation 1, using the average cloud droplet surface area and radius for A and r, respectively (see Figure S1 for average modelled surface values). For computational efficiency, Hg(p) was not explicitly modelled in the aerosol-phase and was treated here as an additional gas-phase species.

We also included the conversion of Hg(p) on aerosols back to gas-phase Hg(II) via volatilization. This process is parameterized based on the empirical gas-particle equilibrium of Amos et al. (2012) and is considered only for fine-mode aerosols (≤2.5 μm), corresponding to the first 6 MOSAIC aerosol size bins. This empirical equilibrium, described by Amos et al. (2012), defines the ratio of Hg(II) and Hg(p) available as a function of particulate matter (PM_2.5) concentration and air temperature. Volatilized Hg(II) gas, denoted as HgX, is treated in the model as HgCl₂, which does not photolyze at tropospheric wavelengths (Saiz-Lopez et al., 2018).

The transformation between Hg(0), Hg(II), and Hg(p) has important consequences on the dry and wet deposition rates of Hg over the Arctic (Zhang et al., 2009). We included dry deposition for the added Hg species, Hg(0), Hg(II), and Hg(p), excluding Hg(I) species as they are thermally unstable and short-lived intermediates. Dry deposition of Hg(0) is calculated based on the Wesley resistance scheme (Wesely, 1989) and is implemented using four species-specific parameters: the Henry’s law constant (H*); the Henry’s law temperature correction factor (DHR); a surface reactivity factor (f₀); and the molecular diffusivity (dvj). For Hg(0), we used values of H* = 0.11 mol m⁻³ hPa⁻¹ and DHR = 4800 K from Clever et al. (1985), f₀ = 1.0 × 10⁻⁵ from Selin et al. (2008), and dvj = 0.071 cm⁻² s⁻¹ based on the inverse square root of the molar mass of Hg(0). Henry’s law constants for Hg(II) species are less well known, and so we assumed a dry deposition velocity of 0.1 cm s⁻¹ for all Hg(II) species based on estimates of reactive gaseous mercury dry deposition to ice surfaces (Lindberg et al., 2002; Skov et al., 2006; Zhang et al., 2009). Measurements of the dry deposition velocity of Hg(p) currently remain limited; however, previous results have shown that Hg(p) deposits several times slower than Hg(II) (Poissant et al., 2004; Zhang et al., 2012). We therefore assumed a dry deposition velocity of 0.01 cm s⁻¹ for Hg(p), but we note that this assumption should be revised in future work as more measurements become available.

We also added to the model wet removal of Hg(II) and Hg(p) considering both in-cloud (washout) and below-cloud (rainout) scavenging. Hg(II) is highly soluble; however, the Henry’s law constants for each individual species remain relatively unknown. We therefore assumed the wet scavenging rates of Hg(II) and Hg(p) to be equal to that of HNO₃, based on similarities in solubility and on previous Hg modelling work (Seigneur et al., 2004; Gencarelli et al., 2014). This assumption is an imperfect simplification and should be revised in future developments of the model.

The fate of deposited Hg(II) to the Arctic snowpack remains a major scientific question, with reported Hg(0) re-emission ranging between 40% and 90% of deposited Hg(II) during AMDEs (Lalonde et al., 2002; Poulain et al., 2004; Brooks et al., 2006; Kirk et al., 2006; Skov et al., 2006; Sommar et al., 2007; Ferrari et al., 2008; Dommergue et al., 2010; Durnford and Dastoor, 2011; Steffen et al., 2013; Douglas and Blum, 2019). The amount of Hg(0) re-emitted can depend on several factors including the chemical composition of snow, the amount of solar radiation, snow temperature, liquid water content of snow, and snowpack ventilation and gas transport (Steffen et al., 2002; Lalonde et al., 2003; Poulain et al., 2004; Ferrari et al., 2005; Dommergue et al., 2007; Steffen et al., 2013; Mann et al., 2015a; Mann et al., 2015b; Mann et al., 2018). Here, we added a description of Hg(0) re-emission from both land-based snow and snow on sea ice. Total reactive mercury (RM), where RM = Hg(II) + Hg(p), that has been deposited to the snow surface is tracked in the model and stored as a surface reservoir. We assumed 60% of deposited RM at the surface is photoreducible and available for re-emission under sunlit conditions, following Holmes et al. (2010) and Fisher et al. (2012). The sensitivity of this value was tested by performing a simulation with 100% of RM in snow available for re-emission. The rate of Hg(0) re-emission from land-based snow and snow on sea ice is parameterized in the model following Equations 2 and 3, respectively:

where $RHg0snow$ and $RHg0sea ice$ are the rates of Hg(0) re-emission (s⁻¹) from land-based snow and snow on sea ice, respectively; $kredsnow$ and $kredsea ice$ are the net reduction rate constants of RM (s⁻¹) from land-based snow and snow on sea ice, respectively; cos(SZA) is the cosine of the solar zenith angle (dimensionless); and f_snow and f_{sea ice} are the fractions of each grid cell (0 − 1, dimensionless) covered by snow or sea ice, respectively. For land-based snow, we used $kredsnow$ = 2.5 × 10⁻⁵ s⁻¹ (0.09 h⁻¹), following Poulain et al. (2004), and for snow on sea ice we assumed a lower value of $kredsea ice$ = 2.8 × 10⁻⁷ s⁻¹ (0.001 h⁻¹) based on a mid-range value from observational estimates (Durnford and Dastoor, 2011), following Fisher et al. (2012). Current estimates of k_red in Arctic snow samples have a large variability, with reported values ranging from 7 × 10⁻⁶ to 0.6 h⁻¹ (Poulain et al., 2004; Dommergue et al., 2007; Mann et al., 2015b; Mann et al., 2018). Part of this variability can be explained by differences in measurement techniques, with some values reporting gross Hg(II) photoreduction rate constants, thereby neglecting the effects of in-snow Hg(0) oxidation (e.g., Mann et al., 2015b; Mann et al., 2018). The values of k_red used here are therefore tentative estimates to be examined further in future studies. A temperature dependence was applied for Hg(0) re-emissions, where re-emission is only active for snow-covered grid cells with a skin temperature below 0°C. We also applied a fractional sea ice cutoff of 75% for Hg(0) re-emission, in which grid cells below this sea ice threshold do not re-emit Hg(0), following similar descriptions of Br₂ emissions in Marelle et al. (2021).

We set up the model for the dates between March 1 and April 14, 2020, during leg 3 of the MOSAiC expedition. In mid- to late April 2020, warm air intrusions were observed reaching the Polarstern, transporting pollution from southern latitudes to the central Arctic (Dada et al., 2022). We therefore excluded this period from our model simulation as these events require a detailed evaluation of pollution transport which is beyond the scope of this study. The first 2 weeks of model output are considered spin-up and are not included in the analysis (see Section 4.2). A horizontal resolution of 100 × 100 km is used to encompass the entire Arctic (model domain shown in Figure 2) with a vertical resolution of 72 levels up to a pressure of 50 hPa. Details about the model setup are given below and are selected based on extensive testing for the Arctic from previous studies (Marelle et al., 2017; Marelle et al., 2021). A summary of the specific model chemistry and physics options selected are given in Table 2.

Initial and boundary conditions for meteorology were obtained from the National Centers for Environmental Prediction Final Analysis (NCEP FNL; National Centers for Environmental Prediction, 2000) and we applied spectral nudging within and above the boundary layer. The boundary layer scheme and land surface model used are the Mellor–Yamada Nakanishi Niino Level 2.5 Scheme (MYNN; Nakanishi and Niino, 2009) and the Noah Land Surface Model (Noah-LSM; Tewari et al., 2004), respectively. These options were selected as they were found to give the best representation of boundary layer dynamics in the Arctic from a number of different WRF dynamics configurations (Marelle et al., 2021). Furthermore, we included an additional sea ice thickness variable from the TOPAZ4b model to better simulate surface temperatures over sea ice (EU Copernicus Marine Service Information, 2022).

Initial and boundary conditions for chemical concentrations were set using data from the global CAM-Chem model (Buchholz et al., 2019; Emmons et al., 2020). We also included an initial and boundary concentration of atmospheric Hg(0), obtained from model output simulated by an ensemble of four chemical transport models (see Section 3.4). More details about the model ensemble and its performance in the Arctic can be found in Dastoor et al. (2022a). Concentrations of RM were initialized as zero and were allowed to reach a natural equilibrium during the model spin-up period.

We considered bromine emissions from snow/ice and aerosols using both surface snow and blowing snow sources of reactive bromine emissions and recycling, as described in Marelle et al. (2021). These are interactive parameterizations within the model that calculate Br₂ emissions online. Specifically, the surface snow mechanism calculates the emission of Br₂ following the deposition of ozone to the snowpack over sea ice (Toyota et al., 2011). This process is accelerated under sunlit conditions compared to dark conditions. The blowing snow mechanism calculates the release of Br₂ and sea salt aerosols from the sublimation of lofted snow under windy conditions (wind speeds >7 m s⁻¹; Yang et al., 2008; Yang et al., 2019). Finally, heterogeneous recycling of bromine on all frozen surfaces is also considered in this work, following Marelle et al. (2021). These interactive bromine emissions were evaluated in Marelle et al. (2021) by comparing with Arctic ozone observations, including in the central Arctic from ice-tethered buoys (O-Buoys; Knepp et al., 2010; Halfacre et al., 2014), showing a good representation of springtime ozone depletion.

Biogenic emissions are calculated online using the Model of Emissions of Gases and Aerosols from Nature (MEGAN; Guenther et al., 2012). Fire emissions were obtained from the Fire INventory from NCAR version 2.5 (FINNv2.5; Wiedinmyer et al., 2011; Wiedinmyer et al., 2023) and anthropogenic emissions are from the global ECLIPSEv6b inventory (Evaluating of the Climate and Air Quality Impacts of Short-Lived Pollutants version 6b; Klimont et al., 2017). ECLIPSEv6b includes revised international shipping emissions compared to previous versions of the inventory, relevant for capturing local Arctic emission sources. We further added anthropogenic emissions of Hg(0), Hg(II), and Hg(p) from the global anthropogenic mercury emissions inventory for 2015 (Steenhuisen and Wilson, 2022). This inventory was prepared as part of the 2018 AMAP/UNEP Global Mercury Assessment (GMA; AMAP/UN Environment, 2019) and groups emissions into four distinct sectors: fuel combustion, industrial sectors, waste from intentional use, and artisanal and small-scale gold mining (Steenhuisen and Wilson, 2019, 2022). To include anthropogenic emissions of gas-phase Hg(II) in the model, we assumed that these emissions are evenly distributed between each gas-phase Hg(II) species in the model.

The MOSAiC expedition took place on board the Polarstern icebreaker, between October 2019 and September 2020, and is the largest scientific exploration of the Arctic to date. A comprehensive suite of measurements were made during drift through the Arctic Ocean to better understand the links between the atmosphere, sea ice, and ocean (Nicolaus et al., 2022; Rabe et al., 2022; Shupe et al., 2022). This extensive dataset provides a unique insight into the chemical behaviour of the atmosphere in the central Arctic. More details about the MOSAiC campaign can be found in Shupe et al. (2022). For the purpose of our modelling study, we used a subset of meteorological and chemical observations made during MOSAiC which are introduced below.

Continuous surface meteorological measurements were made onboard the Polarstern at various locations and heights. Measurements of air temperature and relative humidity were made at a height of 29 m above sea level using a Vaisala HMP155 probe. Wind speed and wind direction were measured at 39 m above sea level using an ultrasonic anemometer. The meteorological measurement data used here to evaluate the WRF-Chem model are available at Schmithüsen (2021).

Radiosondes (model Vaisala RS41-SGP) measuring temperature, humidity, and winds were routinely launched every 6 hours throughout the whole expedition (Maturilli et al., 2021). During synoptic events of special interest, such as major storms, the launch frequency was enhanced up to three-hourly. The launches were performed from the helicopter deck of the Polarstern, approximately 12 m above the sea level, so the data do not capture the lowermost part of the boundary layer. Furthermore, polluted data may exist in the lowermost approximate 100 m, in cases when the sonde flew through the ship’s exhaust fan or was otherwise influenced by the ship’s presence.

As described in Angot et al. (2022b), Hg(0) measurements were performed in the University of Colorado sea-container laboratory using a Tekran 2537B analyzer. Only Hg(0) was collected and analyzed (as opposed to total gaseous mercury) as cation-exchange membranes were used to remove potential divalent Hg species. All instruments located within this container were automatically backflushed with zero-air when wind direction was more than ±130 degrees from the ship bow to prevent any contamination of the sampling line by the ship exhaust. This backflushing explains the presence of gaps in the Hg(0) time series. The Hg(0) dataset can be accessed at Angot et al. (2022a).

Ozone ambient-air mole fractions were monitored in three different sea-container laboratories using commercial instruments. The three individual ozone time series were cross-evaluated and used to generate the hourly averaged merged dataset used in this study. This merged dataset limits gaps in the ozone time series (as opposed to the Hg(0) time series; see above). More information can be found in Angot et al. (2022b). The O₃ dataset can be accessed at Angot et al. (2022c).

The Multi-Axis Differential Optical Absorption Spectroscopy (MAX-DOAS) technique (Lohberger et al., 2004; Plane and Saiz-Lopez, 2006; Platt and Stutz, 2008) was used to make observations of BrO during the MOSAiC expedition. The instrument (Prados-Roman et al., 2015) has an external telescope used to collect scattered sunlight connected to an indoor spectrometer using an optical fibre. The outdoor unit consisted of a telescopic lens (focal length of 200 mm, diameter 50.8 mm), which focused light onto an optical fibre. The optical fibre was connected to the indoor unit consisting of a spectrometer (Princeton Instruments SP500i) and a Charge-Coupled Device detector (camera CCD Princeton Instruments Pixis 400B). The instrument was placed at a height of 15 m from the sea surface and an inclinometer was used to correct the measured elevation angles using the actual pitch and roll of the instrument during the expedition. The solar spectra were then analysed using the QDOAS software (Fayt et al., 2011) to retrieve the absorption due to the oxygen dimer (O₄) and BrO at different viewing elevation angles, using the zenith spectra as a reference (see Figure S2). The resulting Differential Slant Column Densities (DSCDs) are the difference of the slant column densities (SCDs) in the viewing direction and the SCD in the zenith direction. The DOAS retrieval settings are given in Table S1, and an example of the DOAS fit can be found in Benavent et al. (2022). The AC-2 radiative transfer model (Benavent, 2020) was then used to derive mixing ratios of BrO by estimating the path length from the O₄ DSCDs. This is a two-step process where the O₄ DSCDs are used to estimate the light path length and subsequently the mixing ratio of BrO is estimated. Considering that in most cases BrO was above the detection limit only at the lowermost angles (<3 degrees), an average mixing ratio across the first kilometer is computed due to lack of information above it. Photographs of sky conditions were used to filter the data for blowing snow and broken cloud cover to avoid multiple scattering effects, which can lead to the incorrect conversion of DSCDs into mixing ratios.

Measurements of total gaseous mercury at Villum Research Station (Station Nord, Greenland), Zeppelin Observatory (Svalbard), and the Dr. Neil Trivett Global Atmosphere Watch Observatory (Alert, Nunavut, Canada) were performed using Tekran 2537 instruments and are part of ongoing long-term monitoring efforts (e.g., MacSween et al., 2022). Teflon inlet filters were used to remove gaseous oxidized mercury during sampling, but it is possible that small amounts of Hg(II) are measured by the analyzers. This contribution is assumed to be small, relative to the overall Hg(0) signal, and therefore these measurements are henceforth referred to as Hg(0). Hg(II) and Hg(p) measurements at Alert were made using a Tekran 1130/1135 speciation unit, sampled at a 2-hour time resolution and reported on a 3-hour time interval (Steffen et al., 2014). More information on the instrumental setup and quality control procedures can be found elsewhere (Berg et al., 2013; Steffen et al., 2014; Angot et al., 2016; Skov et al., 2020).

Surface ozone data from Villum (Greenland) and Zeppelin (Svalbard) were retrieved from the EBAS database (http://ebas-data.nilu.no/default.aspx). All measurements were performed using commercial UV absorption instruments (detection limit of 1 ppb). Surface ozone data from Utqia͘gvik (Alaska), Summit (Greenland), and the Polar Environment Atmospheric Research Laboratory (Eureka, Nunavut, Canada) were provided by the NOAA Global Monitoring Laboratory (https://www.esrl.noaa.gov/gmd/ozwv/surfoz/data.html). These measurements are also performed using a commercial UV absorption instrument. More information can be found in Platt et al. (2022).

Satellite BrO data of the high-resolution TROPOspheric Monitoring Instrument (TROPOMI) on the Sentinel-5 Precursor (S5P) satellite were used to retrieve total BrO columns by applying the method from Seo et al. (2019). To obtain the tropospheric BrO Slant Column Density (SCD), the stratospheric part was removed from the total column using the stratospheric correction method described in Theys et al. (2011). It requires stratospheric ozone and NO₂ columns from TROPOMI and the tropopause height, which were taken from the NCEP Reanalysis 1 product (Kalnay et al., 1996). The tropospheric BrO Vertical Column Densities (VCDs) were obtained from the tropospheric SCDs using an air mass factor suitable for a surface BrO layer of 400 m thickness over a bright surface. As a result, boundary layer BrO over dark surfaces such as open oceans may be underestimated.

Simulations for the MOSAiC year (Oct 2019 to Sept 2020) were performed with the multi-model ensemble (GLEMOS, GEOS-Chem, GEM-MACH-Hg, and DEHM) used for the recent Global Mercury Assessment (Arctic Monitoring and Assessment Programme/UN Environment, 2019) and Arctic Monitoring and Assessment Programme reports (AMAP, 2021). A description of the four models can be found in Dastoor et al. (2022a), along with a full evaluation of their performances in the Arctic. A subset of this simulation (March–April 2020) was used to provide initial boundary conditions of Hg(0) in the model simulation presented here.

We evaluated the hourly model outputs against meteorological and chemical measurements made onboard the Polarstern during MOSAiC. We used several metrics to provide a comparison between the model and observations, including the Pearson correlation coefficient (r), root-mean-square error (RMSE), and mean bias error (MBE). Performance metrics between the model and the MOSAiC observations are given in Table 3.

Polar boundary layer stability is crucial in modulating atmospheric chemical composition near the surface via impacts on vertical mixing, surface emissions, and chemistry. Accurately modelling boundary layer structure is therefore necessary to investigate surface chemistry and air-snow exchange processes. We evaluated the meteorological configuration of the model by comparing the simulated boundary layer meteorology with observations obtained onboard the Polarstern. Figure 3 shows the comparison between the simulated and observed meteorology for several variables including air temperature, relative humidity, wind speed, and wind direction. Simulated variables were extracted at the nearest model grid cell to the location of the ship. Modelled air temperatures are largely in good agreement with the observations (r = 0.85; Figure 3a), despite a mean positive bias (MBE = 1.95 K) slightly overestimating both warm and cold periods. The model also shows good performance in simulating surface wind speeds (r = 0.97) and wind directions, with the exception of relative humidity (r = 0.31), which is biased low (MBE = −7.69%).

To evaluate the vertical structure of the boundary layer, we compared the average simulated and observed temperature profiles above the Polarstern (Figure 4). Radiosondes were launched from the ship up to four times a day at regular intervals: 05:00, 11:00, 17:00, and 23:00 UTC (to the nearest hour). In total, 119 radiosondes were launched during the simulated period. Figure 4 shows the comparison between the mean observed and simulated temperature profiles. For this comparison, we have interpolated the model output to the hour/location of the sonde releases and averaged all profiles available for each daily release time (05:00, 11:00, 18:00, 23:00 UTC); we show the lowest two kilometers of the atmosphere to focus on the near-surface atmosphere representation in the model. The purpose of this comparison is to provide an indication of overall model performance in simulating the vertical structure of the polar boundary layer. At each time interval, the model shows good agreement with the mean observed vertical temperature profiles, within the standard deviation range of the observations. In the lowest 250 m, we find a positive bias of 1–2°C in the model, consistent with the warm bias in Figure 3a. This warm bias can be associated with the presence of near-surface temperature inversions, common in the Arctic. The difficulty in capturing these inversions is not unique to our model, and is a challenge facing many models and reanalyses (Graham et al., 2019). Similar comparisons between modelled and observed vertical profiles of winds and relative humidity (Figures S3–S5) also show good model performance in simulating the vertical structure of the boundary layer, despite some differences in the lowest 250 m. Overall, we show that the model performs well in simulating the observed springtime meteorology in the central Arctic, necessary for modelling boundary layer chemistry.

To evaluate the model, we compared modelled surface Hg(0), O₃, and BrO with observations made onboard the Polarstern (Figure 5). Performance metrics are shown in Table 3. We also note here that a 2-week model spin-up time was selected based on the time needed to stabilize the chemical state of Hg(0) and O₃ in the central Arctic (grey shaded area in Figure 5). Observations of Hg(0) show extended periods of depletion (Figure 5a), defined as Hg(0) < 1 ng m⁻³ (Angot et al., 2016) and concentrations were frequently recorded below 0.5 ng m⁻³. Similarly, intense depletion of ozone was observed (Figure 5b), where O₃ was regularly below 10 ppb and depletion lasted multiple days, particularly in late March/early April. This level of Hg(0) and O₃ depletion is indicative of enhanced surface chemistry and emissions of bromine from sea ice, capable of sustaining the depletion over an extended period of time. There is also a strong correlation between the observations of Hg(0) and O₃ (r = 0.91), suggesting depletion of both species is driven by bromine, as reported in a previous Arctic measurement study (Wang et al., 2019). Hg(0) concentrations during periods of missing measurement data (e.g., March 18–23 and March 31 to April 3) are therefore likely to be consistent with O₃ concentrations, which exhibit regular depletion during these periods. Our model predicts substantial depletion of both Hg(0) and O₃, in relatively good agreement with the observations (r = 0.81 and r = 0.82 for Hg(0) and O₃, respectively). In particular, the model is able to capture the prolonged periods of Hg(0) and O₃ depletion in late March/early April. However, simulated Hg(0) has a small negative bias (MBE = −0.09 ng m⁻³), potentially due to high bromine-initiated oxidation, and the model also underestimates certain periods of ozone recovery (e.g., March 21–23). Comparing the simulated time series of Hg(0) and O₃, we find a very strong correlation (r = 0.99), indicating oxidation by Br is the major process in the model.

To further assess the extent of bromine-initiated oxidation of Hg(0) and O₃, we compared the average modelled BrO in the lowest 1 km of the atmosphere to MAX-DOAS observations (Figure 5c). The TROPOMI tropospheric column BrO is shown along the MOSAiC ship track as a reference for comparison (right axis in Figure 5c). We highlight here that there is some uncertainty in directly comparing the lowest 1 km of model data with the MAX-DOAS measurements; however, this uncertainty is difficult to quantify. Our model predicts high BrO concentrations above the shiptrack of up to 40 ppt. The observed concentrations are up to a factor of 5 times smaller than the model predicted estimates (r = 0.42, MBE = 9.78 ppt), with many periods of BrO below the detection limit and uncertainty typically below 1 pptv. The average concentrations in the lowest 1 km of our WRF-Chem simulation are, however, consistent with timing and relative abundance of BrO from the total tropospheric column retrieved from TROPOMI (Figure 5c). Putting these measurement data in the context of previous Arctic BrO observations, they are much lower than peak levels of around 50 pptv recorded from satellite and ground-based observations (Wagner et al., 2001; Hönninger and Platt, 2002; Peterson et al., 2017; Simpson et al., 2017). These observations of BrO during the MOSAiC spring season were made mainly during periods when ozone was near complete depletion. That reactive bromine can be present even when there are low BrO concentrations is well known, due to the fact that BrO formation cannot occur when O₃ is depleted.

There are several possible causes for the difference between modelled and observed BrO. The main loss of BrO is via reaction with hydroperoxyl radicals (HO₂), forming HOBr. Inaccurate HO₂ concentrations from VOC oxidation may result in insufficient HOBr production, which could maintain high concentrations of BrO. In addition, our description of bromine activation assumes an infinite bromide reservoir from snow on both first-year sea ice and multi-year sea ice (Herrmann et al., 2022). This assumption may overestimate bromine activation from multi-year ice regions, where measurements of bromide in snow on Arctic multi-year ice have recorded lower concentrations than on first-year ice (Krnavek et al., 2012; Peterson et al., 2019). Uncertainties related to the BrO retrievals that allow for calculation of the BrO concentrations from the MAX-DOAS observations may also contribute to these differences. Finally, heterogeneous recycling on aerosols may sustain reactive bromine concentrations above the surface, independently of snowpack activation (Peterson et al., 2017). An overestimation in modelled aerosol concentrations, or a high reactive uptake probability of bromine, may therefore contribute to an enhancement of bromine recycling and BrO concentrations at altitudes above the surface.

In summary, the model is able to simulate Hg(0) and O₃ depletion in the central Arctic, in good agreement with the observations. There is strong correlation between the Hg(0) and O₃ time series (observed and modelled), indicating oxidation of both species is driven by bromine. However, modelled BrO quantities are overestimated compared to the observations, which may positively bias oxidation of Hg(0) and O₃ by bromine. A recent study, based on reactive bromine and iodine measurements, has attributed a large contribution of ozone destruction to iodine chemistry, on a level comparable with bromine (Benavent et al., 2022). These measurements were made during MOSAiC, at the same time and location modelled in this study, and are the first measurements of iodine monoxide (IO) in the central Arctic. In this work, we did not include iodine chemistry in our model, and our halogen descriptions are limited to bromine and chlorine cycling. Future investigations could aim to extend our model chemical mechanism to include descriptions of this chemistry and to explore its impact on ozone depletion. For this work focused on the mercury cycle, the impact of iodine chemistry is expected to be minimal on mercury oxidation and deposition.

In addition to the model comparison with MOSAiC observations, we further evaluated the model with measurement data from other Arctic locations. Here, we compare the simulated surface Hg and O₃ concentrations to observations at Arctic stations, as well as modelled tropospheric BrO vertical column density with satellite-derived measurements. An additional comparison of meteorology at 3 Arctic stations is also provided in Figures S6–S8.

Figure 6 compares modelled surface Hg(0) to hourly averaged observations at Villum (Greenland) and Zeppelin (Svalbard), and 3-hour averaged observations at Alert (Canada). Overall, we find relatively good agreement between the model and observations in capturing the timing of AMDEs at each measurement site. There is a general negative bias in the simulated Hg(0) time series compared to the observations (Table 4), suggesting either too much net Hg(0) oxidation (via high bromine concentrations) or low/missing sources (including re-emission from snow, open ocean sources, and transport from mid-latitudes).

At Villum (Figure 6a), the model is able to simulate key features of the observed Hg(0) time series (r = 0.50), with some discrepancies discussed below. During the simulated period, six AMDEs were observed (March 27, March 30–31, April 2–4, April 9–10, April 10–13, and April 14). The model is able to simulate most features (timing and intensity) of five of the six events. On April 2, a depletion event lasting approximately 2 days was observed, where Hg(0) concentrations fell below 0.5 ng m⁻³. The intensity and duration of this event is captured by the model, although recovery to background levels in the model lags by approximately 1 day, possibly suggesting low re-emission from snow in the model. The model incorrectly predicts a depletion event between March 24 and March 29 which was not observed. These differences can partly be explained by the coarse horizontal resolution of the model (100 km). Due to the proximity of Villum to sea ice, this grid cell may partly factor in emissions of bromine from sea ice. We therefore also show in Figure 6 the minimum and maximum Hg(0) values within the two nearest neighbour grid cells (shaded regions). The large variability of Hg(0) concentration within neighbouring grid cells highlights the marked differences in emissions and chemistry of Hg(0) over land-based snow and sea ice.

At Alert (Figure 6b), Hg(0) is regularly close to background levels (approximately 1.4 ng m⁻³), with some intermittent periods of Hg(0) depletion. Simulated Hg(0) shows reasonable agreement with the observations (r = 0.54) with a small negative bias (MBE = −0.05 ng m⁻³). Speciated Hg measurements at Alert were also available during this period and are evaluated in more detail in Section 5.2. Observations at Zeppelin (Figure 6c) recorded several depletion events lasting 1–3 days. There is a broad underestimation of modelled Hg(0) compared to the observations (r = 0.53), which could be explained in part by low/missing land-based sources (including re-emission from snow). Limitations in model resolution may also contribute to this negative bias by poorly representing the local mountain meteorology at Zeppelin. Chemistry, emissions, and transport at this site can all be impacted as a consequence. Despite the mean negative bias (MBE = −0.49 ng m⁻³), the timing of depletion and replenishment of Hg(0) are often captured by the model at Zeppelin.

Overall, the model is generally able to reproduce the behaviour of Hg(0) at these sites, with the timing of depletion events often simulated correctly by the model, but underestimating Hg(0) concentrations. The negative bias in modelled Hg(0) concentration is likely a factor of the coarse horizontal resolution of the model, resulting in high bromine oxidation and low Hg(0) re-emission at these sites. Higher resolution model runs are therefore desirable to discern the contributions of bromine oxidation and Hg(0) re-emission from land-based snow and from sea ice at coastal Arctic stations.

Figure 7 compares modelled and hourly averaged observations of surface O₃ at six Arctic stations: Villum (Greenland), Alert (Canada), Zeppelin (Svalbard), Utqia͘gvik (Alaska), Summit (Greenland), and Eureka (Canada). Model performance metrics are listed in Table 4. At Villum (Figure 7a), observations of O₃ are correlated with the Hg(0) time series (r = 0.93) shown in Figure 6a. Modelled O₃ is in reasonable agreement with the observations (r = 0.54), with similar features and discrepancies as previously shown in the modelled Hg(0) time series. In particular, an ODE during April 2–4 was observed coinciding with depletion of Hg(0), which is reproduced by the model. However, the model also incorrectly predicts depletion of O₃ between March 23 and March 31, indicating an overestimation of bromine. Modelled O₃ at Alert (Figure 7b) is relatively well captured by the model (r = 0.64). The observed time series of O₃ and Hg(0) also show very similar features, with a depletion event recorded for both species during March 23–25. At Zeppelin (Figure 7c), multiple ODEs were observed and O₃ concentrations correlate strongly with the Hg(0) measurements (r = 0.88). The timing and intensity of these ODEs are mostly captured by the model (r = 0.56). However, there is an overall mean negative bias in the model (MBE = −7.61 ppb), with a marked underestimation between March 23 and March 26 (as in the Hg(0) comparison). This bias is again suggestive of an overestimation in bromine activation during this period. Comparing the modelled O₃ and Hg(0) time series at these three stations, we find a very strong correlation at each site (Villum r = 0.73, Alert r = 0.88, and Zeppelin r = 0.99), indicating that bromine is also the major Hg(0) oxidant in the coastal Arctic.

Ozone measurements at Utqia͘gvik show several ODEs which the model is only partially able to capture (Figure 7d). In general, the model overpredicts ozone depletion indicating highly active bromine chemistry (r = 0.52, MBE = −11.3 ppb). At Summit (Figure 7e), no depletion events were observed or simulated by the model due to the location of the station (r = 0.64). Summit is situated at a high altitude (3216 m above sea level) in central Greenland and is therefore largely unaffected by air masses originating over sea ice. There is an overall positive bias in modelled O₃ concentrations (MBE = 6.67 ppb), indicating that the difference may be due to local meteorology (e.g., boundary layer stability) rather than differences in chemistry. At Eureka station (Figure 7f), multiple ODEs were observed typically lasting 1–2 days. For most of the simulated period (March 18–April 14), the model performs very well to capture the timing and intensity of ODEs observed at this site (r = 0.52). However, ozone concentrations in the first week of the simulation are overpredicted and the depletion is completely missed by the model.

Figure 8 shows the average surface O₃ concentration during the simulation period. On the whole, we find good representation of ozone in the model at several coastal and non-coastal Arctic stations. Most ODEs (including the timing, intensity, and duration) are captured by the model, with the exception of some events which are missed (e.g., at Eureka) or overpredicted by the model (e.g., at Villum). For coastal sites, there is a general negative bias in O₃ suggesting an overestimation in bromine-initiated oxidation. Increased horizontal resolution could again potentially address some of the disparities between model and observations at all sites.

We compare in Figure 9 the modelled and observed mean tropospheric BrO vertical column densities for the simulation period. Satellite observations from TROPOMI on Sentinel-5p are corrected to exclude the stratospheric contribution (see Section 3.3). The first striking difference between TROPOMI (Figure 9a) and the modelled BrO (Figure 9b) is the different total columns outside the central Arctic over primarily open ocean regions (south-east of Greenland and between Alaska and Russia). This finding can be partly explained by the fact that over open oceans, the sensitivity of satellite observations to boundary layer BrO is largely reduced due to low albedo and cloud effects (Seo et al., 2019), and therefore satellite columns may be too low. The difference also points to the fact that the BrO lifetime in our model may be too long and allows for transport outside the sea ice-covered source region for initial bromine activation or the fact that we have too much bromine recycling on surfaces outside sea ice-covered regions. Peterson et al. (2017) reported observations of BrO in the Arctic boundary layer sustained by heterogeneous recycling on aerosols, which is one way that reactive bromine can be sustained in air away from the surface. Reactive bromine concentrations are sustained on aerosols and transported to higher altitudes in the free troposphere where they can be transported longer distances away from the original emissions sources over sea ice-covered regions. However, the rate of recycling and re-emissions on aerosols remains uncertain. To explore the effect of recycling on aerosols away from the surface and its impact on sustaining BrO activation away from sea ice, we completed a sensitivity run with aerosol bromine recycling turned off. The result of this run (Figure 9c) is in better agreement with BrO VCD observed, with an underestimation in BrO abundances. This result is clearly not a realistic representation of bromine aerosol chemistry, but provides valuable information regarding the role of bromine activation and transport via recycling aerosols in the Arctic. It shows the need for future work on both evaluating modelled aerosol concentrations and/or inaccuracies in the treatment of heterogeneous bromine chemistry.

There are also differences in the central Arctic between TROPOMI (Figure 9a) and the modelled BrO (Figure 9b) VCDs, over the sea ice-covered regions, for example at the North Pole (90°N). Here and over other sea ice-covered regions the model underpredicts the total BrO column compared to the satellite VCD. There are several possible reasons for underprediction, including the fact that the model underpredicts total BrO due to uncertainties in emissions, vertical transport, and recycling on aerosols (discussed in the previous paragraph). A key issue that may control BrO abundance near the surface in the model is the ability to replenish ozone-rich air from above the surface down to the key bromine emission sources near the surface, which are snow on sea ice and sea salt aerosols. When ozone is depleted, BrO cannot form even if activated bromine radicals are present. The underprediction of central Arctic BrO VCDs by the model can also be due to the inaccurate representation of vertical transport of activated bromine away from the surface into the free troposphere, where large quantities of ozone are available for reaction to form BrO. The simulated boundary layer (Figure S9) was often below 700 m, which indicates the volume in which emitted bromine is in direct contact with the surface and can be recycled on surface snow and aerosols. BrO formed near the surface was found to be mixed above the boundary layer height (Figure S10), with ozone depletion extending up to 2 km above the surface (Figure S11). While aloft, recycling of bromine is limited to aerosol surfaces as shown in Peterson et al. (2017), resulting in less frequent events of elevated BrO aloft above the boundary layer height.

In Figure 10, we have plotted the modelled chemical speciation of surface Hg(0), Hg(II), and Hg(p) averaged over the simulation period. Observational averages from Arctic stations (where possible) and the MOSAiC shiptrack are also plotted for the same period. Modelled Hg(0) shows depletion in the central Arctic, due to bromine emissions from sea ice, with higher Hg(0) concentrations over land (Figure 10a). The simulated surface average is in good agreement with the available observations and predicts a latitudinal gradient with increasing Hg(0) concentrations at lower latitudes. Simulated surface gaseous Hg(II) concentrations are in the range of 0–100 pg m⁻³, with the highest concentrations around the coasts and minima over the central Arctic (Figure 10b). The five major contributing species in the model are HgX (gaseous Hg(II) volatilized from aerosols), HgOHOH, HgBrOH, HgBrBrO, and HgBr₂, indicating high bromine oxidation (Table 1). Measurements at Alert report an average Hg(II) concentration of approximately 10 pg m⁻³ during this period, lower than the model average of 48 pg m⁻³. Modelled Hg(p) concentrations (Figure 10c) are close to an order of magnitude greater than Hg(II). During the simulated period, the mean Hg(p) concentration observed at Alert was 125 pg m⁻³, with the model average in close agreement (158 pg m⁻³). The tropospheric budget of Hg (Table 1) shows that approximately 63% of all oxidized mercury in the model is present as Hg(p), indicating high aerosol processing. We note here that we do not include the photoreduction of Hg(p) on organic aerosols to Hg(0), as in previous Hg modelling studies (e.g., Shah et al., 2021; Zhang and Zhang, 2022), as this photolysis rate is still poorly constrained. By including this process we could expect to see a decrease in the modelled average Hg(p) concentration, and an increase in mean surface Hg(0) concentrations. Overall, the mean modelled Hg(0), Hg(II), and Hg(p) concentrations show reasonable agreement with Arctic observations, within the measurement uncertainty.

We next consider these results in the context of previous Hg speciation measurements in the Arctic during spring. Observations have consistently shown greater Hg(p) concentrations than gaseous Hg(II) in early spring, over both tundra and sea ice (Cobbett et al., 2007; Manca et al., 2013; Steffen et al., 2013; Steffen et al., 2014; MacSween et al., 2022). Long-term observations of Hg at Alert (2002–2011) reported average Hg(II) concentrations of 11 pg m⁻³ and 34 pg m⁻³ and Hg(p) concentrations of 137 pg m⁻³ and 150 pg m⁻³ during March and April, respectively (Steffen et al., 2014). These amounts are within the range of our model results, demonstrating a higher fraction of modelled Hg(p) during early spring. From the same long-term dataset, Hg(p) concentrations reached an annual maximum in April, while Hg(II) concentrations peaked during May. The reason for the springtime peak in Hg(p) has been investigated and associated with several variables. Low temperatures have been shown to increase the transition of gaseous Hg(II) to Hg(p), as well as a high fraction of organic and chloride containing aerosols (Amos et al., 2012; Steffen et al., 2014; Xu et al., 2020). In addition, the transport of air pollution from mid-latitudes to the Arctic during winter and spring (known as Arctic haze) brings high aerosol concentrations, contributing to the predominance of Hg(p) in spring. Our results are in agreement with these findings, with maximum Hg(p) simulated in the central Arctic, indicating the influence of high aerosol concentrations and low temperatures on Hg(II)/Hg(p) partitioning. Furthermore, Steffen et al. (2013) reported mean concentrations of 30 pg m⁻³ and 393 pg m⁻³ for Hg(II) and Hg(p), respectively, over Arctic sea ice in spring 2009. Observations of Hg(0) made concurrently over both snow-covered tundra and sea ice showed that Hg(0) was often significantly higher over tundra than over sea ice, indicating higher re-emission (Steffen et al., 2013). This finding is consistent with the observations in 2020 (Figure 10a) as well as the model prediction, where Hg(0) concentrations over continental snow are higher than over sea ice. Hg(0) re-emission over land-based snow and sea ice is discussed further in Section 5.4.

Figure 11 shows the time series of simulated and 3-hour averaged observations of gaseous Hg(II), Hg(p), and total Hg at Alert (Canada). The model overestimates Hg(II) at this site (Figure 11a), and for Hg(p), we find relatively good agreement with the observations during March, but overestimated Hg(p) concentrations in April (Figure 11b). Several factors may contribute to this overestimation. Firstly, uncertainties in the measurements may partially explain some of this difference (Gustin et al., 2015). There is growing evidence that Hg speciation measurements may be biased low by a factor of 1.5–12 (e.g., Gustin et al., 2013; Huang et al., 2013; Gustin et al., 2015; Osterwalder et al., 2021) due to analytical challenges. Whether measurements in high latitudes are biased low or not is, however, still a matter of discussion as relative humidity and ozone levels (both shown to influence the collection efficiency of the denuders) are typically low. Model inaccuracies in Hg gas-particle partitioning and other relevant processes are also likely to contribute. As shown in Figure 6b, modelled Hg(0) concentrations at Alert are in good agreement with the observations. Additionally, the total Hg concentration is well reproduced by the model (Figure 11c), suggesting that high Hg(II) and Hg(p) concentrations are unlikely to be caused solely by an overestimation in Hg(0) oxidation or re-emissions. Consequently, the overestimation of Hg(II) and Hg(p) indicates a potential underestimation in the loss of Hg(II) and Hg(p) via deposition. The relative amounts of atmospheric Hg(II) and Hg(p) have been shown to directly impact snow Hg concentrations in the Arctic (Steffen et al., 2014). Mercury deposition to snow and sea ice is discussed further in Section 5.3. Finally, we once again note the potential role of coarse model resolution on differences between the model and observations, with large variation shown for Hg(II) and Hg(p) in neighbouring grid cells (shaded regions in Figure 6).

Overall, the model is capable of simulating the average springtime speciation of Hg(0), Hg(II), and Hg(p) over snow and sea ice compared to observations. Hg(II)/Hg(p) partitioning at Alert was not fully captured by the model and more work is needed to refine/test the relevant processes including Hg(II) and Hg(p) deposition, sensitivity of Hg(II) uptake to aerosols, photoreduction of Hg(p), and uncertainties in the gas-particle partitioning equilibrium. Recent analysis of trends in Arctic Hg speciation over the past 20 years suggests that the composition and timing of AMDEs are changing (MacSween et al., 2022). Therefore, more observations of speciated Hg at high latitudes are also necessary to better evaluate models and to monitor shifts in springtime Hg chemistry.

Modelled Hg deposition fluxes (dry and wet) are presented in Figure 12. We find that dry deposition (Figure 12a) is the main Hg deposition process over land and sea ice. The majority of dry-deposited Hg in the model is from Hg(II) and Hg(p), accounting for 88% of total Hg dry deposition, with only 12% from Hg(0). Over sea ice, the model predicts an even greater fraction of dry deposition from RM (approximately 98%). The percentage contribution of RM to dry deposition is consistent with the higher fraction of RM during spring, compared with other seasons. For wet removal (Figure 12b), the model predicts maximum wet deposition over the open ocean, consistent with previous Hg modelling studies (Holmes et al., 2010; Horowitz et al., 2017; Travnikov et al., 2017; Shah et al., 2021; Zhang and Zhang, 2022). We also find minimal wet deposition of Hg in the central Arctic and over land. This finding is unsurprising as the Arctic typically exhibits low precipitation rates (snow and rain) during spring, and consequently low wet deposition. To contextualize these results, we next compare our modelled deposition fluxes with previous observational and modelling results in the Arctic.

Annually averaged observations of Hg wet-deposition fluxes in the Arctic and sub-Arctic have been reported in the range of 0–5 g m⁻² yr⁻¹ (Sanei et al., 2010; Sprovieri et al., 2017; Pearson et al., 2019). In Alaska, mean annual Hg wet-deposition fluxes of 2–5 g m⁻² yr⁻¹ were recorded at 5 different locations, over several years (2008–2015; Pearson et al., 2019). Sprovieri et al. (2017) also reported multi-year (2012–2015) mean fluxes of Hg wet deposition of 0.8–1.7 g m⁻² yr⁻¹ at Ny-Ålesund, Svalbard. Furthermore, measurements from two Canadian sub-Arctic sites have shown wet deposition fluxes of 0.5–2.0 g m⁻² yr⁻¹ (Sanei et al., 2010). Direct comparison between the model and observations is not entirely feasible, as the simulation is not temporally consistent with the observations, thus ignoring seasonal variation in precipitation rates, chemistry, and deposition. Measurement challenges also introduce some uncertainty in these values as the collection efficiency of samplers can be affected by the type of precipitation (rain vs. snow; Prestbo and Gay, 2009; Rasmussen et al., 2012). However, these observations can still serve as a good indicator of general model performance.

Compared to the reported measurements above, the model underestimates Hg wet deposition over all these locations, predicting a negligible amount of Hg wet deposition. One explanation is the seasonal differences in precipitation rates, with trends showing an annual peak in the Arctic during summer, likely contributing to greater Hg wet deposition. During the simulated period, our model predicts low precipitation above 60°N, with a cumulative average of 78 mm. Cumulative snowfall measurements from MOSAiC reported estimates of 72–107 mm, between October 31, 2019, and April 26, 2020, giving evidence to the low precipitation rates in the central Arctic (Wagner et al., 2022). Another possibility is the model implementation of Hg(II) wet deposition, which is currently treated by considering the solubility of Hg(II) equal to that of HNO₃. However, this assumption is a broad simplification and should be addressed in future model developments. The model-predicted average of total Hg deposition in the Arctic (>60°N) is 1.65 g m⁻² (Figure 12c). This value is on the lower range of the recent model-ensemble predictions (10.5 ± 5.0 g m⁻² yr⁻¹, between March and May) of Dastoor et al. (2022a). However, our modelled deposition rates are not truly representative of the entire spring season, as we would expect higher deposition rates in late spring as a result of the transition from high Hg(p) to high Hg(II) in the Arctic (Steffen et al., 2014). There is also a large variability in simulated springtime deposition fluxes between models (e.g., Angot et al., 2016; Dastoor et al., 2022a) due to differences in model performance of simulating AMDEs. Further evaluation of modelled Hg deposition fluxes with springtime measurements are needed, particularly at Arctic sites where observations are limited. Alternatively, a year-long model simulation could provide a better comparison with observations; however, additional model developments would be required (e.g., improved description of Hg(II) and Hg(p) wet deposition) before this simulation can be performed.

Figure 13 shows the simulated mean Hg(0) re-emission flux, based on deposited RM to snow and sea ice. Mean Hg(0) re-emission fluxes up to 2.5 ng m⁻² h⁻¹ are predicted over coastal snow-covered regions, whereas re-emission from sea ice is considerably lower (approximately 0.1 ng m⁻² h⁻¹). In the context of Arctic observations, Hg(0) re-emission fluxes from snow have been reported with large variability, ranging from mean net negative fluxes (e.g., Brooks et al., 2006; Cobbett et al., 2007; Steen et al., 2009; Manca et al., 2013) to mean positive fluxes up to 534 ng m⁻² h⁻¹ during spring (Schroeder et al., 2003; Ferrari et al., 2005; Sommar et al., 2007; Ferrari et al., 2008; Steen et al., 2009; Mann et al., 2015b; Kamp et al., 2018). Direct comparison between reported values and the modelled fluxes is difficult due to differences in measurement techniques, sampling locations, AMDE frequency, and time of year, resulting in different atmospheric conditions and snowpack properties. As a broad assessment, however, the simulated re-emission fluxes here (0–2.5 ng m⁻² h⁻¹) are within the ranges reported by many Arctic and sub-Arctic measurement studies (Dommergue et al., 2003; Schroeder et al., 2003; Ferrari et al., 2005; Sommar et al., 2007; Ferrari et al., 2008; Mann et al., 2015b). Caution should also be taken when evaluating mean fluxes as studies have shown large re-emission fluxes immediately following AMDEs, often several times greater than the average springtime re-emission flux (e.g., Ferrari et al., 2008; Manca et al., 2013; Kamp et al., 2018). A subsequent investigation of Hg(0) re-emission following AMDEs in the model would therefore be useful to assess re-emission flux variability during the spring season, and, potentially complement location-based measurement studies. Measurement data of Hg(0) re-emission from sea ice are even more scarce, making it difficult to assess the model values against observations. Additional measurements of snowpack Hg(0) re-emission fluxes over sea ice are particularly desirable to better evaluate and refine the current Hg(0) re-emission parameterization in the central Arctic.

To understand the transfer of Hg between air and snow in the model, we calculated the fraction of deposited Hg re-emitted from snow and sea ice for the entire simulation period. Above 60°N, approximately 40% of deposited Hg is re-emitted back to the atmosphere, with 60% remaining within the snowpack. Over land-based snow, 74% of deposited Hg is re-emitted, whereas only 4% of deposited Hg over sea ice is released back into the atmosphere. We applied a smaller photoreduction rate (k_red) of RM on sea ice than for land-based snow (see Section 2.1.4), based on the hypothesis that Hg(0) re-emission is comparatively smaller over Arctic sea ice than the continental snowpack (Steffen et al., 2013). Hg(0) observations over Arctic tundra have shown large peaks following depletion which were not observed with the same intensity over sea ice, suggesting lower re-emission rates (Steffen et al., 2013). The presence of chloride in snow has also been hypothesized to increase Hg retention via stabilisation of Hg(II) and a suppression of photoreduction (Poulain et al., 2004; Hintelmann et al., 2007; Lehnherr and St Louis, 2009). Recent experimental evidence has supported this hypothesis, showing a negative relationship between chloride concentrations and the amount of photoreduced Hg in snow (Mann et al., 2018). In the context of the Arctic, snow on sea ice is overall more likely to retain mercury than land-based snow as sea ice regions are typically more enriched with chloride (Krnavek et al., 2012; Peterson et al., 2019). The exact mechanism of Hg stabilization by chloride in snow remains unclear, with theories proposing the formation of photostable chlorocomplexes, or by considering increased Hg(p) concentration in snow to be more stable against photoreduction (Brooks et al., 2006; Hintelmann et al., 2007; Poulain et al., 2007). For this study, we assumed deposited Hg(p) to be equally photoreducible as Hg(II); however, this assumption may need to be revisited as future work refines our knowledge of the fate of Hg in snow.

We also assessed the sensitivity of our assumption that 60% of deposited RM to snow and sea ice is available for re-emission by performing a simulation where 100% of deposited RM is assumed to be photoreducible (Figure 14). Results from this sensitivity test show minimal change to the re-emission fluxes from land-based snow, but re-emissions over sea ice increase by approximately 60–70%. The uncertainty associated with the rate and magnitude of RM photoreduction in snow remains large, motivating more observational and modelling studies to improve our understanding of these processes. Future model development could also aim to refine these parameterizations by testing for different seasons, particularly during summer when Hg(0) re-emissions are known reach their annual maximum (Araujo et al., 2022).

This study presents a comprehensive development of the WRF-Chem model, including Hg gas-phase and aerosol chemistry, deposition, and re-emission processes in the Arctic. The recently improved mercury chemical mechanism of Shah et al. (2021) was added to a version of WRF-Chem that includes polar bromine emissions from surface snow and blowing snow together (Marelle et al., 2021), and tested in a dedicated regional Arctic modelling study for the first time. This new model development allowed us to simulate springtime Arctic ozone and mercury depletion on an hourly timescale, with modelled oxidant concentrations calculated online. Model results were then evaluated with simultaneous measurements of Hg(0), O₃, and BrO from the central Arctic during MOSAiC and from Arctic stations. The main results of this study can be summarized as follows:

• Model predictions of Hg(0) and O₃ show extended periods of depletion in the central Arctic during spring, in agreement with observations from MOSAiC. Our model results also indicate that bromine is the major Hg(0) oxidant in the Arctic, driving springtime Hg(0) depletion.

• Oxidized mercury (Hg(II) and Hg(p)) in the model, evaluated to understand the speciation and spatial distribution of Hg, indicated high aerosol processing of Hg during spring, with Hg(p) accounting for approximately 63% of all Hg(II) species. Measurements of speciated Hg remain sparse and uncertain. A focus on reduced uncertainty in Hg(II) and Hg(p) observations, additional measurement campaigns, and long-term observations are needed to better evaluate the model.

• Both deposition fluxes of oxidized mercury and re-emission fluxes of Hg(0) from snow/ice remain uncertain. Despite this uncertainty, the balance of oxidation, deposition, and re-emission predicted here provides reasonable modelled quantities of gas-phase Hg(0) compared to measurements. Based on experimental evidence, we can expect an average systematic uncertainty of approximately 10%, and in extreme cases up to 20%, for Hg(0) measurements (Slemr et al., 2015).

• For our modelled period, only 4% of deposited Hg(II) and Hg(p) over the Arctic Ocean is re-emitted compared to 96% that remains trapped in snow/ice. Over snow-covered land, the percentage of re-emitted Hg(0) is higher (74% re-emitted vs. 26% retained in the snowpack). The implications of these findings for Arctic Hg(0) summertime re-emission (Araujo et al., 2022), and under long-term environmental changes (e.g., sea ice loss), are important to consider in future work.

This work serves as a basis for future studies to explore some key questions regarding Arctic Hg chemistry. For example, the contribution of the Arctic Ocean to summertime Hg(0) re-emission is an outstanding research question which could be tested using this model. Additionally, the broader impacts of future climate scenarios on Hg chemistry, emissions, and deposition in the Arctic may also be assessed. Overall, this work could enable better predictions of the long-term implications of climate change on Hg contamination of Arctic ecosystems. Finally, we note below the key processes for mercury included in the model that remain uncertain:

• Chemical kinetics—Identification of Hg(II) species is a top priority to improve our understanding of mercury redox chemistry. This goal requires additional theoretical, experimental, and modelling studies to reduce uncertainties in the reaction rates of mercury. In particular, photoreduction of Hg(II) and Hg(p) should be investigated further to better constrain the contribution of reduction kinetics to atmospheric Hg chemistry.

• Dry and wet deposition—Model descriptions of Hg dry and wet deposition should be revised in future work to more accurately determine the transfer of Hg to snow and ice surfaces.

• Gas-particle partitioning—Better understanding of the main model parameters (e.g., heterogeneous uptake rate, partitioning coefficient) that control Hg gas-particle partitioning is needed. An evaluation of modelled aerosol concentrations using observations is needed to refine bromine and mercury heterogeneous recycling on aerosols.

• Re-emission fluxes from snow and snow on sea ice—Hg(0) emissions from snow and ice remain difficult to constrain in models due to the complexity of processes that contribute to them and the difficulty of measuring fluxes. Model assumptions (e.g., 60% reducible Hg in snow, photoreduction rate constants from snow and snow on sea ice) need additional testing. A more accurate representation of factors which influence Hg(0) re-emission (e.g., snowpack chloride concentration) will also likely be needed in future work.

The updated WRF-Chem model code used in this study is available on Zenodo at Ahmed et al. (2022), https://doi.org/10.5281/zenodo.7137482. The surface meteorology measurements made onboard the Polarstern can be accessed at Schmithüsen (2021). The MOSAiC radiosonde data are available at Maturilli et al. (2021). The elemental mercury dataset from MOSAiC is available at Angot et al. (2022a). The merged ozone dataset from MOSAiC is available at Angot et al. (2022a). BrO MAX-DOAS measurement data are available at Mahajan (2022). Surface elemental mercury and ozone data at Villum (Greenland) and Zeppelin (Svalbard) are available at the EBAS website https://ebas.nilu.no/. The surface total gaseous mercury measurement dataset at Alert (Canada) is available at https://data-donnees.ec.gc.ca/data/air/monitor/monitoring-of-atmospheric-gases/total-gaseous-mercury-tgm/. The speciated mercury data at Alert (Canada) are available at https://data-donnees.ec.gc.ca/data/air/monitor/monitoring-of-combined-atmospheric-gases-and-particles/speciated-mercury. Surface ozone measurements from Alert (Canada) are provided from the Canadian Air and Precipitation Monitoring Network (CAPMoN) and are available at https://data-donnees.ec.gc.ca/data/air/monitor/monitoring-of-atmospheric-gases/ground-level-ozone/. Surface ozone data at Utqia͘gvik (Alaska), Summit (Greenland), and Eureka (Canada) were provided by the NOAA Global Monitoring Laboratory (https://www.esrl.noaa.gov/gmd/ozwv/surfoz/data.html). Satellite BrO VCD data from TROPOMI are available on Zenodo at Zilker and Richter (2023), https://doi.org/10.5281/zenodo.7711532.

The supplemental files for this article can be found as follows:

Figure S1–S11. Table S1. Available as a merged PDF.

We would like to thank the following individuals for their contributions during MOSAiC: Jacques Hueber, Matthew D. Shupe, Marion Maturilli, and Holger Schmithüsen. We also thank all those who contributed to MOSAiC and made this endeavour possible (Nixdorf et al., 2021). We thank all those involved in the operation and maintenance of the Polarstern research vessel (Knust, 2017). Computer analyses benefited from access to IDRIS HPC resources (GENCI allocation A011017141) and the IPSL mesoscale computing center. Some observations used here were provided by the Atmospheric Radiation Measurement (ARM) Climate Research Facility, a US Department of Energy (DOE) Office of Science User Facility sponsored by the Office of Biological and Environmental Research. We thank the many ARM operators who supported the field observations. We also thank the Dr. Neil Trivett Global Atmosphere Watch Observatory at Alert for technical support.

This work was supported by the Ecole Doctorale Sciences de la Terre, de l’Environnement et des Planètes (ED105) of Université Grenoble Alpes. JLT is funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 101003826 via project CRiceS (Climate Relevant interactions and feedbacks: the key role of sea ice and Snow in the polar and global climate system). JLT and SA also acknowledge support by the CNRS INSU LEFE-CHAT program under the grant BROM-ARC. Some of the observational data reported in this manuscript were produced as part of the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition with the tag MOSAiC20192020, with activities supported by Polarstern expedition AWI_PS122_00. Data acquisition during the MOSAiC expedition was funded by the US National Science Foundation (NSF) (award OPP 1914781 and award OPP 1807163), the Swiss National Science Foundation (grant 200021_188478), the Swiss Polar Institute, the DOE Atmospheric System Research Program (DE-SC0019251), and the US National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Laboratory. Part of this project was also funded by the European Research Council (ERC), H2020 European Research Council (GASPARCON; grant No. 714621), and the Academy of Finland Flagship funding (grant No. 337552). JS holds the Ingvar Kamprad chair for extreme environments research, sponsored by Ferring Pharmaceuticals. The measurements at Villum Research Station was financially supported by the Danish Environmental Protection Agency with Funds for Environmental Support to the Arctic Region (Project No. J.nr. 2021 – 60333). Measurements in Alert are supported by the Northern Contaminants Program (Crown-Indigenous Relations and Northern Affairs Canada). Satellite BrO data analysis was supported by the Deutsche Forschungsgemeinschaft (project no. 268020496 – TRR 172) within the Transregional Collaborative Research Center “ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC)3” in subproject C03. AR at AWI was partly supported by the European Union’s Horizon 2020 research and innovation framework programme under Grant agreement no. 101003590 (PolarRES).

Detlev Helmig (Editor-in-Chief), Hélène Angot (Associate Editor), and Byron Blomquist (Associate Editor) are part of the Editorial board of Elementa: Science of the Anthropocene. They were not involved at any step in the review process of this manuscript.

Conceptualization, design, and initial draft of the manuscript: SA, JLT, ADo, HA, RL, JSt.

Model developments: SA, JLT, LM.

Model simulations: SA.

Data acquisition: HA, SDA, A-MB, BB, IB, LB, MB, NB, JHC, ADa, SD, DHo, DHe, H-WJ, TJ, TL, ASM, KP, KAP, LLJQ, ARin, ARy, ARic, AS, AS-L, GS, HS, JSc, OT, BZ.

Editing and revision of the manuscript: All authors.

Approval of the final submitted manuscript for publication: All authors.