Evaluating atmospheric mercury (Hg) uptake by vegetation in a chemistry-transport model

Mercury (Hg), a neurotoxic heavy metal, is transferred to marine and terrestrial ecosystems through atmospheric transport. Recent studies have highlighted the role of vegetation uptake as a sink for atmospheric elemental mercury (Hg0) and a source of Hg to soils. However, the global magnitude of the Hg0 vegetation uptake flux is highly uncertain, with estimates ranging 1000–4000 Mg per year. To constrain this sink, we compare simulations in the chemical transport model GEOS-Chem with a compiled database of litterfall, throughfall, and flux tower measurements from 93 forested sites. The prior version of GEOS-Chem predicts median Hg0 dry deposition velocities similar to litterfall measurements from Northern hemisphere temperate and boreal forests (∼0.03 cm s−1), yet it underestimates measurements from a flux tower study (0.04 cm s−1vs. 0.07 cm s−1) and Amazon litterfall (0.05 cm s−1vs. 0.17 cm s−1). After revising the Hg0 reactivity within the dry deposition parametrization to match flux tower and Amazon measurements, GEOS-Chem displays improved agreement with the seasonality of atmospheric Hg0 observations in the Northern midlatitudes. Additionally, the modelled bias in Hg0 concentrations in South America decreases from +0.21 ng m−3 to +0.05 ng m−3. We calculate a global flux of Hg0 dry deposition to land of 2276 Mg per year, approximately double previous model estimates. The Amazon rainforest contributes 29% of the total Hg0 land sink, yet continued deforestation and climate change threatens the rainforest's stability and thus its role as an important Hg sink. In an illustrative worst-case scenario where the Amazon is completely converted to savannah, GEOS-Chem predicts that an additional 283 Mg Hg per year would deposit to the ocean, where it can bioaccumulate in the marine food chain. Biosphere–atmosphere interactions thus play a crucial role in global Hg cycling and should be considered in assessments of future Hg pollution.

The cuticular resistance is calculated as: where -+, " ! is a reference cuticular resistance for the land type " which is scaled by the solubility and reactivity factors. The aerodynamic resistance in the lower part of the canopy is calculated as a function of the solar radiation (?): -!-" = 100 × @1 + 1000 ? + 10 A (6) The lower canopy surface resistance is calculated by scaling the resistances of the reference compounds sulfur dioxide (SO2, which is soluble yet unreactive) and ozone (O3, which is reactive yet insoluble) by the solubility and reactivity of Hg 0 : 1 --+ " = : 24 # * : 73 " * 1 --+,73 " " + < 6 24 # 1 --+,3 $ " The ground surface resistance is calculated in a similar manner as the lower canopy resistance: 1 To account for reduced uptake by surfaces under cold temperatures (Wesely, 1989), an additional resistance (-: ) is calculated as a function of temperature: This resistance is used to correct resistances for the ground surface, lower canopy, and cuticular surfaces, with the maximum effect capped to a factor of 2 times the original resistances (Jaeglé et al., 2018): In addition to the temperature correction, the cuticular resistance is scaled by leaf area index:

