Mixed trends in heavy metal-enriched fugitive dust on National Park Service lands along the Red Dog Mine haul road, Alaska, 2006–2017

This study presents the status and trends of long-term monitoring of the elemental concentrations of zinc (Zn), lead (Pb), and cadmium (Cd) in Hylocomium splendens moss tissue in Cape Krusenstern National Monument (CAKR), Alaska, adjacent to the Red Dog Mine haul road. Spatial patterns of the deposition of these metals were re-assessed for the period from 2006–2017 following an identical study that assessed trends between 2001–2006. In contrast to the widespread and steep declines in Zn and Pb levels throughout most of the study area between 2001–2006, this study showed more mixed results for 2006–2017. At distances within 100 m of the haul road, only Pb decreased between 2006–2017. At distances between 100–5,000 m, however, both Zn and Cd decreased between 2006–2017, with high probabilities of decrease and percent decreases of 11–20% and 46–52% respectively. Lead did not decrease in any of the more distant areas. Following earlier work on lichen species richness in the study area, it appears that 2017 Zn levels are approaching those associated with “background” lichen species richness throughout a relatively large proportion of the study area at least 2,000 m from the haul road and several km from the port site. The findings in this study may be used to plan additional mitigation measures to reduce Zn deposition related to impacts on lichen communities.


Introduction
Cape Krusenstern National Monument (CAKR) was created in 1980 [1] as a US National Park Service (NPS) unit.Located above the Arctic Circle north of Kotzebue, Alaska (Fig 1) the monument features low windswept mountains along the Chukchi Sea coast, a series of lagoons that open seasonally to the sea and feature a huge variety of fish and bird habitats, extensive softsediment shorelines, and a wide coastal plain [2].In 1985, Congress authorized an industrial easement through the monument allowing for the trucking of metal ore concentrates from the Red Dog Mine to the Red Dog Port [3].Both the mine and the port are located on Northwest Alaska Native Association (NANA) lands, and are linked via the Delong Mountain Transportation System (DMTS) haul road (referred to hereafter as "haul road").The haul road crosses 32 km of NPS lands.Since 1989, the mine has transported Zn and Pb concentrates (approximately 55% in fine powder form) year-round from the mine to large concentrate storage buildings at the port.In doing so, fugitive dusts enriched with Zn, Pb, and Cd have been dispersed onto NPS lands both inside and outside of the industrial easement [4][5][6].
The moss Hylocomium splendens (Hedw.)Schimp has been used as a biomonitor of heavy metals, atmospheric pollutants, and even anthropogenic substances for over 30 years and mosses in general have been used for this purpose since the 1960s [7][8][9].Dozens of studies have utilized the low cost and high repeatability of moss biomonitoring as an alternative to instrumented monitoring (see [7] for review), and this approach has been effective for both temporal and spatial analysis of air pollution deposition [10,11].While comparisons of elemental concentrations in moss tissue across and between broad regions (e.g., northern Europe, North America) can be partly confounded by environmental and depositional particularities, moss biomonitoring is widely known for high accuracy and consistency on local to smaller regional scales [12].
There is a large volume of literature on biomonitoring of heavy metal pollution related to metal smelting (e.g., [13][14][15][16][17]), but studies on biomonitoring of fugitive dust from mining operations are far fewer.While smelters are often responsible for acute effects on vegetation in a localized area, fugitive dusts can cause chronic effects on a landscape level for hundreds of km 2  or more [18].When moss biomonitoring for contaminant deposition is paired initially with monitoring of biological effects, simple moss biomonitoring may then be used on a regional level to predict biological consequences.
NPS first measured the elemental concentrations of Zn, Pb, and Cd in the moss Hylocomium splendens in 2000 [6] along six transects from the DMTS haul road.Sampling was expanded to a broader area of inference using a spatial model of elemental deposition patterns in 2001 [5].In 2006, NPS simultaneously sampled lichen communities and remeasured elemental concentrations in a new sampling grid using five strata based on distance from the haul road [4,18].The 2001 and 2006 analyses were formally included as a protocol of the NPS Arctic Network Inventory and Monitoring Program, ensuring decadal monitoring over the long term [19].This protocol strives to determine the status and trends of key ecological indicators using a stable sampling design, methodology, and analysis.In 2017, NPS conducted the third measurement of moss elemental concentration in CAKR, co-located with the second sampling of lichen communities.
This paper presents the results of a decadal spatial analysis of contaminant deposition based on the resampling of locations originally sampled in 2006 [4] and fitting the 2006 and 2017 data to the model used previously for the 2001 and 2006 data.Our objectives were to: 1) describe the results of the spatial analysis of 2017 data using a Bayesian geostatistical model of elemental concentrations of Zn, Pb, and Cd, 2) summarize strata-based predictions of elemental concentrations, 3) summarize and map pointwise posterior predictions for elemental concentrations of these 3 metals, and 4) assess trends in elemental concentrations between 2001-2017 by determining the probability of decrease by element for both points and strata.
This article describes the results of a long-term monitoring protocol, and as such, it represents the latest findings in a series of research articles using stable methodology and analyses.For brevity, we recommend readers refer to earlier articles [4,5,18] for greater detail on study design, methods, lab analyses, and statistics.

