Isotopic evidence for acidity-driven enhancement of sulfate formation after SO2 emission control.

After the 1980s, atmospheric sulfate reduction is slower than the dramatic reductions in sulfur dioxide (SO2) emissions. However, a lack of observational evidence has hindered the identification of causal feedback mechanisms. Here, we report an increase in the oxygen isotopic composition of sulfate ([Formula: see text]) in a Greenland ice core, implying an enhanced role of acidity-dependent in-cloud oxidation by ozone (up to 17 to 27%) in sulfate production since the 1960s. A global chemical transport model reproduces the magnitude of the increase in observed [Formula: see text] with a 10 to 15% enhancement in the conversion efficiency from SO2 to sulfate in Eastern North America and Western Europe. With an expected continued decrease in atmospheric acidity, this feedback will continue in the future and partially hinder air quality improvements.

INTRODUCTION

Atmospheric sulfate () has a notable but uncertain impact on the global radiation budget and cloud lifetimes (). Sulfate also accounts for a major component of fine particulate matter mass in urban regions, affecting visibility () and public health (). Since the Industrial Revolution, increased emissions of sulfur dioxide (SO2) have resulted in an increase in sulfate load. The period from the 1950s to the 1970s had increased SO2 emissions leading to high-pollution; however, switching to cleaner technology and fuels decreased these emissions after the 1980s across North America (NA) and Western Europe (WE) (). These decreases successfully lowered sulfate concentrations, which avoided hundreds of thousands of deaths and illnesses from exposure to particulate matter smaller than 2.5 μm (PM2.5) in the United States alone (). However, atmospheric sulfate declined less rapidly than SO2 emissions, especially in wintertime (–). This unexpected phenomenon suggests the existence of feedback processes, which render SO2 emission reductions less efficient than expected for mitigation of sulfate aerosol.

In the atmosphere, the majority (60 to 80%) of sulfate formation occurs through oxidation of SO2 in the aqueous phase (i.e., in clouds), with gas-phase SO2 oxidation by hydroxyl radicals (•OH) accounting for most of the remainder (20 to 40%) (). Upon dissolution in the aqueous phase, SO2 dissociates into S(IV) species (mainly HSO3− and ) that react with hydrogen peroxide (H2O2), ozone (O3), molecular oxygen (O2), and hypohalous acid (e.g., HOBr) to form sulfate (). S(IV) oxidation is thought to be dominated by H2O2 oxidation, which is largely insensitive to cloud acidity, whereas the more minor S(IV) + O3 pathway exhibits strong pH dependence (). Among these chemical processes, several studies have proposed enhancement of S(IV) + H2O2 (, ), S(IV) + O3 (), or SO2 + •OH (, ) oxidation pathways as mechanisms causing the weakened response of sulfate abundances to decreases in SO2 emissions. These conclusions are based on a comparison between observed and modeled sulfate concentrations, and model estimates of atmospheric sulfate production mechanisms. However, observational evidence so far has not pointed to a specific mechanism. These uncertainties limit confidence in the model forecast and hindcasts of management of current and future tropospheric sulfate aerosol and its environmental impacts.

One tool for providing insight into sulfate formation mechanisms is offered by the mass-independent oxygen isotopic composition (Δ17O) () of sulfate (see Materials and Methods). Δ17OSO2 equals 0‰ because of rapid oxygen exchange between SO2 and H2O in the atmosphere (), and thus the value reflects the oxidation pathway of SO2 to . Gas-phase SO2 oxidation by •OH, where Δ17O•OH = 0‰, produces with = 0‰ (). Aqueous-phase S(IV) oxidation by H2O2 and O3 leads to nonzero values through transfer of the nonzero Δ17O from the oxidants. Δ17OH2O2 has been reported to be 1.6 ± 0.3‰ (), and the transfer factor from oxidant to sulfate is 0.5 based on a laboratory experiment (), which yields values for S(IV) oxidation by H2O2 [] of 0.8 ± 0.2‰. For Δ17OO3, the two earliest studies using a cryogenic technique showed a large range over ±10‰ (, ). This range is much greater than expected from the pressure and temperature dependency of Δ17OO3 (–). On the basis of these experimental data, the large variability found in the two early studies using cryogenic techniques would be caused by random errors associated with sampling artifacts (). We therefore exclude these two studies from consideration and instead use the average value of the tropospheric Δ17O with 25.6 ± 1.3‰ originating from the observations using nitrite-coated method among different locations and seasons (–). As a consequence, the values for S(IV) + O3 [] is assumed to be 6.4 ± 0.3‰, using a transferring factor of 0.25 based on a laboratory study (). The S(IV) + O2 reaction catalyzed by trace metal ions (TMI) produces of −0.1‰ by transferring one oxygen atom transferred from atmospheric O2 [Δ17O ≈ −0.3‰ ()]. In-cloud HOBr leads to of 0‰, and primary from natural and anthropogenic sources also has = 0‰ ().