S2. Validation of offline dry deposition model with online results
The offline model (code can be found at: https://github.com/arifein/offline-drydep) for Hg 0 dry deposition velocities is compared with online GEOS-Chem simulations to validate its use in this study. Tests were made at different time resolutions of the offline model to identify a suitable time resolution that balances computational cost of the offline simulation with accuracy for predicting online simulated deposition velocities (Table S1). When using daily resolution of meteorological inputs, the mean error of the offline model over land grid cells is 15.2%. The large error illustrates the importance of accurately capturing the diurnal cycle of Hg 0 dry deposition. When the time resolution of the offline model is set to hourly, the mean error of land grid cells shrinks to 0.1%. Since we use two significant digits to report observed and simulated deposition velocities, we consider this error to be sufficiently small to apply the offline model in the current study. The offline model shows mean errors of similar relative magnitude (~0.1%) for f0 = 10 -5 ( Figure S1) and f0 = 1 ( Figure S2), and for January ( Figure S1) and July ( Figure S3).    is shown in (c). Hourly time resolution was used for the offline model.

S3. Effect of land-surface adjustment on dry deposition evaluation
We follow the approach of Silva and Heald (2018) for ozone to account for the actual land category at the observation sites in the offline model calculations of Hg 0 dry deposition velocities. In the original model predictions, the model grid cell can contain multiple land cover types (e.g., ocean, deciduous forest, coniferous forest, grassland, etc.). Our adjusted predictions use only the observed land cover part of the grid cell to calculate the modelled dry deposition velocity. Figure    illustrating that much of the variability in the original model predictions for f0 = 10 -5 is due to variability in the fraction of grid box covered by forest. Figure S5. Comparison between measured litterfall deposition velocities and modelled Hg 0 dry deposition velocities, without and with land type adjustment for f0 = 10 -5 .

S4. Three-box model tuning of GEOS-Chem photoreduction
Due to its large uncertainty, the photoreduction rate of Hg 2+ in organic aerosol and clouds is used as a tuning parameter in the GEOS-Chem Hg simulation (Horowitz et al., 2017;Shah et al., 2021). Upon changes to the simulation chemical mechanism or ocean setup, the reduction rate coefficient (S) is adjusted to give the best overall agreement with observed atmospheric Hg concentrations. Common practice has been to conduct GEOS-Chem Hg simulations with S adjusted ±10% until optimum agreement is achieved (http://wiki.seas.harvard.edu/geos-chem/index.php/Mercury#K_RED_JNO2).
However, this approach requires multiple tuning runs of GEOS-Chem each time a parameter is changed within the Hg simulation.
To reduce the computational expense of the tuning procedure, we use a three-box model of the  We illustrate our GEOS-Chem tuning procedure through an example from the current study: 1) From the BASE simulation results, calculate the burdens and rate coefficients needed for the three-box model ( Figure S6).
2) In the OBRIST simulation, the dry deposition of Hg 0 is enhanced, reducing the overall lifetime of Hg 0 and the burden. Thus, the Northern Hemisphere Hg 0 burden becomes too low compared to observations ( Figure S9). We want the Northern Hemisphere Hg 0 burden to be roughly equivalent to BASE, which matches well with the annual mean of atmospheric Hg observations ( Figure S9).
3) From the OBRIST simulation, we calculate the new rate coefficients of Hg 0 dry deposition to land (D0SSH and D0SNH) and implement these in the box model. 4) Using the box model, we run simulations at different tropospheric reduction rates of Hg 2+ ( Figure   S7a). We identify the value of the tropospheric reduction rate (102 yr -1 ) needed to yield the same burden of Hg 0 in the Northern Hemisphere as the BASE simulation. We focus on the Northern Hemisphere since the bulk of the atmospheric Hg measurements used to tune the model are located in the Northern Hemisphere ( Figure S9).
5) The actual tuning parameter in GEOS-Chem is S, the reduction rate coefficient. Using three previous Hg runs in GEOS-Chem, we identify a logarithmic relationship between S and the reduction rate in the Northern Hemisphere ( Figure S7b). Using the fitted logarithmic curve, we find required reduction rate (102 yr -1 ) corresponds to S = 0.33. 6) We run the OBRIST_R simulation with S set to 0.33.

7)
We verified that the Northern Hemisphere tropospheric Hg 0 burden in the GEOS-Chem OBRIST_R simulation (1990 Mg) is roughly equivalent to the BASE simulation (1959 Mg).
Also, we verified that the R 2 and bias of the OBRIST_R simulation compared to benchmark atmospheric Hg concentrations is similar to the BASE simulation ( Figure S9). Using results from three previous runs of GEOS-Chem with different S parameters used, we fit a logarithmic relationship. This logarithmic relationship is used to calculate the value of S (0.33) that would yield the required rate of reduction in the Northern Hemisphere (102 yr -1 ).