Ethics statement
This study was approved by the National Park Service-Western Arctic National Parklands interdisciplinary project review team under the provisions of the National Environmental Policy Act [20].The project was assessed to be low in impacts to natural resources, cultural resources and subsistence and was granted a categorical exclusion from further review.

Study area, field and lab methods
Methods for this long-term monitoring protocol are thoroughly described in a previous openaccess publication (Neitlich et al. [4]) and in Swanson and Neitlich [19] including study area, study design, field sampling, and lab analyses.In brief, the goal of this study was to detect trends in Zn, Pb, and Cd concentrations in the feather moss Hylocomium splendens in CAKR adjacent to the DMTS haul road between 2006 and 2017.To obtain a landscape-scale assessment of contaminants, samples of the H. splendens were obtained in 2017 immediately adjacent to the 104 long-term vegetation monitoring plots established in 2006 [4,18] and from 6 other roadside collection points (Fig 1).These 104 vegetation monitoring plots formed the network from which �95% of the samples in both 2006 and 2017 were obtained.
Earlier studies along the DMTS haul road showed that heavy metal deposition declined in a roughly logarithmic curve from the road [5,6].Sampling distances from the road were thus chosen in 2006 to obtain an even spread of values along this deposition gradient.Twelve non-linear "transects" were created according to the distances in Table 1, with points at 10, 50, 100, 300, 1000, 2000, and 4000 m.To help account for spatial autocorrelation in the model, one additional plot (termed an "autocorrelation plot") was established in each transect at a random distance 10 to 20 m from the 1,000, 2,000, or 4,000 m plot in each transect.Each "transect" consisted of a group of points within one of the two most common land cover classes in the study area at specified distance classes.Points were chosen at random in GIS in close geographic proximity to a line, but did not represent a strictly linear feature.Final locations were additionally randomized using a random number generator to select a group of at least 10 contiguous pixels of the target landcover class closest to the transect line at a particular distance.Linear transects were tested in the planning stages, but rejected because it was not possible to find both the correct land cover classes and distances along a straight line.The two most common landcover types accounting for more than 66% of the study area were Upland Moist Dwarf Birch-Ericaceous Shrub and Upland Moist Dwarf Birch-Tussock Shrub [21].These closely-related classes were each dominated by cottongrass (Eriophorum vaginatum), dwarf bich (Betula nana) and mixed Ericaceous shrubs (e.g., Ledum decumbens, Vaccinium uliginosum, Vaccinium vitis-idaea).
Ten reference samples were obtained in southern CAKR at distances of 42,000 to 64,000 m from the DMTS haul road adjacent to the long-term vegetation plots there.To avoid the presence of high dust and mud immediately adjacent to the road, long term vegetation monitoring plots transects started at a distance of 10 m from the road.To capture the higher contaminant values immediately adjacent to the road, therefore, a small number of roadside moss samples unassociated with the long-term monitoring plots were also obtained in both 2006 and 2017.A slightly larger number of samples (131) were collected in 2017 (Table 1) than in 2006 (125).This was due to a greater number of field duplicate samples and two additional stand-alone roadside moss samples.Moss samples were either collected and stored dry in heavy duty chemical sample bags or collected wet and frozen until cleaning and analysis, and were collected and stored double bagged using sterile techniques [4].The concentrations reported from lab analysis are mg/kg of analyte in Hylocomium splendens moss tissue, dry weight.

