Multidecadal declines in particulate mercury and sediment export from Russian rivers in the pan-Arctic basin.

SignificanceRussian rivers are the predominant source of riverine mercury to the Arctic Ocean, where methylmercury biomagnifies to high levels in food webs. Pollution controls are thought to have decreased late-20th-century mercury loading to Arctic watersheds, but there are no published long-term observations on mercury in Russian rivers. Here, we present a unique hydrochemistry dataset to determine trends in Russian river particulate mercury concentrations and fluxes in recent decades. Using hydrologic and mercury deposition modeling together with multivariate time series analysis, we determine that 70 to 90% declines in particulate mercury fluxes were driven by pollution reductions and sedimentation in reservoirs. Results suggest that Russian rivers likely dominated over all other sources of mercury to the Arctic Ocean until recently.

High levels of methylmercury accumulation in marine biota are a concern throughout the Arctic (1), where coastal ocean ecosystems receive large riverine inputs of mercury (Hg) and sediment. Recent estimates indicate that riverine export of Hg to the Arctic Ocean (36 to 44 Mg⋅y−1) (2–4) is comparable to inputs from coastal erosion (15 to 47 Mg⋅y−1) and nearly half that of atmospheric inputs from direct deposition and via sea ice and snow meltwater (45 to 98 Mg⋅y−1) (5, 6). Particulate Hg (PHg) accounts for roughly half of the Hg export in northern rivers (3), in part because of the strong binding affinity between Hg and organic matter in suspended solids (7). River discharge and total suspended solids (TSS) are thus strong predictors of Hg export to the coastal Arctic Ocean (4), where much of the riverine Hg is thought to be resequestered in marginal sediments or reduced and evaded to the atmosphere (8). Superimposed on these controls on Hg cycling, climate warming and intensifying hydrologic cycles are strengthening land–freshwater linkages at northern high latitudes (9). These effects are reflected by increases in riverine export of weathering ions and dissolved nutrients across the circumpolar north in recent decades (10–12) and regional intensification of Hg release into fluvial networks (13, 14). Yet, Hg concentrations in fish harvested from major Russian rivers, which account for 70% of riverine Hg inputs to the coastal Arctic Ocean (4), have declined since the late 20th century (15, 16). No published temporal data exist to evaluate trends in Russian river Hg against other hydrochemical constituents which have exhibited striking multidecadal increases.

Amid recent and ongoing environmental change, long-term trends in riverine Hg are likely to be controlled by factors that influence sediment and pollution loading to watersheds. First, natural and/or anthropogenic changes in river discharge (e.g., reservoirs) may alter fluvial energy and thus erosion of sediment Hg from river banks and beds (17). Dams for energy and water resources are widespread across Russia (18), and associated reservoirs are likely to reduce net sediment export by increasing sedimentation (19). Additionally, industrial activities are historically prominent local and regional sources of pollution across Russia and in 2010 released an estimated 10 Mg of Hg via effluent directly into freshwaters (17). Polluting activities in major Russian watersheds include Chlor-alkali–based sodium hydroxide production to bleach paper products (20), coal- and gas-fired heat and power generation (21), gold mining (22), and municipal waste processing (23). From the late 1980s to late 1990s, economic output in Russia declined by 44% as the former Union of Soviet Socialist Republics (USSR) collapsed (24). This decreased industrial activity may have reduced Hg emissions to the atmosphere and curbed direct releases into aquatic ecosystems. Finally, the onset of Hg emission controls in the mid-20th century (25) reduced atmospheric Hg deposition in the northern hemisphere (26). Atmospheric Hg deposited to soils may be mobilized into fluvial networks via leaching and runoff over decades to centuries (27). The scarcity of riverine Hg data available from watersheds in the Russian sector of the pan-Arctic basin hinders a more concrete understanding of controls on northern high-latitude Hg cycling in recent decades.

Here we present multidecadal observations on TSS and PHg concentrations for eight major Russian rivers which discharge into seas bordering the Arctic Ocean: the Onega, Northern Dvina, Mezen, Pechora, Ob’, Yenisey, Lena, and Kolyma (Fig. 1). Using 1,300 paired measurements of TSS and PHg concentration and daily river discharge made from the 1970s to the late 2010s as part of a long-term Russian biological monitoring program, we investigate spatiotemporal trends in TSS and Hg concentration and export with the aim of better understanding recent changes in the Arctic Hg cycle. We then assess potential controls on changes in riverine Hg using observations on TSS, inventories of reservoir establishment (18), records of industrial Hg release (17), climate reanalysis data (28), and estimates of historical atmospheric Hg deposition simulated using a global atmospheric chemical transport model for Hg species (GEOS-Chem) (29). Our findings help to conceptualize future trajectories of Hg cycling at northern high latitudes, where rapid environmental change is expected to accelerate Hg release into ecosystems (13, 14, 30).

Fig. 1.

The eight rivers examined in this study span the longitudinal breadth of Russia. Their watersheds encompass nearly 10,000,000 km2 and contain varied histories of industrial pollution and reservoir establishment associated with dams. The Yukon and Mackenzie Rivers are shown here and referenced in the main text but were not sampled as part of this study.

Results and Discussion

Multidecadal Declines in TSS and PHg Export.