Consequently, is solely determined by the proportions of different sulfate formation pathways that yield nonzero values of ()where the f (x = O3, H2O2, TMI, and zero) terms indicate the respective fractions of SO2 oxidized by O3, H2O2, TMI-catalyzed O2, and oxidants which have = 0‰, and fO3 + fH2O2 + fTMI + fzero = 1 (see Supplementary Text). To date, this approach using has enabled observation-based quantification of atmospheric sulfate formation in different regions and time periods, including glacial-interglacial cycles (), stratospheric volcanic eruptions (), preindustrial biomass burning (), and transition after the Industrial Revolution (). In particular, comparison of observed and modeled has enabled the recognition and quantification of sulfate formation mechanisms often ignored in models, including S(IV) + O2 oxidation catalyzed by TMI (, ), sulfate formation by ozone oxidation on sea salt aerosol (), S(IV) oxidation by hypohalous acids (, ), and heterogeneous reactions in extreme haze events (, ). However, there is no record of that can provide information on changes in sulfate formation pathways in response to the reduction in air pollution following the implementation of governmental reduction policy such as the U.S. Clean Air Act of 1970.

Here, we present the first observations of changes in Northern Hemisphere between 1959 and 2015, based on a continuous and high-resolution ice core record from a high-elevation dome site in southeast Greenland called SE-Dome. The record from the past 60 years covers the high-pollution decades of the 1950s to 1970s as well as the substantial SO2 emissions-reduction period from the 1980s to the present (). The SE-Dome ice core preserves atmospheric aerosols that originate mainly from NA and WE, with no notable change in the air mass origin over the period (). We reconstructed in 3 to 6 years resolution with an accuracy of dating better than 2 months, which is given by precise age-depth scaling with the oxygen-isotope matching method ().

RESULTS

Increase of ice core over the past 60 years

The ice core ranges from 1.0‰ to 1.7‰ and shows a substantial increase throughout the record (Fig. 1A), with a ~0.4‰ difference (P < 0.05) between the average of 1960 to 1970 (1.14 ± 0.05‰, n = 4) and that of 2005 to 2015 (1.51 ± 0.19‰, n = 3). This increase in clearly indicates that the sulfate formation pathways responsible for sulfate preserved in the SE-Dome have changed from the 1960s to the present. Given that the signatures for the sulfate formation pathways are all lower than 0.8‰ except for the S(IV) + O3 pathway [ = 6.4 ± 0.3‰], the increase can reasonably be interpreted as the reflection of an increase in the relative importance of the S(IV) + O3 pathway.

Fig. 1

and chemical fluxes at the SE-Dome and GEOS-Chem model results during the last 60 years.

Open black circles represent ice core record, the closed black circle represents data from snow pit, and colored symbols represent model results for given years for TRAJ, ENA, and WE regions. The observed chemical fluxes and neutralization ratio were obtained from Iizuka et al. (). The thin lines represent observed data for each year, and the open circles with thick lines represent the weighted average flux data corresponding to sample resolution (see Supplementary Text). (A) record. Open black circles represent ice core record and the closed black circle represents data for shallow snow with 1σ uncertainty shown as error bar. Colored symbols represent annual-mean, mass-weighted average of tropospheric for given years. The shaded area for the modeled indicates the 1σ uncertainty. (B) flux and modeled annual-mean concentration normalized to 1973 (C) NH4+ flux and modeled annual-mean concentrations of NH3 + NH4+ normalized to1973, (D) neutralization ratio: NH4+/(2 + NO3−) for observation, and (NH4+ + NH3)/[2 + (NO3− + HNO3)] calculated from modeled, annual-mean tropospheric concentrations, and (E) H+ flux and modeled tropospheric annual-mean, cloud liquid water weighted, bulk cloud pH.