Posterior predictions modeling
We analyzed the 2006 and 2017 chemical concentration data using a geostatistical model that accounts for spatial correlation among observations.We let y ij denote the observed chemical concentration at spatial location i = 1,. ..,n for measurement j = 1,. ..,J i and x i denote a vector of predictor variables (including 1 as the first element corresponding to an intercept parameter) measured at location i.Our analysis included distance to road (natural logarithmic scale, referred to as "log"), an indicator for south side of the road, and their interaction as predictor variables.We formulated our geostatistical model as where β is a vector of regression coefficients.The α i terms represent site-level random effects that account for the repeated measurements at sampled locations.We assumed these terms were independent across sites and were normally distributed with mean zero and variance s 2 1 .To account for spatial correlation among observations, we also included a spatial random effect η i for each location.Defining η�(η 1 , η 2 ,. ..,η n ) 0 , we assumed η~N(0, S) where the covariance among spatial effects was specified as a function of the distance between the corresponding locations.Specifically, we used the exponential covariance function so that where s 2 2 and θ are parameters to be estimated and d i ĩ is the distance between locations i and ĩ.The measurement errors � ij were assumed to be independent and normally distributed with mean zero and variance s 2 0 .Separate geostatistical models were fit for each combination of element (Cd, Pb, and Zn) and year (2006 and 2017).
We used Bayesian methods to fit these models and specified weakly informative prior distributions for each parameter.The prior distribution for β was specified as multivariate normal with mean equal to a zero vector and covariance matrix 100I.For s 2 0 ; s 2 1 , and s 2 2 we assumed Gamma(2,2) (shape and rate parameterization) priors.We specified a Half-normal(0,5000) prior distribution for the spatial range parameter θ based on the scale of the observed distances between sampled locations.Note that these prior distributions differed slightly from those assumed by Neitlich et al. [4].However, for the 2006 data, inferences were not sensitive to the change in prior distributions based on comparisons to previous results (e.g., see tables in Results section).All statistical analyses were conducted in R [22].Each model was fit with Markov chain Monte Carlo (MCMC) using Stan [23] called from R with the rstan package [24].Our Bayesian geostatistical models were fit using 4 chains of 4,000 iterations each where the first half (2,000 iterations) of each chain was discarded as adaptation/burn-in and the second half was used for inference.
Posterior predictions of chemical concentrations were made across a grid of 2,374 total locations for 200 random posterior draws (to reduce computation time).The prediction grid was established by Neitlich et al. [4] and consists of five strata that were defined based on proximity of points to the road.For posterior draw k, let ŷðkÞ year;l denote the predicted concentration (on the original scale) for a particular year and prediction grid point l.We created maps summarizing the posterior predictive distribution at each grid point based on the posterior mean, 2.5% posterior quantile, and 97.5% posterior quantile.Additional maps were created to summarize the posterior distribution of percent change and the posterior probability of decrease at each prediction grid point.Percent change associated with posterior draw k and point l was calculated as

: ð3Þ
The posterior probability of decrease was calculated based on the proportion of posterior draws where ŷðkÞ 2006;l was larger than ŷðkÞ 2017;l .Additional summaries of the posterior predictive distributions were made by strata.Following Neitlich et al. [4], we defined the average concentration for strata t and posterior draw k as where m t is the total number of prediction grid points in strata t.We then summarized the posterior predictive distributions for the change in average concentration, the percent change in average concentration, and probability of decrease in average concentration.