Our multidecadal hydrochemical observations, together with flux estimates using weighted regressions of concentrations on time, discharge, and season (WRTDS) software (31, 32), reveal large declines in TSS and PHg fluxes in major Russian rivers in the pan-Arctic basin between the 1970s and 2010s (Table 1). Our trend analyses leverage WRTDS flow normalization, which controls for the influence of variability in discharge to better understand potential effects from disturbance (e.g., pollution and reservoirs) on constituent concentration and flux. Trend significance results for flow-normalized concentration and flux from WRTDS are equivalent to two-sided P values and we report them following ref. 32: highly likely (≥0.95 and ≤1.0), very likely (≥0.90 and <0.95), likely (≥0.66 and <0.90), about as likely as not (>0.33 and <0.66), and unlikely (≤0.33). Declines in flow-normalized TSS flux were highly likely in the Ob’ (19.2 Tg⋅y−1), Lena (17.4), Yenisey (5.6), Pechora (4.3), and Mezen (0.3) and likely in the Northern Dvina (0.9). For these rivers, relative change in flow-normalized TSS flux ranged from –23% (Northern Dvina) to –81% (Pechora), whereas increased flow-normalized TSS flux was likely in the Onega and highly likely in the Kolyma. Declines in flow-normalized PHg flux were very or highly likely for the Yenisey (27.7 Mg⋅y−1), Ob’ (6.5), Northern Dvina (3.3), Pechora (1.9), and Mezen (0.3). For these five rivers, relative declines ranged from 68% (Ob’) to 92% (Yenisey), while flow-normalized PHg flux in the Onega likely did not change and flux could not be estimated for the Lena or Kolyma due to a more limited number of observations in these two rivers (Table 1 and ).

Table 1.

Absolute (total) and relative (percent) changes in mean annual flow-normalized concentrations and total annual flow-normalized fluxes of TSS and PHg, discharge, peak discharge magnitude, mean annual air temperature, and total annual precipitation

Parameter	River	n	Start (t1)	End (t2)	Value at t1	Value at t2	Change (total)	Change, %	Trend
TSS concentration, mg⋅L−1	Onega	149	1980	1992	12.8 (11.2–14.4)	16.1 (12.4–22.2)	3.3	25.8	L
	Northern Dvina	1,062	1979	2001	21.0 (15.5–23.8)	16.6 (13.3–19.4)	−4.5	−21.0	L
	Mezen	378	1980	2018	16.8 (13.1–22.4)	8.2 (7.0–10.5)	−8.7	−51.2	HL
	Pechora	415	1983	2018	22.7 (16.3–34.4)	6.8 (5.1–8.8)	−15.9	−70.0	HL
	Ob'	319	1974	2008	53.8 (48.4–58.7)	22.2 (17.1–25.6)	−31.8	−58.7	HL
	Yenisey	122	1980	2011	9.4 (7.2–10.5)	5.5 (2.4–7.9)	−3.9	−41.5	HL
	Lena	58	2004	2018	40.8 (18.0–99.9)	13.9 (0.4–18.6)	−26.4	−65.9	L
	Kolyma	55	1982	1987	19.8 (15.2–25.9)	39.7 (33.0–53.7)	19.8	100.5	HL
PHg concentration, µg⋅g−1	Onega	155	1980	1992	0.44 (0.41–0.48)	0.33 (0.28–0.36)	−0.11	−25.0	HL
	Northern Dvina	304	1980	2001	0.97 (0.89–1.11)	0.41 (0.37–0.49)	−0.57	−57.7	HL
	Mezen	223	1980	2002	0.67 (0.61–0.72)	0.30 (0.26–0.34)	−0.37	−55.2	HL
	Pechora	182	1983	2002	0.50 (0.47–0.55)	0.27 (0.23–0.31)	−0.23	−46.0	HL
	Ob'	222	1979	2008	0.59 (0.50–0.75)	0.21 (0.18–0.28)	−0.39	−64.4	HL
	Yenisey	146	1980	2011	2.0 (1.9–2.2)	0.57 (0.37–0.66)	−1.46	−71.8	HL
	Lena	40	2004	2011	—	—	—	—	—
	Kolyma	63	1982	1987	0.80 (0.68–0.90)	1.09 (0.97–1.27)	0.29	36.3	VL
TSS flux, Tg⋅y−1	Onega	149	1980	1992	0.26 (0.20–0.31)	0.48 (0.29–0.93)	0.21	84.6	L
	Northern Dvina	1,062	1979	2001	4.0 (2.5–4.8)	3.1 (2.1–3.8)	−0.85	−22.5	L
	Mezen	378	1980	2018	0.62 (0.44–0.99)	0.29 (0.20–0.40)	−0.34	−53.2	HL
	Pechora	415	1983	2018	5.3 (3.3–8.5)	1.0 (0.7–1.2)	−4.3	−81.1	HL
	Ob'	319	1974	2008	31.3 (28.3–34.7)	12.7 (10.2–15.1)	−19.2	−59.4	HL
	Yenisey	122	1980	2011	9.9 (7.0–11.8)	4.3 (1.7–5.1)	−5.6	−56.9	HL
	Lena	58	2004	2018	31.0 (15.6–37.6)	13.5 (2.3–17.9)	−17.4	−56.5	HL
	Kolyma	55	1982	1987	5.6 (3.9–7.0)	12.2 (9.4–16.8)	6.6	117.9	HL