By applying a simple isotope mass balance method (Eqs. 2 and 3), we calculated the maximum and minimum contribution of oxidation by the S(IV) + O3 pathway (fO3, max and fO3, min)

The fO3, max is estimated from the two end-members mixing between the S(IV) + O3 reaction and the S(IV) + O2 reaction catalyzed by TMI, which has the lowest end-member of (i.e., −0.1‰). The fO3, min, on the other hand, is calculated based on the mixing between the S(IV) + O3 pathway and the S(IV) + H2O2 pathway. Using these equations, the lowest fO3, min and fO3, max values were, respectively, calculated to be 3.2 ± 0.9% and 16.6 ± 1.9% for years 1975 to 1977 ( = 0.98 ± 0.10‰; sample 6 in table S1). These fO3, min and fO3, max values, respectively, increase to 15.5 ± 4.2% and 27.2 ± 1.9% at maximum of 1.67 ± 0.09‰ (years 2011 to 2014, sample 15 in table S1). The increase of the relative importance of the S(IV) + O3 pathway can be explained by (i) enhancement of the S(IV) + O3 pathway itself and (ii) decrease of other oxidation pathways that have lower values.

First, we consider the possible reasons that could cause an enhancement of the S(IV) + O3 pathway. Given that high pH conditions promote the S(IV) + O3 pathway (), acidity changes may be an important factor accounting for the increase in . The best correlation between and other ice core measurements is found for the neutralization ratio, NH4+/(2 + NO3−) (Fig. 1D), where r = 0.80 (P < 0.01; table S1). In addition, the H+ flux (Fig. 1E), the ice acidity indicator, also shows a strong correlation with (r = −0.71, P < 0.01; table S1). The increase in the neutralization ratio and the decrease in the acidity result from the simultaneous decrease in flux (Fig. 1B) and increase in NH4+ flux (Fig. 1C). This is consistent with previous studies that observed an increase in the pH of precipitation in NA and WE after the 1970s () and in a Greenland ice core (), mainly because of the mitigation of SO2 emission and the simultaneous increase in NH3 emission from agricultural and industrial sectors (). Changes in O3 concentrations over the period may also contribute to changes in sulfate formation pathways. Tropospheric O3 concentrations in the free troposphere [~3 km above sea level (a.s.l.)] increased by 1 to 3 ppbv (parts per billion by volume) decade−1 from the 1970s to 2000s, but there is no significant subsequent increase (). Thus, an increase in tropospheric O3 might partially contribute to the increase of the S(IV) + O3 pathway before 2000 but is not consistent with the substantial increase of for the post–year 2000 period (Fig. 1A).

Second, we discuss the possibility of inhibition of other oxidation pathways that have low , because the decrease of other sulfate formation pathways could also increase the relative importance of the S(IV) + O3 pathway. Influence from changes in the contribution from SO2 + OH due to changes in tropospheric •OH concentrations is thought to be minimal, given that no significant decrease in •OH is observed between 1980 and 2010 (). Greenland ice core shows a post-1950 increase of H2O2 (), but this H2O2 increase is not consistent with the observed increase, given the values of only 0.8‰ from S(IV) + H2O2 reaction. An observed 0.6‰ decrease in in a Greenland ice core from 1900 to 1980 was partially attributed to an increase in the TMI-catalyzed S(IV) + O2 oxidation pathway resulting from increases in anthropogenic metal emission after the Industrial Revolution (). A decrease in anthropogenic metal emissions since 1980 would tend to decrease the importance of the TMI-catalyzed S(IV) + O2 oxidation pathway, resulting in an increase in as observed.