Spatial model
Summaries of the posterior distributions for the regression coefficients are shown in Table 2.
Based on these results, we can infer that log-concentrations decrease as log-distance to the haul road increases for all three elements and both years because the credible intervals for β 1 include only negative values.For Cd, the increase in the intercept parameter from 2006-2017 suggests that log-concentrations increased at locations close to the road.However, the posterior mean of the road distance coefficient decreased from 2006-2017, indicating that Cd log concentrations decreased more rapidly as distance from road increases.The posterior means of regression coefficients for Pb in 2017 show that log concentrations are estimated to be lower closer to the road but decrease less rapidly as distance from road increased relative to 2017.
The posterior distributions of the regression coefficients for Zn are similar for 2006 and 2017.
To better explore the differences in estimated log concentrations across years, we summarize the posterior predictive distributions for these models in Table 2.

Spatial patterns of deposition
Strata-based summarization.The strata used in this study followed Table 3 and Fig 3 in Neitlich et al. [4] and are based on distance from the DMTS haul road as follows: Stratum 1, 0-100 m; Stratum 2, 100-2,000 m; Stratum 3, 2,000-4,000 m; Stratum 4, 4,000-5,000 m, Stratum 5, 30,000-64,000 m.Each stratum represents a distance on both the north and south side of the road with the exception of Stratum 5, which occurs only on the south side of the road.Strata-based summarization is presented for Zn, Pb, and Cd independently according to side of road as well as grouped.
In Stratum 1, change in Zn concentrations in Hylocomium splendens moss between 2006-2017 was highly uncertain with a 56% probability of decrease (Table 3).In Strata 2, 3, and 4, the probabilities of a decrease in Zn between 2006-2017 were >90% with mean decreases of 12-15 mg/kg (representing a percent decrease between 11-20%).In Stratum 5, there was a 78% probability of a decrease from 40 to 31 mg/kg.The main difference between the 2001-2006 period and that of 2006-2017 is that Zn concentrations did not decrease in 2017 in Stratum 1 as they did for every other stratum in 2017 (Fig 2).
Changes in Pb levels between 2006-2017 were considerably different from that of Zn (Table 4).There was 100% probability of a decrease of 80 mg/kg (from 198 to 118 mg/kg) in Stratum 1.Other strata had only moderate probabilities of decrease (29-82%) with very small change quantities (0-1.7 mg/kg).These findings do not point to any meaningful change in Pb levels over this time period at distances >100 m, but a very sizeable decrease close to the The change in Cd during this time period was yet again different from that of Zn and Pb (Table 5).In Stratum 1, there was an 80% probability of a decrease of 0.60 mg/kg from 5.64 mg/kg to 4.64 mg/kg.In all other strata, there was a 97-100% probability of decrease, with percent decreases ranging from 41-51%.This level of change represents a much more sizeable decrease than was recorded between 2001-2006 [4], during which meaningful change occurred in Stratum 1 but was either small in quantity or a possible increase in more distant strata (Fig 2).
Point-based mapping.Point-based mapping provides an opportunity to examine spatial patterns in finer detail than strata-based mapping.In contrast to the widespread and substantial decrease in concentrations of Zn, Pb, and Cd in moss tissue from 2001-2006, the period from 2006-2017 was more variable.For example, while a decrease in Zn of �25% was predominant in the DMTS corridor from 2001-2006 (Fig 8 in [4]), the predominant changes in Zn ranged from +50% to -50% between 2006-2017 (Figs 3-5).Similarly, the posterior probability of decrease was much higher over the corridor as a whole in the earlier monitoring period.Between 2001-2006, a dominant portion of the DMTS corridor out to 2,000 m showed �75% probability of decrease (Fig 11 in [4]), compared to a patchy mosaic of 40-80% probability of decrease between 2006-2017 (Fig 6).On a landscape level, the largest change in Zn concentrations between 2006-2017 was the decrease from the 60-100 mg/kg class to the 13-60 mg/kg class over a large portion of the north side of the road (Figs 4 and 5).This decrease is significant as the provisional threshold for effects on lichen species richness is approximately 58-70 mg/kg of Zn in Hylocomium moss tissue, meaning this area may once again fall below levels that cause lichen mortality [18].Less clear is whether recovery can occur at these levels given both legacy Zn in soils [25] and ongoing lower levels of Zn deposition.The areas closest to the haul road (�100 m) showed the persistence of Zn concentrations in the 400-1200 mg/kg range and a predominance of increase of up to 25%.The widespread probabilities of decrease between 40-80% pointed to a lack of clear trends in this zone of most elevated concentrations.This is corroborated by the 56% probability of increase in Stratum 1 (Table 3).The concentrations of Pb in moss tissue between 2006-2017 changed in a bimodal pattern.The area immediately adjacent to the road decreased predominantly by �-25%, while more distant areas north of the road tended toward increases of �50%.Large portions south of the road increased by 25-100%.(Figs 7-9).The probabilities of decrease (Fig 10 ) were strongest immediately adjacent to the DMTS haul road (predominantly 80-100%), mixed on the north side of the road (20-80%) and predominantly low on the south side of the road (20-40%, i.e., 60-80% chance of increase).The high probability of decreases adjacent to the road (predominated by the 0-25 and 25-50% classes) mirrors the 100 percent probability of a decrease of 39% (or approximately 80 mg/kg) in Table 4.The probability of decrease within the 1,000 m haul road buffer zone was predominantly in the 40-60% range with low overall percent change values (0-25%), indicating little change.Pb concentrations stayed at �10 mg/kg in moss tissue in this zone and commonly in the 200-3,200 mg/kg range adjacent to the haul road.There appeared to be a 60-80 percent probability of decrease in CAKR adjacent to (and especially on) the north side of the Port Site, with decreases of up to 25%.The high dominance of the 20-40% probability of decrease (i.e., 60-80% probability of increase) throughout the south side of the corridor at distances >1,000 m from the road is suggestive of an increase, though certainly not definitive.Table 4 provides insight that regardless of the probabilities, the mean change in these more distant southern areas is quite small (<2 mg/kg Pb in moss tissue).
As also seen in Fig 2 and Table 5, change in Cd was stronger and more certain at distances beyond the 100 m of the immediate DMTS haul road corridor .Almost the entire study area outside of the immediate corridor decreased by 25-75% with about equal amounts in the 25-50% and 50-75% classes.The probabilities of decrease were predominantly in the >95% and the 80-95% classes.Along the road, trends were not as clear, as there was at the  [14,18,26,27].This study documents modest levels of decrease (10-15 mg/kg in moss tissue) of Zn in areas beyond 0-100 m between 2006-2017.Neitlich et al. [18] reported a preliminary Zn threshold of 60-70 mg/kg Zn in moss tissue for effects on lichen species richness (LSR) in CAKR.In that study, Zn deposition at levels higher than this resulted in a 50-100% reduction of LSR on 11 km 2 and 25-100% reduction of LSR on 55 km 2 .A total of 157 km 2 showed reductions in LSR of 5-100%.Di Meglio [28] reported that lichen species richness continued to decline in 2017 despite the two previous periods of decrease in Zn concentrations.
While the effects of Zn on depressing LSR and causing injury and mortality to lichens is well-known, we know much less about how lichens recover following pollution reduction.Several studies [29][30][31][32][33] report recolonization of impacted areas following reductions in sulfur dioxide or copper from smelters, however there are no direct analogs of the current situation in the literature.Most studies have used epiphytic lichens to study recolonization, thus have not had to consider the potential effects of residual contaminants in soil pore water on terrestrial lichens such as occur in CAKR.Melby [25] reported considerable pools of Zn, Pb, and Cd in the organic horizon pore water in the study area, and it is likely that this pore water intermittently comes in contact with lichens reaching deeper into the organic horizon (e.g., Cladonia stygia, Cetraria laevigata).Most recolonization studies concluded that reasonable recovery took decades after pollution stopped due to slow lichen growth and limitations on nearby Table 5. Modeled concentrations (mg/kg dry weight) of cadmium (Cd) in Hylocomium splendens moss tissue in 2006 and 2017, average change, percent change, and probability of decrease for 5 strata.Results were tabulated for the whole data set as well as separately for the north and south sides of the road.Results are presented for the posterior mean, standard deviation (SD), and the lower and upper bounds representing the 2.5% and 97.5% endpoints of the 95% credible interval.Pr.Dec. presents the posterior probability of decrease in concentration between 2006 and 2017.propagule sources [32,33].Lichens and bryophytes along the DMTS haul road face reduced but ongoing pollution from both deposition and pore water exposure, which could conceivably slow recovery additionally.The Zn concentrations in 2006 that caused the decline in lichens and bryophytes along the DMTS haul road averaged between 77-566 mg at distances from 10m to 4000 m from the road.Between 2006-2017, the Zn decreases in strata for which there was fairly high probability of decrease were �15mg/kg, suggesting that the bulk of the pollution signal responsible for causing the effects reported by Neitlich et al. [18] were ongoing in 2017.Still, a reduction of up to 15 mg/kg Zn could place the cleaner parts of the study area, those with higher LSR, near or below the preliminary effects threshold of 60-70 mg/kg.Other biological effects including reduction in lichen cover, reduction in bryophyte cover, reduction in lichen heights, reduction in bryophyte frond width have been reported from the study area [18,27,34], but thresholds for these effects have not yet been determined.After thresholds for Zn by itself or in concert with Pb and Cd are defined for various biological effects, the results of this study may be used to model these effects spatially at the landscape level in CAKR.Moreover, because moss biomonitoring is relatively cheap and efficient relative to detailed biological field work, thresholds analysis may permit effective future monitoring simply via elemental analysis and spatial modeling of the elemental concentrations of moss samples.