PHg flux, Mg⋅y−1	Onega	126	1980	1992	0.10 (0.08–0.11)	0.12 (0.06–0.19)	0.02	20.0	N
	Northern Dvina	296	1980	2001	4.5 (2.7–5.7)	1.2 (0.9–2.0)	−3.3	−73.3	HL
	Mezen	183	1980	2002	0.37 (0.24–0.58)	0.05 (0.02–0.11)	−0.32	−86.5	HL
	Pechora	152	1983	2002	2.3 (1.4–4.0)	0.45 (0.19–0.73)	−1.9	−80.4	HL
	Ob'	222	1979	2008	9.7 (6.7–16.3)	3.1 (1.8–6.3)	−6.5	−68.0	HL
	Yenisey	122	1980	2011	30.0 (20.0–38.8)	2.3 (0.5–3.8)	−27.7	−92.3	HL
	Lena	37	—	—	—	—	—	—	—
	Kolyma	43	—	—	—	—	—	—	—
Discharge, km3⋅y−1	Onega	—	1980	1992	14.8	16.4	1.6	10.6	0.45
	Northern Dvina	—	1979	2001	99.4	107.2	7.8	7.9	0.65
	Mezen	—	1980	2018	19.1	20.5	1.4	7.2	0.92
	Pechora	—	1983	2018	132.4	102.4	−29.9	−22.6	0.02
	Ob'	—	1974	2008	390.2	408.2	18.1	4.6	0.55
	Yenisey	—	1980	2011	567.8	653.6	85.7	15.1	< 0.001
	Lena	—	2004	2018	639.3	629.7	−9.6	−1.5	0.44
	Kolyma	—	1982	1987	85.7	118.2	32.5	37.9	0.46
Peak discharge magnitude, km3	Onega	—	1980	1992	3.7	4.1	0.4	11.5	0.24
	Northern Dvina	—	1979	2001	28.1	32.8	4.8	17.0	0.18
	Mezen	—	1980	2018	7.0	6.9	−0.06	−0.9	0.74
	Pechora	—	1983	2018	38.5	30.7	−7.8	−20.2	0.32
	Ob'	—	1974	2008	57.4	62.0	4.6	7.9	0.27
	Yenisey	—	1980	2011	168.5	166.0	−2.5	−1.5	0.89
	Lena	—	2004	2018	145.0	157.6	12.6	8.7	0.66
	Kolyma	—	1982	1987	8.8	36.2	27.4	313.2	0.81
Mean annual air temperature, °C	Onega	12	1980	1992	1.0	2.1	1.1	112.1	0.84
	NorthDvina	21	1980	2001	0.5	1.1	0.6	133.3	0.88
	Mezen	22	1980	2002	−1.3	−1.4	−0.2	−13.4	0.34
	Pechora	19	1983	2002	−6.5	−7.7	−1.2	−18.3	0.94
	Ob	29	1979	2008	−0.7	1.6	2.3	333.3	0.05
	Yenisey	31	1980	2011	−5.7	−3.3	2.4	42.9	0.28
	Lena	7	2004	2011	−9.3	−7.1	2.2	23.2	1.00
	Kolyma	5	1982	1987	−12.3	−13.3	−1.0	−8.2	0.81
Total annual precipitation, mm	Onega	12	1980	1992	567	700	133	23.5	0.95
	NorthDvina	21	1980	2001	696	674	−22	−3.2	0.38
	Mezen	22	1980	2002	567	664	97	17.1	0.26
	Pechora	19	1983	2002	897	893	−4	−0.4	0.44
	Ob	29	1979	2008	680	572	−108	−15.9	0.81
	Yenisey	31	1980	2011	591	566	−25	−4.2	0.56
	Lena	7	2004	2011	503	471	−32	−6.4	0.23
	Kolyma	5	1982	1987	185	159	−26	−14.1	0.46

Trends in TSS and PHg were evaluated using WRTDS and bootstrapped 90% confidence intervals are shown in parentheses. WRTDS modeling requires n > 50 and was not done for the Lena (PHg concentration and flux) or Kolyma (PHg flux). Trend significance for WRTDS results are equivalent to two-sided P values and are reported as highly likely (HL ≥.95 and ≤.0), very likely (VL ≥.90 and < 0.95), likely (L ≥ 0.66 and < 0.90), about as likely as not (n ≥ 0.33 and < 0.66), or unlikely (≤0.33). Trends in discharge and watershed climate were evaluated using the Mann–Kendall trend test (significant P ≤ 0.05). Decreases are preceded by minus signs.