The GEOS-Chem model reproduces the observed increase of

To examine the relative importance changes in sulfate formation mechanisms on observed ice core , we simulate global tropospheric chemistry using the GEOS-Chem chemical transport model (version 12.5.0, DOI: 10.5281/zenodo.3403111) with an updated “offline” tagged-sulfate aerosol simulation (, , ). We performed model simulations using meteorology and anthropogenic emissions for the years 1960, 1973, 1986, 1999, and 2013 (see Materials and Methods). This model calculates bulk cloud pH based on the local concentrations of inorganic species and CO2 (g) and considers the impact of heterogeneity in cloud water pH within a cloud (). Because the main source region of sulfate aerosol to the Arctic region is transported from Eastern North America (ENA) and WE, we compare the ice core observations with modeled tropospheric in ENA (−90° to −60°E, 30° to 60°N, n = 56) and WE (15° to 40°E, 40° to 70°N, n = 72) (fig. S2). We also compute modeled in the region −90° to 30°E and 35° to 90°N based on computed 10-day back trajectories [n = 300, Trajectory (TRAJ) region hereafter] (fig. S2) ().

In Fig. 1A, the modeled, annual-mean, mass-weighted average of tropospheric for five periods is shown. The model calculates a decrease of from 1960 to 1973 and an increase of from 1973 to 2013 for all three regions for TRAJ, ENA, and WE. Note that the SO2 emission is still increasing from 1960 to 1975 (Fig. 1B), and the in the model shows a decrease (Fig. 1A). This decrease in is not clearly observed in the ice core record, but one data point covering the year 1975 (sample 6, table S1) shows the lowest with 0.98 ± 0.10‰ within the study period (Fig. 1A). Overall, the observed falls within modeled values in all years, consistent with a mixture of sulfate originating in ENA and WE. The increase in modeled is due to an increase in fO3 between 1973 and 2013 from 8.0 to 18.0% for the ENA region and 16.3 to 23.7% for the WE region, respectively (fig. S3), consistent with the observation-based estimate of an increase in fO3 from 3.2 to 16.6% (1960 to 1970s) to 15.5 to 27.2% (present).

GEOS-Chem reproduces not only the observed trend, but also other trends for modeled sulfate (Fig. 1B), NH4+ (Fig. 1C), neutralization ratio (Fig. 1D), and bulk cloud pH (Fig. 1E) between 1973 and 2013. Modeled bulk cloud pH increased by 0.5 pH units between 1973 and 2013 (Fig. 1E) and neutralization ratio also shows substantial increases between 1973 and 2013 (Fig. 1D), with good agreement with modeled changes in for ENA and WE (r > 0.75, P < 0.01; table S2). In addition, no increase in the modeled is found when cloud pH is set constant to pH 4.5 in the model (fig. S6). The TMI-catalyzed S(IV) + O2 pathway, on the other hand, decreases by ~10% in the model between 1973 and 2013, yielding a 0.1‰ increase in (fig. S3). Although this goes in the right direction, it is not large enough to explain the observed 0.7‰ increase in . Modeled O3 concentrations also increase over this same time period, but the increase from 1960 to 1986 is higher than from 1986 to 2013 (fig. S4) and correlation between O3 concentrations to the modeled is not significant (table S2). Overall, we conclude that the modeled increase in is mainly driven by an increase in the cloud water pH over the same time period.

DISCUSSION

Enhancement of sulfate formation efficiency

We use the model to further investigate how changes in sulfate formation pathways alter the conversion efficiency (η) from SO2 to using the following metricwhere P() indicates the tropospheric sulfate production rate from the oxidation of SO2 by oxidant x [i.e., OH, H2O2, O3, O2 (TMI), and hypohalous acids (H)], and S(SO2) represents total tropospheric source SO2 calculated from SO2 emission, photochemical SO2 production, and net transport (import − export). Although there is regional variability (fig. S7), the common feature in all regions is that η decreases until 1973 and increases from 1973 to 2013. The decreases in total η from 1960 to 1973 are caused mainly by decreases in ηO3 and ηH2O2 for ENA, and caused by decreases in ηO3 but compensated by increases in ηOH for WE, respectively (Fig. 2, A and B). By contrast, the 10 to 15% increases in total η from 1973 to 2013 found in ENA and WE are caused mainly by increases in ηO3 and ηH2O2 particularly in ENA, but are partially compensated by decreases in ηTMI and ηOH (Fig. 2, A and B). The modeled increase of η in ENA after 1999 wintertime (fig. S8B) is consistent with a weakened response of reduction of U.S. SO2 emission due to combination of a weakening H2O2 limitation on the S(IV) + H2O2 pathway () and cloud pH–induced promotion of S(IV) + O3 pathway at low SO2 ().