Cd
There have been hundreds of studies in the past 5 decades using both lichens and bryophytes as bioindicators of air pollutants including N, S, semi-volatile organic compounds, and heavy metals [7,35,36].Like this study, the majority have focused on quantifying levels of pollutants in lichen or moss tissue, and to a lesser extent, mapping deposition patterns [10].The National Park Service monitoring project underlying the current research [19] is unique in examining both spatial patterns of heavy metals deposition on a sub-regional scale and linkages between pollution levels and mortality of lichens and bryophytes.In the context of the current literature, research on contaminants from mining-related fugitive dusts and their biological effects are uncommon compared to biomonitoring of urban and industrial point sources, smelters, regional air pollution, and soils polluted by mining wastes.The research comprising nearly 20 years of monitoring and modeling of fugitive dusts from the DMTS haul road in CAKR [4-6, 18, 27, 37, 38] represents an intensive case study on this distinct issue.Mitigation.Monitoring between 2001-2006 showed that the Red Dog Mine's mitigation during this period helped to reduce levels of Zn, Pb, and Cd in fugitive dusts significantly [4][5][6].Major mitigations in during this period included replacing tarp-covered concentrate haul trucks with a fleet of trucks using hydraulically-covered lids, implementing a truck rinse station for use in summer months, containing escapement from the concentrate storage buildings via sealed loading and unloading facilities, enclosing the conveyor system that moves concentrates from the concentrate storage buildings to vessels at the Port Site, applying dust palliatives to the roadbed, and segregating road traffic to reduce concentrate tracking by vehicles [27].The period between 2006-2017 showed improvements in some areas and slight worsening or  [18,27,34], additional mitigations appear to be necessary.Study of surface peat soils along the DMTS haul road showed elevated levels of Zn, Pb, and Cd between 10-100 m from the road, with means of 484-785 mg/kg of Zn in organic soils [25].Additional mitigation would slow the accumulation of heavy metals in these upper organic soil horizons.study area close to (or below) the preliminary impact threshold (approximately 60 to 70 mg/kg of Zn) for lichen species richness.More work on the contaminant thresholds associated with different biological parameters (e.g., lichen species richness, lichen cover, nonvascular plant height and vigor) is needed.Lichen recovery is conceivable in areas that have dropped below the preliminary Zn threshold.However, the levels of Zn at which recovery would begin and the timeline for recovery will require ongoing field study of lichen communities and contaminant levels.The next National Park Service remeasurement of both vegetation and contaminants in moss tissue is scheduled for 2027.