Underlying these trends in total annual fluxes, mean annual flow-normalized concentrations of TSS (milligrams per liter) and PHg (micrograms per gram) declined on average by 51% ([TSS]) and 53% ([PHg]) in all rivers except for the Onega (only PHg declined) and Kolyma (Fig. 2). Measurements of [PHg] from this study (mean = 0.7 µg⋅g−1, range = 0.1 to 5.3, n = 1,484) overlapped with observations for global rivers (0.4 [0.01 to 3.7], n = 168) (). All declines in flow-normalized [PHg] and most for flow-normalized [TSS] were highly likely (Table 1). Except for a 15% increase in annual discharge in the Yenisey and a 23% decrease in the Pechora, declines in flow-normalized [TSS] and [PHg] were not accompanied by significant changes in discharge or peak discharge magnitude in any other rivers, the latter representing a potential driver of fluvial erosion (33). Results derived for the Yenisey and Pechora using generalized flow normalization, which accounts for long-term change in discharge (see ), indicate that declines in flow-normalized TSS and PHg fluxes were mostly (91 to 93% and 91 to 96%, respectively) due to changes in concentration. Hence, declines in TSS and PHg flux in recent decades appear to be primarily linked to diminished inputs of sediment and Hg to these rivers and/or increased sedimentation within fluvial networks (see below).

Fig. 2.

(A–H) Mean annual measured concentrations (± SE) of PHg (primary y axes) and TSS (secondary y axes) by river, shown as points superimposed on flow-normalized concentrations (smoothed line) from WRTDS (see ). Shaded regions are the 90% confidence interval around flow-normalized concentrations, derived using a WRTDS bootstrap approach for estimating uncertainty. Uncertainty could not be estimated for years with insufficient records of daily discharge (Onega and Northern Dvina) or low sample size (PHg, Lena).

These trends extend previous assessments of mid- to late-20th-century sediment fluxes in major Russian rivers, which showed decreased TSS flux in the Yenisey and variability in the Lena and Ob’ (34–36). Our observation that flow-normalized TSS flux in the Kolyma increased within only 6 y from 1982 to 1987 may reflect interannual variability due to a shorter sampling duration but agrees with long-term observations (ca. 1940 to 1988) of increased TSS flux (34, 35). Although our estimates of PHg flux extend to ca. 2010, independent measurements of total Hg (PHg + dissolved Hg) and TSS flux made from 2012 to 2017 in the Ob’, Yenisey, and Lena (3, 4) align with trajectories of declining flow-normalized PHg flux that we document in preceding decades (). Our measurements bridge observations between the mid-20th century and 2010s, revealing large and previously undocumented declines in TSS and PHg concentration and flux that are indicative of striking changes in Arctic fluvial sediment regimes and Hg cycling during recent decades.

Spatial Variability in PHg and TSS Export and Drivers.

Normalizing constituent fluxes by watershed area (i.e., yield) facilitates assessment of controls on constituent export from underlying terrain differences (e.g., surficial geology), hydrology, and anthropogenic impacts across diverse landscapes (37). We investigate TSS and PHg yields as a first step in characterizing potential drivers of declining Hg export among watersheds (continued in the next section). Fluxes we used to calculate yields were estimated by WRTDS Kalman filtering, which accounts for serial autocorrelation and provides non-flow-normalized flux estimates that better capture temporal variation and discharge-driven patterns of constituent export (38). Mean annual TSS yields ranged from 0.7 to 17.2 Mg⋅km−2⋅y−1 and varied significantly among watersheds (ANOVA: F1,7 = 135, P < 0.001). Higher TSS yields in the Kolyma (17.2 ± 4.4 Mg⋅km−2⋅y−1, mean ± SE) and Lena (9.1 ± 0.8) than in the western watersheds (6.0 ± 1.4) have been attributed to relatively erodible Quaternary sediments bearing interstratified layers of ice in eastern Siberia (39). In the western Russian watersheds, sedimentation within expansive downstream lowlands is thought to sequester substrate derived from erosion in mountainous upper reaches, reducing TSS export to the coastal ocean (39). These trends highlight a broad longitudinal gradient in topography and surficial geology underlying regional variation in TSS yields across the Russian north. Superimposed on this gradient, large reservoirs constructed on the Ob’, Yenisey, and their tributaries (Fig. 1) in the mid-20th century reduced TSS export in the decades following by trapping large amounts of sediment (34), consistent with global reductions in riverine sediment export due to reservoir establishment (19). Our observations of large reductions in late-20th-century TSS flux (Table 1) extend these earlier trends and implicate reservoirs as a primary sediment sink in Russian fluvial networks during the mid to late 20th century.

Mean annual PHg yields varied significantly among watersheds (F1,5 = 56, P < 0.001) and overall were relatively low (0.4 to 6.4 g⋅km−2⋅y−1), which is typical of watersheds at northern high latitudes (2), where a cold climate tempers fluvial erosion for much of the year. In comparison to TSS, PHg yields were highest in the Yenisey and Northern Dvina (5.1 ± 0.7 to 6.4 ± 0.7), lowest in the Mezen (0.4 ± 0.1), and intermediate in the other rivers (2.1 ± 0.2 to 4.5 ± 0.7) (there are no data for Lena and Kolyma due to shorter periods of record). While regional variability in TSS yields appeared to be driven by longitudinal variability in terrain conditions and the distribution of reservoirs, relatively high PHg yields in the Yenisey and Northern Dvina correspond with records of historical Hg contamination from regional industry. For instance, major pulp and paper mills along the lower Northern Dvina accumulated large amounts of Hg waste adjacent to their facilities (20) from manufacturing of bleaching agents for paper production. Chlor-alkali production in headwaters of the Yenisey released more than 1,200 Mg of Hg into the environment during operation from 1973 to 1998 (40) and soil Hg concentrations in contaminated regions reach 125 µg⋅g−1 (41). In addition to likely proximate sources of Hg to the rivers in this study, rivers integrate biogeochemical effects of anthropogenic disturbance from across their watersheds. In 2010, an estimated 14 Mg of Hg was released to Russian freshwaters from anthropogenic sources (effluent released directly and indirectly via soil erosion and leaching from contaminated sites) and 7 Mg from background sources (soil erosion and runoff of natural and atmospherically deposited Hg) (17). From a mass budget perspective, these recent inventories together with our hydrochemical observations estimate comparable magnitudes of Hg input to aquatic environments from anthropogenic sources (14 Mg⋅y−1 in 2010) and transported within Russian rivers (∼19 Mg⋅y−1 ca. 2000s; Table 1). The spatial patterns and overall large declines in PHg export we document surround a large decrease in industrial activity from the late 1980s to 1990s as the Soviet Union collapsed (24), which may have partly reduced Hg loading to rivers via declines in industrial effluent and emissions to the atmosphere.