Fig. 2

The conversion efficiency (η) from SO2 to sulfate () for each formation pathway (colored bars), annual source SO2 [S(SO2)] (black dots), and annual production rate of sulfate [P()] (red dots) calculated in GEOS-Chem.

(A) ENA; (B) WE.

It is worth noting that there is a regional difference of changes in η and corresponding processes between ENA and WE after SO2 emission control. In contrast to ENA, the WE region has relatively high coal combustion and thus higher metal emissions, yielding a higher contribution from the TMI-catalyzed S(IV) + O2 oxidation in wintertime (fig. S8C). Increases of ηO3 and ηH2O2 for WE between 1973 and 2013 are 10 and 8%, respectively, but these increases were compensated by a decrease of ηTMI by 5% (Fig. 2B). Notably, the main contribution for total η for WE was switched from the TMI-catalyzed S(IV) + O2 oxidation in 1973 to the S(IV) + O3 oxidation in 2013 (Fig. 2B). An increase in ηO3 in WE thus plays a greater role for the increase of η compared to ENA, because of relatively high cloud pH conditions (Fig. 1E) and more limited oxidant concentrations at higher latitude. The relatively small decreases in ηOH and increase in ηH in both regions are likely driven by an increase in cloud water pH, which enhances SO2 solubility in clouds (). As a consequence, the increase in η over the ENA and WE regions kept sulfate reduction slower than the reduction of source SO2 during the past 60 years (Fig. 2, A and B) by increasing solubility of SO2 and promoting the acidity-dependent S(IV) + O3 pathway.

The increases in sulfate production efficiency (η) have so far been partially compensated by the reduction of other pathways, particularly in WE (Fig. 2B), but this offset may not be significant in the future. The importance of the S(IV) + O3 pathway, on the other hand, will continue to increase in the future because of the increase of cloud water pH by expected future growth of anthropogenic NH3 emission () and possible future usage of NH3 as hydrogen storage/transport medium in hydrogen energy technologies (). SO2 emissions have increased in other parts of the world (e.g., China and India) over the past decades, but they then have been decreasing more recently and are expected to continue to decrease in the coming decades (, ). Thus, the acidity-driven enhancement of atmospheric sulfate formation causing nonlinear sulfate responses to SO2 is expected to occur for the wide area in the world as long as NH3 emissions are also not simultaneously controlled.

In addition to the SO42− burden, changes in sulfate formation pathways have implication for aerosol climatic effects [i.e., size distributions, cloud condensation nuclei (CCN) concentrations and aerosol radiative effect]. In-cloud sulfate production is a potentially important source of accumulation mode–sized CCN due to chemical growth of activated Aitken particles and the enhanced coalescence of processed particles (). Thus, all-sky aerosol radiative forcing (difference between the 1970–1974 and 2005–2009 periods) over the North Atlantic is modeled to be +1.2 W m−2 with constant cloud water pH condition at 5.0, but the forcing increases to +5.2 W m−2 if pH is assumed to increase by 1.0 unit over this period (). A recent climate model shows that reducing SO2 emissions at high CO2 concentrations will significantly enhance atmospheric warming, which is important to consider within the context of the 1.5°C/2°C target in the Paris Agreement (). Given that there is regional difference for the chemical feedback process for sulfate formation over the time period in this study (Fig. 2), the future feedback driven by atmospheric acidity for sulfate formation should be predicted locally and globally to design effective mitigation policies for air quality and climate change.