S5. Alternative solutions for matching Amazon dry deposition observations
The Hg 0 dry deposition velocity in the Amazon is underestimated by BASE and OBRIST simulations, which use a single parameter for the Hg 0 biological reactivity (f0) in the dry deposition scheme (Section 3.1). We investigated whether the model also shows a bias in terms of ozone dry deposition velocities in the Amazon. In Figure S8, we compare the offline modelled ozone dry deposition velocities for 2015 with measurements in the Amazon from Fan et al. (1990) and Rummel et al. (2007), compiled in Silva and Heald (2018). The model underpredicts ozone dry deposition velocities by 39%, showing mean velocities of 0.49 cm s -1 compared to the observed 0.80 cm s -1 . For Hg 0 , we compare model estimates using f0 = 3 × 10 -5 (OBRIST setting) with the mean of ten Amazon litterfall measurements compiled by Fostier et al. (2015). The mean modelled dry deposition velocity estimate is 0.098 cm s -1 , which is 53% below the mean measured value of 0.21 cm s -1 . Thus, the model appears to show a negative bias, of similar relative magnitude, for both Hg 0 and ozone dry deposition in the Amazon. It is important to note, however, that only three measurement campaigns are available for the Amazon ozone dry deposition velocity estimate.
Given the similarity between both ozone and Hg 0 dry deposition biases, we explored resistance parameter solutions that could improve model-measurement agreement for both compounds. Initial tests of the GEOS-Chem parametrization identified two influential parameters for the ozone dry deposition velocity: the internal stomatal resistance (RI, called -"* " ! in Eq. 4) and the cuticular resistance (RLU, called in Eq. 14). For the rainforest land class in GEOS-Chem, the parametrization uses RI = 200 and RLU = 1000. Figure S8 illustrates two scenarios that match the measured mean ozone dry deposition velocity: when RI = 65 or when RI = 100 and RLU = 500. When these parameter combinations are used to calculate the mean Hg 0 dry deposition velocities in the Amazon, the model predicts values of 0.12 and 0.11 cm s -1 , below the observed mean of 0.21 cm s -1 . By decreasing RI to 1, the Hg 0 dry deposition velocity (0.19 cm s -1 ) approaches the observed value. However, the predicted ozone dry deposition velocity from this parameter scenario (2.1 cm s -1 ) strongly overestimates the observed value (0.80 cm s -1 ).

Figure S8.
Comparing observed deposition velocities in the Amazon rainforest with offline model calculations using different estimates for internal stomatal resistance (RI) and cuticular resistance (RLU).

Markers indicate mean value over observation locations and error bars indicate maximum and minimum
values. Amazon ozone measurements are taken from Fan et al. (1990) and Rummel et al. (2007) and Amazon Hg litterfall measurements are compiled in Fostier et al. (2015). Elemental mercury dry deposition velocities are calculated assuming f0 = 3 × 10 -5 .
In conclusion, although resistance parameters can be tuned to match the Hg 0 litterfall dry deposition velocities in the Amazon, this approach may yield unreasonable dry deposition velocities for other chemical compounds. Therefore, we proceeded with our approach of creating an additional compound-specific parameter (f0) for the Amazon rainforest region. The need for a new parameter may suggest that Amazon tree species take up Hg 0 at a higher rate than other trees, which would not be surprising given inter-species differences in Hg 0 uptake found in a European study (Wohlgemuth et al., 2022). Additional Hg 0 uptake measurements in the Amazon region would be needed to assess whether a species-specific vegetation uptake scheme could be applied in GEOS-Chem. Ongoing improvements to model dry deposition schemes (e.g., Clifton et al., 2020;Lin et al., 2019) may require future retuning of the Hg 0 biological reactivity to available vegetation uptake measurements.

S6. Additional model evaluation plots and table
Additional comparisons between simulations and observations are shown in Figure S9 for global surface TGM concentrations, Figure S10 for the overall global wet deposition comparison, Figure S11 for North American Hg wet deposition, and Figure S12 for European wet deposition. Table S2 lists the global Hg burdens and major fluxes for all simulations.

S7. Influence of meteorological year on modelled deposition
To test the sensitivity of our offline dry deposition results to the choice of year for the meteorological and LAI data, we ran the offline model under the AMAZON_U scenario using data from three different years (2014)(2015)(2016). The median model predictions for the different datasets do not show substantial changes between the three meteorological years. Due to these sensitivity tests, we have used a single simulation year (2015), since comparable input data for leaf area input is not available over the entire observational time period. This result agrees with other studies that have determined GEOS-Chem modelled ozone dry deposition variability due to meteorology to be generally small (<5%, Silva and Heald, 2018). However, the simulated dry deposition variability has been found to be smaller than observed dry deposition variability, which could be due to the lack of moisture availability as a limiting factor in the dry deposition scheme (Clifton et al., 2020;Lin et al., 2019).