Atmospheric Hg Deposition to Russian Watersheds.

The onset of Hg emission controls during the mid-20th century and associated declines in deposition from primary anthropogenic emissions (42) are thought to have reduced Hg loading to Russian ecosystems and also riverine Hg export (16). To test this hypothesis, we modeled decadal trends in deposition of primary anthropogenic Hg emissions from regional (i.e., Europe and USSR) and global sources to the eight focal watersheds, using a global chemical transport model (GEOS-Chem v.12.7) (29) and published inventories of Hg emission magnitudes (25, 43, 44). Global primary anthropogenic Hg emissions declined from a peak of 3,020 Mg⋅y−1 in 1970 to a minimum of 1,960 Mg⋅y−1 in 1990, before rising to 2,390 Mg⋅y−1 in 2015. Within Europe and the former Soviet Union, primary anthropogenic emissions declined by 1,070 Mg⋅y−1 (from 1,230 to 160) between 1970 and 2015. During this time, mean annual Hg deposition from regional emissions to all focal watersheds declined by 85 Mg⋅y−1. The decline in Hg deposition from global emissions, which includes regional sources, was smaller (79 Mg⋅y−1), largely due to increased emissions from eastern Asia. Declining primary anthropogenic emissions in recent decades are consistent with the trends we observed in Hg deposition to the focal watersheds (Fig. 3 and ).

Fig. 3.

(A–H) Total deposition of Hg [as dry Hg(0) + wet and dry Hg(II)] to the eight focal watersheds from 1960 to 2015, from global sources (solid line) and regional sources (Europe [EU] and Union of Soviet Socialist Republics [USSR]; dashed line). Deposition was calculated using regional Hg emissions totals from the literature (43) as inputs to a global three-dimensional model of atmospheric chemistry (GEOS-Chem) (29). Gray shaded regions span the years for which flow-normalized PHg concentration was estimated using WRTDS.

Hg uptake within terrestrial ecosystems and eventual export of Hg from soils to freshwaters can occur over decades to centuries (27). This export lag may be accelerated in northern watersheds if winter and springtime Hg deposition is preserved on the surface and flushed into fluvial networks via the freshet snowmelt pulse. Thus, comparing deposition and riverine Hg mass fluxes among watersheds facilitates assessment of potential contributions from changing emissions to fluvial export. In the four smaller western watersheds and also in the Ob’, declines in atmospheric Hg deposition (–15 to –17 Mg⋅y−1) exceeded changes in riverine PHg flux (–12 Mg⋅y−1) (Table 1). Normalized by watershed area, average peak deposition in the eastern watersheds (22 g⋅km−2⋅y−1; Ob’, Yenisey, Lena, Kolyma) was smaller than in the four western watersheds (46 g⋅km−2⋅y−1), where a larger proportion of Hg deposition originated from proximate regional urban and industrial centers rather than from global sources (Fig. 3). Together, these trends in mass fluxes and areal deposition indicate relatively greater effects from declines in regional primary anthropogenic Hg emissions on riverine Hg export in the smaller western watersheds. In the Ob’ and Yenisey, large decreases in flow-normalized PHg flux relative to deposition implicate reservoirs and reductions in local industrial pollution as primary drivers of declines in riverine PHg export. Together, our findings indicate that multidecadal declines in riverine PHg export across Russia were driven by decreases in [PHg] associated with varied environmental and anthropogenic factors.

Potential Drivers of Trends in PHg Concentration.