Importance of atmospheric acidity on sulfate formation

In this study, we present ice core record providing the first observational evidence that atmospheric acidity has driven changes in sulfate formation pathways as well as enhanced sulfate production rates after SO2 emissions control in NA and WE. Such identification of key processes for the sulfate formation pathways caused by changes in atmospheric acidity provides confidence in model’s forecast for the global sulfur cycle and its relation to climate change, which has not been provided based on only the comparison between observed and modeled sulfate concentrations.

Because atmospheric acidity has a central role for aqueous chemistry in general, cloud/aerosol pH have implications for atmospheric chemistry including not only sulfate but also nitrate, ammonium, halogen, metals, and organic compound reactions that affect air quality, human health, ecosystem, and climate change [() and references therein]. However, even for sulfate formation, until the 1990s, cloud pH was considered too low to allow the S(IV) + O3 reaction to occur, leading modeling studies to often neglect this S(IV) + O3 oxidation mechanism for sulfate formation (, ). Although recent studies have proposed the importance of cloud water pH promoting the S(IV) + O3 pathway as a possible factor explaining the sublinear response of concentrations to the decrease in SO2 emission (, ), even the latest study () prescribes a constant cloud pH in the estimation for radiative forcing effect via aerosols. The relative contribution of the S(IV) + O3 pathway to total production remains uncertain, ranging from 1 to 50% among model results (), and suffers from a lack of observational evidence. In this context, the values obtained from both atmospheric observation and ice core records provide the validation of changes in sulfate formation pathways and its climatic effects in response to the reduction in air pollution such as the U.S. Clean Air Act of 1970.

Although the increasing trends in were found in both observations and model after SO2 emission control, several nonnegligible uncertainties for the modeled remain (see Supplementary Text). In addition, there are several factors controlling in the atmosphere, indicating that it is not possible to use as a simple proxy for reconstructing cloud water pH. Similarly, as for the pH assumption in chemical transport models, spatial and temporal variability in clouds, which are challenging to predict and represented differently across models, could contribute to some of this model variability (). Long-term records of are rare and so far limited, but it is possible to reconstruct temporal variations from various ice cores from Arctic and Antarctic ice sheets and from mountain glaciers in various regions. Moreover, while we only considered the semi-volatile species for the interactive pH calculation in this study (see Materials and Methods), an improved pH calculation in GEOS-Chem including contributions from dust alkalinity, sea salt aerosol alkalinity, and carboxylic acids shows consistency with annual mean precipitation pH observed at wide regions in NA, WE, and East Asia (). Given the heterogeneous distribution of cloud pH (, ) and short lifetime (~4 to 5 days) of sulfate aerosols (), improvements of the model constrained by spatiotemporal variations of are desirable in the future for model validation of sulfate formation pathways, sulfate burden, and prediction of aerosols-influencing climate change.

MATERIALS AND METHODS

Samples

This study is based on the 90.45-m-deep ice core obtained at a southeastern Greenland dome site (67.18°N, 36.37°W, 3170 m a.s.l., SE-Dome hereafter) (). The age-depth scale was determined by the oxygen-isotope matching method with an estimated age error of 2 months (). The ice core samples were stored in the cold room (−50°C) of The Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan.

The ion fluxes of the species used in Fig. 1 are from a previous study (). To match the resolution of ion fluxes and , we reanalyzed the data with our resolutions as summarized in table S1. In addition to ion fluxes, we calculated the neutralization ratio (N) and H+ flux [F(H+)] determined by ion balancewhere F(X) (nmol m3 year−1) represents the flux of X as determined by Iizuka et al. (), and nss refers to non–sea salt component.

Since sea salt sulfate aerosols () are of little importance to atmospheric sulfur oxidation processes (i.e., = 0‰), Δ17O values were corrected for their sea salt sulfate () component to obtain their corresponding content, using Eq. 7 below.

The molar ratio of []seawater/[Na+]seawater in seawater is 0.06 ().

Ice core cut