Fig 2 .
Fig 2. Changes in mean modeled Zn, Pb, and Cd concentrations in Hylocomium splendens between 2001 and 2017.Posterior means that showed >90% probability of decrease from one sampling period to the next are marked with a red asterisk.Error bars show the high and low bounds of the 95% credible interval.The 2001-2006 change data are from Neitlich et al. [4] using exactly the same modeling procedures as the 2006-2017 analysis.https://doi.org/10.1371/journal.pone.0297777.g002

Fig 3 .
Fig 3. Modeled 2006 Zn concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g003

Fig 4 .
Fig 4. Modeled 2017 Zn concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g004

Fig 5 .
Fig 5. Percent change in Zn concentrations in moss tissue in CAKR along the DMTS haul road, 2006-2017.Percent change between 2006 and 2017 Zn concentrations in moss tissue along the DMTS.The 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations are shown in smaller maps.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g005

Fig 6 .
Fig 6.Probability of Zn decrease in moss tissue in CAKR along the DMTS haul road, 2006-2017.Posterior probability of decrease between 2006 and 2017 of Zn concentrations in moss tissue at each prediction point.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g006

Fig 7 .
Fig 7. Modeled 2006 Pb concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g007

Fig 8 .
Fig 8. Modeled 2017 Pb concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g008