Dynamic regression models (DRMs) we developed for each river (see for multivariate time series analysis) provide evidence for potential drivers underlying the near-universal declines in flow-normalized [PHg]. Covariate coefficients indicated that changes in atmospheric Hg deposition from regional sources were most strongly associated with [PHg] in the four smaller western watersheds (Onega, Northern Dvina, Mezen, and Pechora) (Table 2). For these watersheds, deposition lagged at 10 to 15 y may reflect effects on aquatic Hg loading from both rapid (subdecadal) responses in direct deposition to surface waters and also delayed responses via runoff (45) and also capture the period from ∼1970 to 1990 of relatively rapid declines in deposition from regional sources which contributed more to Hg deposition (Fig. 3 and “yields” results). Global Hg deposition lagged at 5 y was most strongly associated with [PHg] trends in the Ob’ watershed, where deposition to surface waters in vast wetlands (46) may accelerate the response of aquatic Hg cycling to changes in declining atmospheric deposition. Climate covariates were included in the most parsimonious models for most of the five western rivers, despite showing negligible change over the study period (Table 1) and generally appearing less-associated with flow-normalized [PHg] than atmospheric Hg deposition and flow-normalized [TSS]. The association between flow-normalized [TSS] and [PHg], an indicator of sediment–Hg coupling, appeared to vary among rivers and was strongest in the Yenisey, where the top two DRMs included terms for flow-normalized [TSS] and lagged cumulative reservoir capacity, respectively. In the most parsimonious DRM, omission of reservoir capacity suggests the covariate for [TSS] captured reservoirs’ effects on sediment concentrations across the Yenisey watershed (Fig. 1). Although the short periods of record precluded assessment of [PHg] trends and drivers in the Lena and Kolyma, accelerating permafrost thaw is likely to emerge as a key driver of fluvial PHg export in these watersheds (47).

Table 2.

Results from DRMs used to evaluate potential drivers underlying trends in flow-normalized PHg concentration

River	BIC	Covariate	Coefficient	SE
Onega	−122.8	Regional Hg dep (15-y lag)	0.009	0.002
		[TSS]	0.007	0.003
		TAP	−0.001	0.002
		Intercept	−0.011	0.001
		AR1	−0.393	0.304
		σ	0.001
Northern Dvina	−144.6	Regional Hg dep (15-y lag)	0.020	0.008
		[TSS]	−0.013	0.003
		MAAT	0.001	0.001
		AR1	0.985	0.020
		σ	0.006
Mezen	−197.4	Regional Hg dep (10-y lag)	−0.042	0.019
		TAP	0.009	0.002
		[TSS]	−0.006	0.003
		AR1	0.865	0.106
		σ	0.002
Pechora	−171.4	Regional Hg dep (10-y lag)	0.008	0.006
		[TSS]	0.008	0.002
		AR1	−0.008	0.245
		σ	0.002
Ob'	−234.1	Global Hg dep (5-y lag)	0.003	0.003
		TAP	−0.002	0.006
		[TSS]	0.002	0.002
		AR1	0.960	0.035
		σ	0.003
Yenisey	−157.4	[TSS]	0.047	0.016
		MAAT	0.003	0.002
		Intercept	−0.041	0.005
		AR1	0.466	0.197
		σ	0.015

For each river, a suite of preliminary models contained all possible combinations of covariates: atmospheric Hg deposition (dep) (global vs. local with 0-, 5-, 10-, and 15-y lags), the flow-normalized concentration of total suspended solids ([TSS]), watershed total annual precipitation (TAP), watershed mean annual air temperature (MAAT), and (for the Ob’ and Yenisey only) cumulative reservoir capacity (20- and 25-y lags). The most parsimonious model was selected as the model with the lowest Bayesian information criterion (BIC) score that did not violate model assumptions. AR1 (regression residuals) and σ (model residuals, calculated as √σ2) are components of the ARIMA error. The ARIMA structure was (1,2,0) for the Mezen and Pechora, and (1,1,0) for all other rivers (see ).

Other processes known to influence riverine PHg, such as local and regional industrial pollution, were not quantified well enough on a watershed basis to include as covariates in our models, but their effects are implicated by historical records and reflected by DRM residuals. Nevertheless, our modeling results suggest that reductions in atmospheric Hg deposition and likely also industrial polluting activities within watersheds drove declines in [PHg] which, together with sediment trapping by reservoirs (especially in the Yenisey), contributed to the large declines in Russian river Hg export in recent decades.

Implications for Northern High-Latitude Hg Cycling.

Trajectories of declining PHg export together with estimates of flux magnitude enable a first assessment of historical riverine PHg fluxes. Summing the initial timepoint values from WRTDS results (Table 1) shows that riverine flow-normalized PHg fluxes from the Russian rivers (excluding the Lena and Kolyma) were ∼47 Mg⋅y−1 around 1980. This value is 7 to 27% higher than estimates of 37 to 44 Mg⋅y−1 for contemporary (2012 to 2017) riverine total Hg flux upscaled to the entire pan-Arctic watershed (3, 4). Total Hg fluxes in the late 20th century were likely greater still, due to fluxes of dissolved Hg for which we are unable to account. Contemporary fluxes of dissolved Hg (<2.2 µm passing fraction) in these Russian rivers are considered significant relative to PHg (49% ± 2%) (3), suggesting that—under a potential scenario of limited historical changes in the proportion of the Hg pool in the dissolved versus particulate phase—total Hg flux by Russian rivers approached 100 Mg⋅y−1 some four decades ago. Major sources of Hg to the Arctic Ocean and contemporary estimates of their inputs include ocean currents (21 to 48 Mg⋅y−1), coastal erosion (15 to 47 Mg⋅y−1), and atmospheric deposition and snowmelt (45 to 98 Mg⋅y−1) (5, 6, 48). The latter may have declined since the onset of Hg controls in the 1970s (25). Therefore, we estimate that Russian rivers were potentially the single largest source of Hg to the Arctic Ocean in the mid-20th century. Further, contemporary total annual Hg flux from the Yukon and Mackenzie Rivers is ∼10 Mg⋅y−1 (4) and has potentially increased in recent decades (15). Thus, changes to Arctic fluvial Hg cycling during the last half century were driven by diminishing terrestrial Hg inputs within Russian rivers in the pan-Arctic basin.

