Mercury in soils of the conterminous United States: patterns and pools

Soils account for the largest global mercury reservoirs, but observations are sparse in many regions. The accumulation and turnover of mercury in soils determines whether they act as an atmospheric source or sink. Here, we present a spatial analysis of soil mercury from a large soil survey (three horizons, ∼4800 sites) across the conterminous United States conducted by the U.S. Geological Survey. Soil mercury pools were calculated for 11 layers, cumulatively representing the top 1 m of soil, and totaling 158 ± 2 Gg (±SD) of mercury (20.3 ± 0.2 mg m−2). Mercury areal density was greatest in mixed forest (27.3 ± 0.5 mg m−2), cropland (25.3 ± 0.3 mg m−2), and deciduous forest (25.6 ± 0.5 mg m−2) ecosystems and lowest in barren (13.5 ± 0.3 mg m−2) and shrubland (12.6 ± 0.2 mg m−2) ecosystems. Assessment of the provenance of soil mercury using bedrock titanium normalization suggests that 62%–95% of soil mercury is unexplained by parental sources.


Introduction
Mercury is a ubiquitous and toxic metal released from natural and anthropogenic sources. The lifetime of mercury in biologically relevant environments is extended by the re-emission of previously deposited mercury back to the atmosphere, where cycling can continue until mercury is eventually sequestered in soils or marine sediments [1,2]. Within the mercury cycle, soils function as a critical medium, by serving as both an atmospheric sink and source via deposition and re-emission, and as an aquatic source via leaching [3,4]. Soils have received considerable attention in the scientific literature [5], but there remains a paucity of robust global and regional soil reservoir estimates. Given the importance of soil pools within the broader global mercury cycle, comprehensive and quantitative estimates of soil pools are needed.
The scarcity of regional soil pool estimates is largely due to the heterogeneity of soils and the laborious nature of their sampling, both of which pose challenges to obtaining accurate pool estimates, especially as spatial scale increases. Consequently, adequately sampled soil data, which are critical for scientific studies, policy making, and remediation efforts, are rarely available beyond the local level, and many studies either extrapolate from smaller datasets or rely on more complex modeling approaches to estimate soil data [2,[6][7][8][9][10][11]. Estimates are most often based, at least partially, on soil organic carbon (SOC) data, which correlate strongly with soil mercury concentrations and are more widely available [10][11][12][13]. Other predictors often employed include climatic factors, ecosystem properties, latitude, grain size, and proximity to anthropogenic mercury sources [13][14][15][16]. The complex interactions among these characteristics can often be more easily captured in broader terms of landcover type or ecoregion [13,15], making them useful for assessing larger soil spatial patterns.
Soil pool estimates can be difficult to compare because they are often delineated differently among studies. Previous studies have delineated pools by horizon (organic vs mineral), mercury source (lithogenic vs excess), turnover time (fast, slow, armored), and depth. Global pools are estimated between 240 and 1088 Gg, but corresponding soil depths vary by study from 10 to 40 cm [6,7,9,11,15]. Robust soil pool estimates for the European Union [16] and the boreal/arctic [8,10,11] help to better constrain current global estimates, but high quality, regional estimates are needed for the United States and other regions.
Beginning in 2007, the U.S. Geological Survey (USGS) collected soil samples from 4857 sites across the conterminous United States (CONUS) and analyzed them for 45 elements, including mercury [17]. The large sample size and randomized sampling design make the dataset useful for investigating the quantity and spatial distribution of mercury in soils across CONUS. The main objectives of this study were to: (a) quantify how mercury distributions vary with depth, land cover, and ecoregion, (b) develop an empirical model using soil and climate data to predict mercury concentrations in the top meter of CONUS soils, and (c) use the predicted concentrations to quantify the soil mercury reservoir in the top meter of CONUS. We hypothesize that external inputs are important in the supply of soil mercury and that strong associations between mercury and SOC yield higher mercury concentrations and pools in the surface soils of ecoregions and land use types with greater carbon inputs (e.g. forested ecosystems).

Data collection
A detailed description of the methodologies used for soil sampling and chemical analysis can be found in Smith et al [17]. Briefly, a generalized random tessellation stratified design was used to select 4857 sites (approximately 1 site/1600 km 2 ; figure S1). Soil was collected from the top 5 cm (∼2 kg), the A horizon (∼2 kg), and the C horizon (∼1 kg) at each site from 2007 to 2010. The top 5 cm samples excluded fresh or recognizable litter layers (L or F horizons) but could include both organic and mineral soil. In cases where the C horizon occurred deeper than 100 cm, a surrogate sample was collected from 80 to 100 cm instead [17]. The mean depth interval (i.e. mean of the top and bottom depth of each sample) was ∼0-17 cm for the A horizon and was ∼65-87 cm for the C horizon. Samples were air-dried, sieved to <2 mm, crushed to <150 µm, and analyzed for mercury, organic carbon, and titanium (Ti), among other elements. Mercury was measured using cold vapor atomic absorption spectrometry in accordance with US Environmental Protection Agency (USEPA) method 7471b (revision 2). The lower limit of determination (LLD) for mercury was 10 ng g −1 and of the 14 434 samples, there were 1505 non-detects. Organic carbon was calculated for the A and C horizon samples by measuring total carbon using an automated carbon analyzer and subtracting inorganic carbon, which was calculated using the carbonate mineralogy determined in the survey (see Smith et al). Organic carbon was not calculated for the top 5 cm. The organic carbon LLD was 0.01% and of the 9587 samples, 124 were below detection. Ti was analyzed via inductively coupled plasma-atomic emission spectrometry and had an LLD of 0.01%. Of the 14 434 Ti samples, four samples were below the LLD for Ti. An overview of the relevant quality assurance/quality control (QA/QC) data from Smith et al [17] can be found in the supplementary material.
Bulk density and soil organic matter coverages (3 km × 3 km) for CONUS were generated from the United States Department of Agriculture (USDA) Gridded National Soil Survey Geographic Database (gNATSGO) database. SOC layers were calculated by multiplying the soil organic matter coverages by 0.58, the mass fraction of organic matter assumed to be carbon [18]. Rock fragment volume data (>2 mm, 3 km × 3 km) were obtained from the Soil Information for Environmental Modeling and Ecosystem Management (Pennsylvania State University), which provided coverages generated from USDA State Soil Geographic Database (STATSGO) data [19]. Coverages consisted of 11 layers, cumulatively representing the upper 1 m of soils (0-5 cm, 5-10 cm, 10 cm intervals thereafter). Precipitation and temperature data (0.25 • × 0.3125 • ) were obtained from the GEOS-5 FP meteorology product from the NASA Global Modeling and Assimilation Office [20,21]. Leaf area index (LAI) data were from MODIS Collection 5 LAI (2005-2016) reprocessed according to Yuan et al [22]. USGS Elevation data (30 m × 30 m) were obtained through ESRI ArcGIS Image service [23]. Level 1 ecoregion spatial data were obtained from the USEPA. Landcover data was obtained from the USGS 2011 national land cover database (NLCD). See the Data Availability section for data access.

Horizon, ecoregion, and landcover statistical analysis
Mercury concentrations and pools were calculated for each horizon, ecoregion, and landcover type. Soil horizons were included as they are commonly used delineators of soil development and were the basis of sampling for the Smith et al [17] survey. Ecoregion and landcover type were included as they capture broad physical and ecological characteristics that likely influence soil mercury distribution. Level I ecoregions are the broadest category in spatial scale and were chosen to match the predictive scale of this study. While both ecoregion and landcover type provide information on vegetation and habitat type, the former also accounts for other qualities, such as soil type, biological communities, climate, etc. For brevity, only landcover results will be presented in the main text; all ecoregion results can be found in supplementary material. Wetland and aquatic landcover types were excluded from landcover analysis, as they are poorly represented in the soil coverages used for pool estimates. The ecoregions and landcover types used in the present study can be found in the supplementary material (ca. tables S1-S4, figures S8 and S9).
Descriptive and inferential statistics for soil horizons, ecoregions, and landcover types were calculated using the NADA package in R (non-detects and data analysis for environmental data) [24]. This approach was previously used by Obrist et al [13] on a subset of these data to account for the relatively high LLD and consequent high number of non-detects. We used a similar approach of maximum likelihood estimators, though outliers were not removed prior to statistical analyses as they were in Obrist et al [13]. Differences in soil mercury concentrations among soil horizons, ecoregions, and landcover types were assessed using the NADA package equivalent of the Peto-Peto test, a non-parametric test for comparing empirical cumulative distribution functions with left-censored data. For those factors with significant effects on mercury concentrations, post-hoc analysis was performed using the same statistical test with a Bonferroni correction to account for multiple comparisons.

Depth profile interpolation
The ancillary soil data (i.e. bulk density, rock fraction, SOC) used for pool calculations were uniformly gridded (3 km × 3 km). However, the soil sampling methodology used in Smith et al [17] was partially based on depth (0-5 cm) and partially on qualitative characteristics (A and C horizon). To develop corresponding mercury and carbon values that were compatible with the other soil data, we used equal area quadratic smoothing splines (EAQSS) to create continuous mercury profiles at each sampling site. EAQSS have been shown to be effective at producing soil profiles from discrete soil collections [25,26]. Splines were calculated using the R package 'GSIF' and applied for every sample for mercury and organic carbon. For mercury, samples with overlapping top 5 cm and A horizons were amended so that the A horizon began below the top 5 cm. Following the construction of continuous mercury and organic carbon profiles, mean concentrations were calculated for each of the 11 depths present in the STATSGO coverages (0-5 cm, 5-10 cm, 10-20 cm, 20-30 cm, etc).

Mercury concentration modeling and soil pool calculation
Mercury concentrations for each soil layer were modeled using a generalized additive model (GAM) following equation (1): where Hg is soil mercury concentration (µg g −1 ), α is the model intercept, Lat is latitude (degrees) Lon is longitude (degrees), and ε is the error term. GAMs are akin to generalized linear models but use the sum of smoothing functions (here, thin plate regression splines; denoted by s i ) in place of a linear predictor for each covariate, allowing for flexible but interpretable modeling of nonlinear covariate effects. Potential covariates were selected from those which are well characterized in the model domain and have previously been found to influence soil mercury content [7,14,16], which included SOC, latitude, longitude, precipitation, elevation, temperature, and LAI. Variables which fit both criteria but added no additional explanatory power to the model or resulted in high model concurvity were excluded. The model was generated for each soil layer using layer-specific mercury and SOC values. The models were then applied to the predictive coverages to generate mercury concentration surfaces for each soil layer. Outliers, defined as mercury concentrations exceeding 0.25 µg g −1 , were removed prior to modeling. These values may represent local areas of direct anthropogenic mercury release or natural enrichment that the applied model lacks the predictors to recreate (figure S2).
Mercury and SOC pools were calculated for each soil layer (L) following equation (2a): where P L is the mercury or carbon pool (Gg, Pg), C L is the concentration of mercury or SOC (µg g −1 , %), BD L is bulk density (g cm −3 ), RC L is the fraction of soil that consists of rock n > 2 mm in size, and D L is the depth of the layer (cm). Pools calculated for each of the 11 soil layers (P L ) were then summed to produce a composite estimate (P) for the upper 1 m of soils The scale of applicability for the calculated pools was 3 km × 3 km, reflecting the coarsest predictive coverage (i.e. bulk density, SOC). For purposes of analysis and regional comparison, pools were also calculated individually by level I ecoregion (USEPA, figure S8) and landcover (USGS 2011 NLCD, figure S9). Note that Smith et al [17] and the NLCD 2011 data use the same primary landcover categories, but different secondary landcover categories. Uncertainty in mercury pool estimates was assessed through a bootstrap analysis in which sites (n = 4671) were randomly resampled with replacement (i.e. a sample could be selected more than once per resampling), and pools were estimated as described above. This process was repeated 10 000 times to produce quantile, standard deviation, and standard error estimates. Additional uncertainty exists surrounding the physiochemical soil properties used to calculate soil mercury pools, namely SOC, bulk density, and rock fraction. These are highly heterogeneous properties and at the scale of the current work (3 km × 3 km) there is considerable uncertainty in their values. These limitations and their subsequent effect on application of the soil mercury reservoirs presented here are discussed further in section 3.7.

Soil mercury sourcing
Three candidate approaches were considered to characterize the source of mercury in soil pools (for details on the candidate approaches and their assessment, see supplementary material). After assessing the suitability of each approach, bedrock Hg/Ti ratios were selected to calculate the amount of parent-derived mercury expected in the soil profile, based on parent rock weathering alone. This approach assumes that Ti is primarily derived from parent material (regolith) and conserved throughout the soil profile following weathering. This approach has been previously applied for source apportionment of soil mercury [27][28][29].
The bedrock Hg/Ti ratios used in this study were from Richardson et al [29], who analyzed Hg/Ti in deep regolith from six United States Critical Zone Observatories. The 10th and 90th percentile ratios from samples greater than 300 cm were used to calculate background pools; ratios from Eel Creek were not included, as the Hg content is substantially enriched in the local geology [29]. The Ti concentration observed in each horizon (Ti obs ) was scaled by the bedrock Hg/Ti ratio to estimate parent-derived mercury concentrations (Hg parent ): Parent-derived mercury concentrations were then subtracted from the measured mercury concentration at each site (Hg obs ) and the remainder was deemed excess (non-parental or exogenic; Hg non-parental ): Hg non-parental = Hg obs − Hg parent .
Excess mercury values for the three horizons were then splined and used to generate excess mercury soil pools, as detailed above (see section 2.4); zero and negative values produced through this process were replaced with 0.001 µg g −1 to allow for logtransformation without substantially contributing to pool magnitude estimates.

Patterns in measured soil mercury concentrations with depth and landcover
The distribution of mercury in soils is a function of both the magnitude of mercury inputs (e.g. wet and dry deposition, regolith weathering) and the climatic, vegetation, and soil processes governing mercury retention and release. Across CONUS, soil mercury concentrations were highly variable and decreased significantly with depth (χ 2 = 273, df = 2, p < 0.05).
Concentrations (mean ± SE) in the upper 5 cm (0.035 ± 0.0005 µg g −1 ) did not differ significantly with the A horizon (0.034 ± 0.0005 µg g −1 ; χ 2 = 4.1, df = 1, p < 0.04), but both were significantly higher than in the C horizon (0.027 ± 0.0004 µg g −1 ; χ 2 = 237, 180, df = 1, p < 0.01). This enrichment of surface soils supports the premise that external inputs such as atmospheric deposition drive spatial variation in soil mercury concentrations. Landcover has been shown to affect both deposition and re-emission of soil mercury and therefore has important implications for soil mercury concentrations [30,31]. We found landcover had a significant effect on the concentration of mercury in soil samples for all three horizons (χ 2 = 1171 for the top 5 cm, 932 for the A horizon, 284 for the C horizon, df = 6, p < 0.05). Excluding the non-natural woody landcover type (0.011 µg g −1 ), which was highly variable due to low sample count, the barren (0.014 µg g −1 ) and shrublands (0.016 µg g −1 ) landcover types had the lowest median concentrations of mercury for all horizons (supplementary material table S1). This pattern is likely due to low organic carbon input, low quantities of precipitation, and high incident solar radiation, which drives photoreduction in surface soils. Conversely, forested uplands had a significantly higher median mercury concentration than all other ecoregions for the top 5 cm (0.043 µg g −1 ) and A horizons (0.040 µg g −1 ); the median C horizon mercury concentration in forested uplands (0.024 µg g −1 ) was significantly higher than all ecoregions except in the developed landcover type (0.025 µg g −1 ; supplementary material table S1).
Compared to several non-forested landcover types (i.e. shrublands, herbaceous uplands, barren, non-natural woody), mercury concentrations were two times greater in forest surface soils. The enrichment of forested soils emphasizes the important role that forest canopies play in mediating mercury deposition through foliar uptake, litterfall and throughfall [31][32][33], and control of re-emissions [34]. In addition to magnitude, forested uplands showed the greatest difference in median concentration between upper and lower soils. Both patterns may be explained by increased mercury input through litterfall and throughfall, decreased photoreduction in surface soils by canopy cover, and from higher forest soil SOC content-which is shown to retain mercury in soils.
While forested ecosystems have been previously recognized for their high soil mercury concentrations [13], developed (i.e. urbanized) and planted/cultivated landcover types have received less attention and thus their relatively high concentrations, particularly at depth, are notable (supplementary material table S1). Both landcover types are defined by anthropogenic alterations and are exposed to very different conditions than 'natural'  [17] and were modeled based on latitude, longitude, and soil organic carbon. Note: Color scaling is not at regular intervals. Scale of applicability is 3 km × 3 km. soils (e.g. tilling, construction, impervious surfaces), which affect the magnitude and distribution of soil mercury. Median mercury concentrations in developed soils ranged from 0.030 µg g −1 in the top 5 cm to 0.025 µg g −1 in the C horizon and were significantly higher for all horizons than in planted/ cultivated soils, in which concentrations ranged from 0.023 µg g −1 in the top 5 cm to 0.019 µg g −1 in the C horizon. One explanation for the elevated mercury concentrations in these landcover types is their proximity to population centers, and thus mercury-releasing processes (e.g. coal combustion, waste incineration). Additionally, contemporary and legacy inputs, including direct inputs of mercurycontaining products (e.g. fungicides, fertilizers), may be accompanied by alteration of the physical and chemical characteristics of soil, further impacting mercury retention and release dynamics.

Spatial modeling of mercury and patterns in modeled values
Mercury concentration modeling via GAM produced concentration surfaces for all 11 soil layers (see figure 1 for top 5 cm). Mean (±SD) concentrations ranged from 0.033 µg g −1 (±0.032) in the top 5 cm to 0.021 µg g −1 (±0.010) in the bottom 10 cm (supplementary material figure S4). Model fitness decreased with depth, from an adjusted R 2 of 0.58 for the top layer (0-5 cm) to an adjusted R 2 of 0.37 for the bottom layer (90-100 cm). This pattern is due in part to decreasing sample size with depth, but also because the surface processes driving the model are increasingly less explanatory with depth (supplementary material figure S5). The controls on mercury in deeper soils are likely related more to the physical and chemical characteristics of the soil than to the processes governing mercury input. However, spatial data detailing these characteristics were unavailable and could not, consequently, be incorporated into the empirical model.
Modeled mercury concentrations were variable among landcover types and relative patterns generally matched those of SOC concentrations, with forested regions being highest and shrub/barren lowest (figure 2). Mixed and coniferous forests soils had higher surface concentrations than deciduous forests, consistent with observations from Ballabio et al [16] and Richardson and Friedland, which the later suggests could result from differences in litter quality, microbial degradation rates, and SOM sorption capacity [35]. For all landcover types, concentrations of mercury were highest in surface soils and decreased with depth, most notably within the top 20 cm. Below this depth, soils concentrations exhibited low variability, suggesting that while litter and throughfall are important determinants of mercury concentration in surface soils, soil characteristics are more important for deeper soil mercury concentrations [35]. Croplands were unusual in this regard because organic carbon, and consequently mercury, concentrations were largely unchanging within the top 20 cm. This pattern matches the observations of Sulman et al [36], who found that tilled soils often have lower surface carbon concentrations, but deeper and more gradually decreasing concentrations, than non-tilled soils. Physical and elemental soil profiles by secondary landcover type for the conterminous United States. Soil organic carbon, bulk density, and fine fraction were generated using USDA soil data, retrieved from gNATSGO and PSU SIFEMEM. Mercury data were obtained from Smith et al [17] and were modeled based on latitude, longitude, and soil organic carbon. Landcover type was calculated using data from the NLCD 2011 dataset. All developed landcover types were combined for easier visualization. These data summarize the whole of the top 1 m of soil of the conterminous United States and are therefore dominated by mineral soil. Data is presented by soil layer, which is variable in depth.

Generation of pools for mercury and SOC
Using the modeled values, mercury pools were generated for all 11 layers representing the top 1 m of soil. Per unit depth, mercury mass is highest in surface soils (median ± SD = 11.8 ± 0.2 Gg in the top 5 cm) and decreases steadily with depth before a marked decrease from 80 to 100 cm (supplementary material figure S6). Cumulatively the top 1 m of soil contains 158.0 ± 1.8 Gg (median ± SD) of mercury in CONUS ( figure 3). The mean areal density of mercury is 1.5 ± 0.02 mg m −2 in the upper 5 cm and 0.09 ± 0.002 mg m −2 in the lowest 10 cm, with an overall mean areal density of 20.3 ± 0.2 mg m −2 for the top 1 m. Generally, the eastern half of the United States is elevated in mercury mass compared to the western half, with the exception of the west coast, which includes some of the highest mercury values calculated. Many areas of higher mercury concentration (see figure 1) are unintuitively lower in mercury mass than their surrounding areas (e.g. Everglades, Appalachian Mountains, Ozark Mountains); this contradiction is likely a result of the physical soil factors (e.g. low bulk density, high rock fraction) and are further explored below.
The concentration and mass of SOC in soil layers show similar patterns with depth as soil mercury (figures 2 and 3; supplementary material figures S3, S4 and S6). Both concentration and mass decrease with depth, though below the top 5 cm SOC content decreases more sharply than mercury. Mean (±SD) SOC concentrations range from 4.18% (±9.11) in the top 5 cm to 0.77% (±3.88) in the lowest 10 cm. SOC mass ranges from 9.8 Pg in the top 5 cm (1.3 kg m −2 ) to 0.2 Pg in the bottom 10 cm (0.03 kg m −2 ). Summed, the total mass of SOC in the top 1 m of soil of CONUS is 67 Pg. The areal density of SOC in the top meter is 8.63 kg m −2 , which agrees with the estimate of 8.56 kg m −2 reported by Bliss et al [18].

Mercury sourcing
The provenance of mercury in soils was investigated using bedrock Hg/Ti ratios from Richardson et al [29] and indicate that 62%-95% (10th and 90th percentile, respectively) of soil mercury is from nonparental sources. These sources may include natural and anthropogenic inputs as well as soil translocation, though the current methodology does not allow for further inquiry. The fraction of excess mercury is highest in surface soils and decreases throughout the profile, which is consistent with surface inputs being an important driver of excess mercury. While these estimates provide insight into the provenance of CONUS soil mercury, more data on bedrock mercury concentrations are needed to better constrain future estimates.

Patterns in mercury pools by landcover
Mercury pools vary by landcover type, though the patterns are markedly different than concentration due to variations in the physical characteristics of the soil (figure 2). Bulk density is negatively correlated with organic carbon concentrations, and consequently the high mercury input associated with leaf litter is often associated with low bulk density. This consideration results in a less pronounced enrichment of mercury in carbon-rich landcover types than mercury concentration alone might suggest. The rock fraction of soils also has a notable impact on pool calculations. For example, coniferous forests, which often reside at higher elevation, have a much lower fine fraction than the other forest types, resulting in smaller mercury areal densities compared to landcovers with lower mercury concentrations (e.g. developed). Conversely, croplands have the lowest rock fraction due to historical tilling and therefore have a greater mercury areal density than some landcovers with higher mercury concentrations (e.g. coniferous forests). Overall, mercury areal densities are highest in mixed and deciduous forests, pastures and croplands, and developed land types, while grasslands, shrublands, and barren landcover types are the lowest.

CONUS mercury pools in the global context
The present study helps to clarify the magnitude and uncertainty surrounding current global soil mercury reservoir estimates [10,15,16] as it is one of only a handful of studies that has calculated regional pools based on a large number of measurements collected from multiple depths using a randomized sampling design. Numerous studies estimate that global soils contain ∼1 Tg of mercury [1,9,37] or less [6,7]. Based on these estimates, the 158 Gg estimated in the present study would result in CONUS containing greater than 15% of the global soil mercury reservoir while representing only ∼5% of Earth's land area. This difference largely reflects distinctions in how soil pools are delineated. Many studies focus on more rapidly cycling, organic-rich surface pools in upper horizons (10-30 cm), as these are most relevant to human and ecosystem health. In contrast, studies that include deeper and less dynamic soil pools, which as this work demonstrates can be sizable, produce much larger soil reservoir estimates. Considered alongside other 0-100 cm regional estimates by Schuster et al (755 Gg) [10] and Olson et al (408 Gg) [8] the top meter of global soils pools may contain substantially more than the estimated 1 Tg Hg contained in rapidly-cycling surface pools. The approach adopted in the present study allows for pool calculation at variable depth. Methods for distinguishing between organic and non-organic soil mercury fractions would be helpful for mass budgets and pool estimates to increase consistency across studies.
Ballabio et al [16] provides a valuable comparison for the results presented here, as the study used a different methodology to estimate soil mercury stocks with a large dataset over an area with comparable climate and land cover. Ballabio et al calculated that the top 20 cm EU soils contain approximately 45 Gg of mercury. This estimate is similar to the top 20 cm estimate presented here for CONUS (44 Gg), despite the former being roughly half the area (∼4.4 million km 2 vs 7.6 million km 2 ). Two factors may contribute to the difference in areal density of soil mercury between the US and EU study areas. First, Ballabio et al did not account for the rock fraction when calculating the soil mercury pool, which may elevate their estimates in areas with substantial rock content. If the rock fraction were not considered for CONUS, the estimated soil mercury pools would be 50.9 Gg. Second, while the first quartile and median are similar between the two studies (0.0135 vs 0.0136 µg g −1 and 0.0233 vs 0.0240 µg g −1 , for the EU and US, respectively), the mean and 3rd quartile are higher in EU soils (0.0554 vs 0.0341 µg g −1 and 0.0492 vs 0.0423 µg g −1 , respectively) suggesting that enriched areas in the EU may have greater mercury contamination than those in CONUS. This potentially reflects the different emission and industrial histories and the greater population density (∼1.5×) of the EU or differences in soil characteristics between the two survey areas. Overall, the two estimates are relatively similar despite the different approaches employed.
Estimates of soil mercury pools now exist for much of Europe, the United States, and the Arctic. Although there are still many regions with limited soil mercury data, each new soil mercury pool estimate generated using largescale, field-collected datasets helps to improve current estimates of the global soil mercury reservoir. Moving forward, soil survey information from underrepresented areas (e.g. the global south) will greatly improve global estimates.

Limitations
The mercury reservoir estimates presented here have several limitations that should be noted. Modeling and pool calculation relied on available soil coverages, namely SOC, bulk density, and rock fragment volume. These data were all sourced from USDA soil databases (e.g. gNATSGO, STATSGO) and have their own limitations, which propagate to the final soil mercury pool estimates. Soil sampling efforts are biased towards mineral soil estimates and away from some soil types (e.g. wetlands are not well represented). As a result, the mercury soil pool estimates presented here are biased towards the mineral soil and do not represent organic rich landcover types. Additionally, the scale of bulk density and SOC measurements used here was relatively coarse (9 km 2 ) and likely does not capture the heterogeneity of these metrics. Consequently, the reservoir estimates presented here should be used at the regional to continental level. Despite these limitations, this work presents the most robust estimate of the soil mercury reservoir for CONUS currently available.

Data availability statement
All data used in this study are free and publicly available. The data that support the findings of this study are openly available at the following URL/DOI: https://pubs.usgs.gov/ds/801/. Precipitation data are available through the GEOS-5 FP meteorology product from the NASA Global Modeling and Assimilation Office (https://gmao.gsfc.nasa.gov/GEOS_ systems/geos5_access.php). Bulk density and soil organic carbon data are available through the USDA gNATSGO database (www.nrcs.usda.gov/wps/portal/ nrcs/detail/soils/survey/geo/?cid=nrcseprd1464625). Rock fraction volume data are available through the Soil Information for Environmental Modeling and Ecosystem Management (www.soilinfo. psu.edu/index.cgi?soil_data&conus). Elevation data are available through the USGS (www.usgs.gov/ core-science-systems/ngp/tnm-delivery/gis-datadownload). Landcover and Land use data are available through the USGS National Land Cover Database (www.usgs.gov/centers/eros/science/nationalland-cover-database?qt-science_center_objects=0# qt-science_center_objects). Ecoregion shapefiles are available through the USEPA (www.epa.gov/ecoresearch/ecoregions).