Fig 9 .
Fig 9. Percent change in Pb concentrations in moss tissue in CAKR along the DMTS haul road, 2006-2017.Percent change between 2006 and 2017 Pb concentrations in moss tissue along the DMTS.The 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations are shown in smaller maps.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g009

Fig 10 .
Fig 10.Probability of Pb decrease in moss tissue in CAKR along the DMTS haul road, 2006-2017.Posterior probability of decrease between 2006 and 2017 of Pb concentrations in moss tissue at each prediction point.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g010

Fig 11 .
Fig 11.Modeled 2006 Cd concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g011

Fig 12 .Fig 13 .
Fig 12. Modeled 2006 Cd concentrations in moss tissue in CAKR along the DMTS haul road.Smaller maps show the 2.5 th and 97.5 th percentiles (lower and upper bound of the 95% credible interval) of the modeled concentrations.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g012

Fig 14 .
Fig 14.Probability of Cd decrease in moss tissue in CAKR along the DMTS haul road, 2006-2017.Posterior probability of decrease between 2006 and 2017 of Cd concentrations in moss tissue at each prediction point.DMTS haul road and distance-based buffers at 1 km, 2 km, and 4 km are shown along with NANA lands and the CAKR boundary.https://doi.org/10.1371/journal.pone.0297777.g014

Table 1 . Number of plots and tissue samples by distance from the DMTS haul road, 2017. Standard
and autocorrelation plots sampled both vegetation and moss tissue.Autocorrelation plots were established to provide autocorrelation estimation in the posterior predictions model.A

Table 2 . Summaries of posterior distributions for regression coefficients in Eq 1.
The means and standard error of the posterior distributions are provided along with the lower 2.5% and upper 97.5% endpoints of the 95% credible interval.

Table 3 . Modeled concentrations (mg/kg dry weight) of zinc (Zn) in Hylocomium splendens moss tissue in 2006 and 2017, average change, percent change, and probability of decrease for 5 strata.
Results were tabulated for the whole data set as well as separately for the north and south sides of the road.Results are presented for the posterior mean, standard deviation (SD), and the lower and upper bounds representing the 2.5% and 97.5% endpoints of the 95% credible interval.Pr.Dec. presents the posterior probability of decrease in concentration between 2006 and 2017. https://doi.org/10.1371/journal.pone.0297777.t003

Table 4 . Modeled concentrations (mg/kg dry weight) of lead (Pb) in Hylocomium splendens moss tissue in 2006 and 2017, average change, percent change, and probability of decrease for 5 strata.
Results were tabulated for the whole data set as well as separately for the north and south sides of the road.Results are presented for the posterior mean, standard deviation (SD), and the lower and upper bounds representing the 2.5% and 97.5% endpoints of the 95% credible interval.Pr.Dec. presents the posterior probability of decrease in concentration between 2006 and 2017.
https://doi.org/10.1371/journal.pone.0297777.t004same time a low probability of decrease (40-80%) and percentages of increase between 0-100%.It is unclear why Cd decreased markedly during the 2006-2017 period at outer distances but only decreased in Stratum 1 between 2001-2006.Biological effects.A large volume of literature describes the toxic effects of Zn deposition on lichens and bryophytes