Future permafrost thaw is expected to intensify Hg cycling (13, 30). Across the circumpolar north, diverse terrain conditions (e.g., relief, ground ice, and Hg content) are likely to drive regional variability in the effects of permafrost thaw on fluvial export and biogeochemical transformation of Hg (49). Particularly in watersheds with limited flow regulation and large stores of soil Hg, permafrost degradation and increased sediment transport associated with melting of regional glaciers (7, 50) are likely to emerge as prominent sources of fluvial Hg as the Arctic warms (30). In North America, the most severe trajectories of climate warming and associated permafrost thaw project a doubling of fluvial Hg flux to nearly 10 Mg⋅y−1 in the Yukon River Basin by 2100 (30). In Canada, where large areas are susceptible to hillslope thermokarst activity (51), fluvial propagation of thawed substrate may redistribute Hg across watershed scales for centuries to millennia (14, 52). Permafrost degradation in western Siberia is predicted to double fluvial PHg flux in the coming decades (13). The degree to which these effects may counteract our observations of decreased riverine TSS and Hg export will likely be determined by anthropogenic activities (e.g., industrial activity and dam construction) that modify sediment transport and Hg loading to watersheds across the north.

In this study, we found significant declines in sediment and Hg concentration and export in major Russian rivers that were driven by pollution reductions and reservoir establishment, resolving a decades-long period of uncertainty in the magnitude and drivers of Hg cycling across a vast portion of the Arctic. As anthropogenic polluting activities continue to change and permafrost thaw intensifies land–freshwater linkages (13, 14, 30), hydrochemical and biological monitoring programs building on continuous, long-term data collection will be imperative to capture changes in Hg cycling at northern high latitudes (53). Future research should also prioritize 1) investigating the fate of PHg within northern reservoirs, lakes, and coastal regions, where large accumulations of Hg in sediments (54) may influence methylmercury dynamics and accumulation within food webs; 2) modeling of permafrost region Hg cycling (30) which incorporates potential effects from regional human activities on Hg mass budgets; and 3) tracing sources of fluvial Hg, to discriminate between permafrost and nonpermafrost Hg sources and to elucidate Hg cycling (55).

Materials and Methods

Sampling and Analysis.

River water samples for TSS and PHg concentration were collected approximately monthly and sampling duration varied across all rivers from 5 to 36 y ( and Table S2). River water was collected at midchannel and at the left and right banks from near the surface (0- to 0.5-m depth) and from midchannel only in November, December, February, March, and April during frozen conditions. Samples were collected using a metal-free sampler, stored in precleaned 2- to 5-L borosilicate glass bottles (Teflon-coated bottles in the year 1980) following clean sampling techniques, sealed, and kept chilled until analysis. Laboratory methods for processing of TSS and PHg samples followed clean techniques from previous studies (15, 16, 56–59) and was consistent throughout the duration of this study. Within 2 h of collection, up to 2 L of water was filtered using acid-cleaned polycarbonate Nuclepore filters (0.45 μm). Filters were dried at 105 °C and weighed to obtain TSS (detection limit = 3 mg⋅L−1) (60). For PHg, additional sediment was filtered, digested, reduced, and quantified using cold vapor atomic absorption spectrometry (PerkinElmer Coleman Mas-50 Mercury Analyzer or Spectrophotometer VARIAN AA-475) (61). Certified reference materials for total Hg were obtained from the National Research Council of Canada, the National Institute of Standards and Technology, and the International Atomic Energy Agency, Monaco (certified reference material range: 0.061 to 4.64 μg⋅g−1). PHg (nanograms per liter) used in flux estimates was calculated as the product of PHg (micrograms per gram) and TSS (milligrams per liter). Additional details on quality assurance and control, along with watershed characteristics related to reservoirs and industrial pollution, are available in .

River Discharge and Watershed Climate.

Observations of daily discharge were obtained from ref. 62. We calculated annual peak discharge magnitude (33) as the sum of daily discharge in the 95th percentile for each year, to assess annual trends in fluvial energy, implications for erosion and sediment transport, and potential effects from anthropogenic activity (e.g., river flow regulation) in the context of changing northern hydrology (9). We tested temporal trends in annual discharge, peak discharge magnitude, watershed mean annual air temperature, and watershed total annual precipitation using the nonparametric Mann–Kendall and trend prewhitening to account for serial correlation (63) using the zyp package (64) in R software v.3.6.2 (65). (See for information on climate data sources.) Additional details on discharge measurements are available in .

Hg and Sediment Flux Modeling.