The ice sample was cut by a band saw in a cold room (−20°C) and decontaminated by cutting the outside ice with a ceramic knife to ~70% of the sample’s original weight in a class 10,000 clean booth. The sample (approximately 1 to 2 kg) was then melted in a precleaned bottle and shipped frozen (~−20°C) to the Tokyo Institute of Technology, Yokohama, Japan. The sample was then stored in a freezer kept at −30°C until the time of the experiments described below. The samples and corresponding years covered are summarized in table S1.

Oxygen isotope analysis of sulfate

Δ17O is the deviation from the linear approximation of δ17O = 0.52 × δ18O and defined as Δ17O = δ17O − 0.52 × δ18O (), where δ17,18O = [(17,18O/16O)sample/(17,18O/16O)reference − 1], and reference represents the composition of Vienna Standard Mean Ocean Water.

The measurement system for follows Savarino et al. (), with modifications described in our previous study (). was separated from other ions using ion chromatography and converted to H2SO4. One micromole of H2SO4 was then chemically converted to Na2SO4, and 1 ml of 30% H2O2 solution was subsequently added and dried up. Next, the Na2SO4 was converted to silver sulfate (Ag2SO4) using an ion exchange resin (). This Ag2SO4 powder was transported in a custom-made quartz cup, which was dropped into the 1000°C furnace of a high-temperature conversion elemental analyzer (TC/EA; Thermo Fisher Scientific, Bremen, Germany) and thermally decomposed into O2 and SO2. Gas products from this sample pyrolysis were carried by ultrahigh-purity He (>99.99995% purity; Japan Air Gases Co., Tokyo, Japan), which was first purified using a molecular sieve (5 Å) held at −196°C (). The O2 and SO2 gas products were carried through a clean-up trap (trap 1) held at −196°C to trap SO2 and trace SO3, while O2 was transported to another tubing trap (1/16 inch outer diameter) with a molecular sieve (5 Å) (trap 2) held at −196°C to trap O2 separately from the other gas products. The O2 was purified using a gas chromatograph, with a CP-Molsieve (5 Å) column (0.32 mm inner diameter, 30 m length, and 10 μm film; Agilent Technologies Inc., Santa Clara, CA, USA) held at 40°C, before being introduced to the Isotope-ratio mass spectrometry (IRMS) system to measure m/z = 32, 33, and 34. As discussed by Schauer et al. (), this method results in oxygen isotope exchange between the O2 products and the quartz cups having Δ17O of 0‰ (), as well as the quartz reactor, which shifts δ17O and δ18O, and thus Δ17O measurements. The shift in values was corrected by estimating the magnitude of the oxygen isotope exchange with quartz materials, whose Δ17O value is assumed to be approximately 0‰ (). Note that the SO42− δ17O and δ18O values here are relative values to our O2 reference gas. The shift in due to exchange between quartz cup and samples was corrected by replicate analyses (n = 12) of the standard B ( = 2.4‰) with three independent experimental batches of this study. The four standard B were measured in the run number of 1, 5, 8, and 10 for each batch. In this correction for isotopic analysis, SD (1σ) for the corrected values for standard B was 0.07‰, and this 1σ uncertainty is considered for the error of the isotopic measurement for this study.

Equation 8 is the isotope mass balance equation between ss- and with = 0‰where “total” is the quantity measured by ion chromatography, corresponding to the sum of ss- and components. The data for and are summarized in table S1. The overall observational error for is calculated to be at or smaller than ±0.1‰ for ice core samples (table S1), based on propagation of the errors for ion concentrations of [Na+] and (±5%) and isotopic analysis (±0.07‰).

GEOS-Chem model description and simulations

GEOS-Chem is a global three-dimensional model of atmospheric composition (www.geos-chem.org) originally developed by Bey et al. (). In this study, we use GEOS-Chem (version 12.5.0, DOI: 10.5281/zenodo.3403111) driven by assimilated meteorological fields from MERRA-2 reanalysis data product from NASA Global Modeling and Assimilation Office’s GEOS-5 Data Assimilation System. We simulate aerosol-oxidant tropospheric chemistry containing detailed HOx-NOx-VOC-ozone-BrOx chemistry (–). The model was run at 4° ×5° horizontal resolution and 47 vertical levels up to 0.01 hPa. The model was spun up for 1 year before each of the 5 years simulated. In the model, sulfate is produced from gas-phase oxidation of SO2 (g) by OH, aqueous-phase oxidation of S(IV) by H2O2, O3, HOBr, and O2 catalyzed by TMI, and heterogeneous oxidation on sea salt aerosols by O3 (, ).