Estimates of annual TSS and PHg concentrations and fluxes were made using the United States Geological Survey weighted regressions on time, discharge, and season (WRTDS) modeling approach (31) in the R package Exploration and Graphics for RivEr Trends (EGRET) v.3.0 (66). A major benefit of WRTDS for modeling long-term trends in hydrochemical constituent concentrations and fluxes is that it allows for flexibility in the relationship between discharge and constituent concentration, resulting in reduced bias and improved accuracy of flux estimates (31). WRTDS outputs include flow-normalized concentration and fluxes, while an extension (WRTDS_K) provides Kalman filter (i.e., non-flow-normalized) flux estimates (38). We show mean annual measured concentrations (Fig. 2) and total annual Kalman fluxes () to illustrate interannual variation in TSS and PHg, and we use Kalman fluxes to estimate constituent yields (flux/watershed area), which facilitates assessment of controls on constituent export from underlying terrain differences (e.g., surficial geology), hydrology, and anthropogenic impacts across diverse landscapes (37). Watershed areas upstream of each discharge monitoring station were obtained from (67) and the ArcticRIMS Project (https://rims.unh.edu). A bootstrap method was used to estimate uncertainty in flow-normalized concentrations and fluxes and to test for statistical significance and strength of trends (32). Trends in river discharge were also tested over the PHg concentration interval and, for rivers exhibiting a significant change in discharge (Pechora and Yenisey only), flow-normalized concentrations, fluxes, and their bootstrap uncertainties were derived using a generalized flow normalization procedure detailed for v.3.0 of EGRET (66) (http://usgs-r.github.io/EGRET/index.html). This procedure also determines the change in flow-normalized concentration and flux that is attributable to discharge versus concentration. Total absolute and relative (percent) changes in TSS and PHg concentrations and fluxes were given as outputs from WRTDS. No issues in model performance were detected during visual inspection of model fit and bias in residuals. Prior to TSS and PHg modeling, corrections for spatiotemporal offsets between chemistry sampling and discharge monitoring stations for the Northern Dvina, Mezen, Pechora, Yenisey, Lena, and Kolyma Rivers were applied following ref. 67.

Atmospheric Hg Deposition 1960 to 2015.

Simulations of atmospheric Hg deposition trends from local (Europe and USSR) and global sources to the eight focal watersheds were conducted using the GEOS-Chem global chemical transport model for Hg (version 12.7; https://geos-chem.seas.harvard.edu/) (29). The model calculates the transport, chemistry, and deposition of atmospheric Hg species and has been extensively evaluated against all available observations (29). Atmospheric transport is driven by assimilated meteorological data from the NASA Global Modeling and Assimilation Office’s MERRA-2 reanalysis product (68). All simulations were conducted at the 2° × 2.5° horizontal resolution, with 47 vertical levels. The model is forced by primary anthropogenic emissions for the period 1960 to 2015 based on published inventories (25, 43, 44). These inventories account for the role of emission control technologies in changing the speciation of power plant emissions, allowing us to better estimate effects on regional deposition which is, to a large degree, determined by Hg(II) emissions. Emissions from terrestrial ecosystems, ocean evasion, and geogenic emissions were held constant at the modern value of 5,260 Mg⋅y−1 in the GEOS-Chem benchmark budget. We used magnitudes of atmospheric Hg deposition together with inventories on industrial Hg release to assess potential contributions from pollution sources to riverine Hg export. Additional details on reconstructing historical deposition trends are available in .

Multivariate Time Series Analysis.

To assess potential drivers underlying trends in PHg concentration in each river, we developed DRMs using the fable package (69) in R. We accounted for time series autocorrelation in each model by including an autoregressive integrated moving average (ARIMA) error structure, which contains the components (p,d,q). Components were determined using a Durbin–Watson test (for p, the autoregressive order) and a Kwiatkowski–Phillips–Schmidt–Shin test (for d, the degree of first-order differencing to achieve stationarity); q, the moving average component, is often used for modeling data with intra-annual (e.g., seasonal) trends (70) and was thus set to 0 for our time series of annual data. ARIMA error structures were either (1,1,0) or (1,2,0). We included predictor covariates for changes in atmospheric Hg deposition, suspended sediments, reservoirs, and climate. Potential delayed effects on PHg concentration from changes in Hg deposition (27) and reservoir establishment (19) were modeled by including lagged values of atmospheric Hg deposition (lags of 5, 10, and 15 y) and cumulative reservoir capacity (20 and 25 y). Lag spans were limited by data availability and reservoir effects were assessed only for the Ob’ and Yenisey, because datasets were not available to evaluate long-term (>10 y) trends in PHg for the Lena and Kolyma. Mean annual air temperature and total annual precipitation were derived from ERA5 monthly reanalysis data (0.25° × 0.25°) (28) available from 1979 to 2015 using Google Earth Engine (71). Reservoir locations, volumes, and dates of establishment were obtained from the Global Reservoir and Dam Database (18) and the literature (). The preliminary suite of DRMs included all possible combinations of predictors, with one term each for deposition and reservoirs to differentiate between contemporary and lagged effects. The model with the lowest Bayesian information criterion was selected as the most parsimonious. A Ljung–Box test and plots of model residual autocorrelation and distribution were used to confirm that residuals were indistinguishable from white noise (i.e., uncorrelated, homoscedastic, and mean of zero), and models violating these conditions were not considered.

Additional Statistics.

Mean annual yields of PHg and TSS among rivers were compared using a Welch’s heteroscedastic F test from the R package onewaytests (72). Bartlett’s test was used to first confirm that variances were not homogenous.