The parameterization of the metal-catalyzed S(IV) oxidation is described in Alexander et al. (). We consider Fe and Mn, which catalyze S(IV) oxidation, in the oxidation states of Fe(III) and Mn(II). Dust-derived Fe ([Fe]dust) is scaled to the modeled dust concentration as 3.5% of total dust mass, and dust-derived Mn is a factor of 50 times lower than [Fe]dust. Anthropogenic Fe ([Fe]anthro) is scaled as 1/30 of primary sulfate and anthropogenic Mn ([Mn]anthro) is 10 times lower than that of [Fe]anthro. In the model, 50% of Mn is dissolved in cloud water as Mn(II) oxidation state, and 1% of [Fe]dust and 10% of [Fe]anthro are dissolved in cloud water as Fe(III) oxidation states. This parameterization might underestimate the anthropogenic Fe and Mn, especially for the U.S. region, which is discussed in Supplementary Text.

For pH-dependent S(IV) partitioning, bulk cloud water pH is calculated iteratively using concentrations of sulfate, total nitrate (HNO3 + NO3−), total ammonia (NH3 + NH4+), SO2, and CO2 = 380 ppmv (parts per million by volume) based on their effective Henry’s law constants and the local cloud liquid water content as described in Alexander et al. (). We use the Yuen et al. parameterization to account for the effect of heterogeneity of cloud water pH on S(IV) partitioning and subsequent aqueous phase sulfate formation (). Sulfate formed from each oxidation pathway was treated as a different “tracer” in the model to calculate Δ17O as described elsewhere (, ).

Given that meteorological fields from MERRA-2 are not available before 1979, we use meteorological fields and nonanthropogenic emissions such as biogenic VOCs (), soil NO (), lightning (), and stratospheric sources () from the year 1986, and set anthropogenic emissions, biomass burning emissions, and CH4 concentrations to specific years. For anthropogenic emissions, we use the Community Emissions Data System (CEDS) inventory (http://globalchange.umd.edu/ceds/) (). Emission species for CEDS include aerosol [black carbon (BC) and organic carbon (OC)], aerosol precursors, and reactive compounds [SO2, NO, NH3, CO, and non-methane volatile organic carbon (NMVOC)]. We use the biomass burning emissions from the CMIP6 (BB4CMIP) inventory for each individual year (). Emission species for BB4CMIP include the following species: BC, CH4, CO, NH3, NMVOC, NO, OC, SO2, and HCl. We prescribe latitudinal CH4 concentrations for historical simulations. For years after 1979, CH4 concentrations are based on the National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA/ESRL) Global Monitoring Division flask observations (http://esrl.noaa.gov/gmd/ccgg/trends_ch4/), and for years before 1979, the CMIP6 monthly mean surface CH4 is used ().

The simulations were performed for years 1960, 1973, 1986, 1999, and 2013, which enables us to isolate the impact of anthropogenic emissions on historical changes in sulfate formation pathways and based on Eq. 1. We considered the error in the modeled by propagating the uncertainties of and , yielding 1σ uncertainty of smaller than 0.1‰ (table S3). In addition to the five model years with calculated cloud pH, we test the same model but assume cloud water pH is constant (pH 4.5) for 1973 and 2013, to examine the importance of changes in bulk cloud pH for modeled over the period.

In addition to the calculation of the transported sulfate in the model, we also extracted online diagnostics to calculate sulfate production efficiency (η) using Eq. 4. S(SO2) was calculated by summing the emission of SO2 from anthropogenic and volcanic activity, photochemical production of SO2 from dimethyl sulfide (DMS) oxidation in the atmosphere, and the net import/export budget of SO2 by transportation. The production of [P()] was also calculated online as the sum from all aqueous, heterogeneous, and gas-phase production pathways.