Sulfur fertiliser use in the Midwestern US increases as atmospheric sulfur deposition declines with improved air quality

Sulfur, as an essential nutrient for plant growth, has increasingly been used in fertiliser applications for many crops. This increase is coincident with declines in atmospheric sulfur deposition in response to air quality improvements in the United States and Europe. Here, we evaluate trends in sulfur fertiliser sales by mass, as a proxy for fertiliser applications, and estimate total atmospheric sulfur deposition across the Midwestern United States. Crop acreage, yield and sulfur fertiliser application substantially increased between 1985 and 2015, coincident with declines in atmospheric sulfur deposition. The increase in sulfur fertiliser has outpaced the relative rate of change in other major nutrient fertilisers including nitrogen, phosphorus and potassium, by approximately 7-fold prior to 2009, and 29-fold after 2009. We suggest that there is a critical need to develop sulfur management tools that optimize fertiliser applications to maintain crop yields while minimizing the consequences of excess sulfur in the environment.


M
anagement of nitrogen (N) and phosphorus (P) fertilisers has long been a concern in agricultural systems around the world, including the United States.Consequently, substantial agronomic research has focused on the timing and amounts of N and P applications to maximize yields while minimizing environmental impacts, including greenhouse gas emissions, soil degradation, and eutrophication of surface waters [1][2][3][4] .Like N and P, sulfur (S) is both a nutrient required for plant growth 5 and can have notable environmental consequences, as demonstrated by ecosystem studies of acid rain in the 1970s and 1980s, including acidification of soils and surface waters, mobilization of toxic metals, and others 6 .Until recently, crop S demand has largely been met by historically high atmospheric S deposition 7 .Yet with air quality changes due to regulation and shifts in energy generation, the free supply of S has declined in the U.S., Europe, and elsewhere (Fig. 1), requiring application of S fertilisers and/or release of legacy S stored in soils-the dynamics and magnitude of which are unknown-to meet plant demand (Fig. 1).This fundamental shift in human manipulation of the S cycle, from diffuse atmospheric inputs as a component of acid rain to targeted applications, necessitates a critical analysis of agricultural S loads, soil S dynamics, and associated environmental consequences 8 .
The Midwestern region of the U.S., which provides approximately one-third of global maize and soybean production 9 , also experienced some of the highest rates of atmospheric S deposition derived from fossil fuel emissions prior to the Clean Air Act and Amendments 10 .On one hand, elevated atmospheric S depositionat its highest rate of ~20 kg S ha −1 yr −1 through the mid-1980smet or exceeded the demand of major crops, like maize [Zea mays (L.)] (Fig. 1).However, increases in both the footprint and yields of maize and soybean [Glycine max (L.) Merr.] in the region have driven increased crop S demand that must be met by sources other than atmospheric deposition 11 (Fig. 2).While there have been many efforts to study and improve the timing and efficiency of N and P applications across the Midwestern U.S. and other regional crop systems [12][13][14][15] , comparable attention has not been given to characterizing the fate or optimization of agricultural S inputs.
One challenge in determining the environmental and human health consequences of large-scale S applications in crop systems has been a lack of publicly available data from individual farmers and governments 8 .While some U.S. states, such as California, have mandatory reporting of S pesticide applications 16 , they are the exception.The U.S. Department of Agriculture (USDA) conducts periodic surveys of S and other fertilisers, but these are neither conducted annually nor are comprehensive across all crops.Here, we present a compiled dataset of S-containing fertilisers from the Association of American Plant Food Control Officials (AAPFCO), available from 1985-2015, and use it as a proxy for S inputs to Midwestern U.S. croplands.While the data are reported by mass sold and at the county level, we provide an analysis of trends across the Midwest region (12 states; Fig. 2a).The results are presented in aggregate, as we recognize that in some cases, fertiliser sales and application location are an imperfect match; large-scale growers may purchase products in one county (reported) and use them in another (unreported).We compare the trends in S fertiliser products, atmospheric S deposition, and estimated S exported in maize tissues using wet S deposition data from the National Atmospheric Deposition Program (NADP) 17 , dry S deposition from the Clean Air Status and Trends Network (CASTNET) 18 , and precipitation data from PRISM 19 (both NADP and CASTNET data are available starting in 1987); fertiliser product sales from AAPFCO; and annual crop acreage and yield data from USDA 20,21 .In addition, we compared the trends in S fertiliser products to those of N, P, and K for the study region (data from AAPFCO).

Results and discussion
From 1987 to 2019, atmospheric S deposition declined at a rate of −0.12 kg S ha −1 yr −1 , from 4.7 to 1.1 kg ha −1 yr −1 , averaged across the Midwestern U.S. region (R 2 = 0.96, p < 0.0001; Fig. 1b).Today, atmospheric S deposition is close to background (i.e., preindustrial) levels 22 .Comparatively, all S-containing fertiliser products have increased at a rate of 0.10 kg S ha −1 yr −1 , from 1.3 to 4.9 kg ha −1 yr −1 between 1985 and 2015, averaged over Midwestern croplands (R 2 = 0.82; p < 0.0001; Fig. 1b).Because S is added to crops not only as a primary fertiliser, but also as a carrier anion for other fertilisers (e.g., N), and the target crop is not reported in the AAPFCO dataset, we also evaluated the change in sulfur-S over the same period.When S was clearly the target Areal loads are calculated as the total amount of S normalized by the total land area of the study region (atmospheric S deposition), the total cropland in the study region (all S fertiliser products), and area of maize planted (maize S).Best fit lines are slope = −0.12kg S ha −1 yr −1 , R 2 = 0.96 (atmospheric S deposition), 0.10 kg S ha −1 yr −1 , R 2 = 0.82 (fertiliser S inputs), and slope = 0.06 kg S ha −1 yr −1 , R 2 = 0.46 (maize S).The linear models are statistically significant (p < 0.001).The slopes of all regression lines are significantly different from zero (p < 0.0001).The mass balance estimate (i.e., atmospheric S deposition and fertiliser inputs less maize S export) is shown for reference.A negative value indicates that atmospheric S deposition plus fertiliser S application exceed losses from maize export.A positive value shows that atmospheric S deposition plus fertiliser S application does not compensate for losses from maize export.
nutrient added to crops (at > 85% S content), the increase over time was also significant (R 2 = 0.71; p < 0.0001; see Data Availability).
The rate of change in average S fertiliser load is important to consider in the context of changes in other major fertilisers over the same period: N, P, and potassium (K).Indeed, the percent change shown as the anomaly relative to 1985 (Fig. 3) illustrates that all S-containing products, as well as sulfur-S, are increasing at a significantly faster rate than other major fertilisers (p < 0.0001).Over the period of record, N and P increased, yet the rates of change in all nutrients other than S were not significantly different from zero (p > 0.05; Fig. 3).Note the rate of change in S calculated with a best-fit linear regression underestimates the marked increase in S loads since 2009.For the period from 2009-2015, the relative change in all S products was 29.6% yr −1 and sulfur-S was 21.1% yr −1 (p < 0.006 and 0.04, respectively).This uptick in use may reflect farmers' recognition of the need to augment S inputs to account for the decline in atmospheric S deposition and insufficient S released from the soil organic pool.
In addition to the decline in atmospheric S deposition, there are three other important factors that likely contribute to the increase in use of S fertiliser products across the Midwestern U.S. First, while total crop acreage has fluctuated from 1985-2020, maize and soybean have increasingly dominated the agricultural footprint of the region and expanded into the Upper Midwest 23,24 ;Fig.2a.Maize yield has increased significantly (p < 0.004; Fig. 2b), creating both higher nutrient export in crop tissues and demand for more intensive fertiliser inputs.Other crops in the region include hay, winter wheat, sorghum, and oats 20 .Sulfur fertiliser is added at notable rates to alfalfa hay, at ~45 kg ha −1 yr −125 .However, by 2015, the combined footprint of hay and other secondary crops was <25% of the total crop acreage 20 , indicating that they are likely to be a minor influence on fertiliser trends in this region.
Yet, productivity of maize and soybean alone cannot explain the regional trends in S fertiliser use.Variability in S fertiliser sales records at finer, local scales may be driven by a second important factor: the distribution of soil S pools and spatial heterogeneity in S cycling and availability to plants (e.g., via mineralization of organic S).Along the Eastern Seaboard of North America, stream export of stored, legacy S has been higher in the northeastern U.S., which was previously glaciated, than in the southeastern U.S., which was not 26 -a pattern controlled by soil properties associated with glaciation.We estimate that soil S pools in the Midwestern U.S. range from 0.002-2.3kg S m −2 (median value of 0.23 kg S m −2 ) in the upper 0.9 m; 27,28 see Supplementary Note 2, suggesting that there is a substantial S reservoir.While our quasi-mass balance of S in crop systems (Fig. 1b) indicates that more recently, inputs may meet crop S demand, this estimate is conservative.It does not account for S export in soybean (often double-cropped with maize) or riverine export of sulfate-S, which occurs at a rate higher than estimated S export in crop tissues 8 .Elevated S losses from the Midwestern region compared with total S inputs highlights a critical need to investigate internal soil S cycling processes that regulate mobility and residence time.Factors like organic matter content, soil type and texture, and rates of mineralization will be an important component of sustainable S management plans.Currently, there does not exist temporal or spatially explicit data on soil S storage or process rates in regional crop systems like the Midwestern U.S (Fig. 1a).7% change in all S products yr −1 , R 2 = 0.82; and slope = 6.6% change in sulfur-S yr −1 , R 2 = 0.71.For the period from 2009-2015, the slopes of the trendlines for all S-containing products and sulfur-S are 29.6% yr −1 (R 2 = 0.91 and p < 0.006) and 21.1% yr −1 (R 2 = 0.55 and p < 0.04), respectively.All linear models are statistically significant (p < 0.001) except K (p = 0.52).The slopes of N, P, and K trends are significantly different from those of all S-containing and sulfur-S products (p < 0.0001); the trend in sulfur-S is not significantly different from that of all S containing products (p = 0.093).Only the slopes of all S products and sulfur-S are significantly different from zero (p < 0.0001).
Finally, decision-making of individual farmers and management companies with respect to treating crop S deficiencies and other practices (e.g., tilling, buffer strips, irrigation, and use of other agriproducts) is another source of variability in regional trends.Our estimate of ~5 kg S ha −1 yr −1 (fertiliser S) is averaged over the entire crop acreage of the region (Fig. 2a).However, S inputs by individual farmers are variable, with far higher applications in some areas than others.In addition, trends in other fertiliser use, such as switching from ammonium nitrate to ammonium sulfate to supply N-a substantial nutrient addition to maize, for example, could affect S loading (see Supplementary Note 1 and Supplementary Fig. 1).Thus, the explanation for the broad trend in increasing S fertiliser use across the Midwestern U.S. is likely a combination of factors: declines in atmospheric S deposition, yield dynamics driven by climate and soil type, local soil S cycling rates, and management decisions.An integrated effort to quantify these components of the current, altered S cycle is an important next step to inform sustainable S management at local to regional scales.
While our multiple lines of evidence suggest that use of targeted, S fertiliser applications was not widespread prior to the Clean Air Act and Amendments, we do not conclude that dirtier air is better for large-scale agriculture.Recent research 29 demonstrates that the benefits of air quality improvements to maize and soybean yields in the Midwest are substantial; namely, reduction in ozone, particulate matter, sulfur dioxide, and nitrogen dioxide.Rather, we argue that with air quality regulation and high agricultural productivity continuing as priorities not only in the U.S., but also in many parts of the world, the pressure to add S fertilisers will continue to increase.This shift in human manipulation of the S cycle-from fossil fuel emissions to agricultural inputs-requires a concerted effort to understand the dynamics of S released from the soil pool, as well as the long-term consequences of S applications for local and adjacent ecosystems 8 .Lessons from studies of excess N and P fertilisers in the environment provide motivation to proactively investigate and address how to sustainably manage S additions in agriculture.Such efforts could yield benefits for major crops, ecosystems, and people around the world.

Methods
Compilation and analysis of fertiliser sales data.We compiled fertiliser sales data from the Association of American Plant Food Control Officials (AAPFCO) for all years currently available (1985-2015) 30 .Sales data (in pounds of each product purchased) are reported annually by farmers and aggregated at the county level.We separated all fertilisers into groups containing nitrogen (N), phosphorus (P), potassium (K), and sulfur (S), regardless of whether each nutrient was the target addition (e.g., N in ammonium sulfate) or a carrier (e.g., S in ammonium sulfate); we used this approach to calculate a total load of each nutrient for our study region.In some cases, the percentage N, P, K, and/or S was provided, while in others it was not.In the absence of reported nutrient content data, we calculated a range of possible values (low, average, and high) based on publicly available product information and/or published literature (see Data Availability).By multiplying the elemental content of a product by its total product mass sold per county, we calculated a mass of N, P, K, and S for each product per county per reporting period.It is possible that some products containing trace amounts of N, P, K, and S were not included in our analysis.However, these trace amounts would be a small amount of the total load, which is normalized over all cropland in the Midwestern U.S. region.
Ultimately, products for which N, P, K, and/or S content was not reported in the AAPFCO dataset did not have appreciable differences in the calculated loads across the range of possible stoichiometric values and the trends were the same.Thus, we reported the loads based on the average reported content of each element in each fertiliser product.In addition, we found that the county-level data reporting was imperfect in some instances.For example, sometimes the data were reported for a state, but the county (FIPS) identifier was absent.We also recognized that large producers might purchase fertilisers in a different location (i.e., county) from where it was applied.Thus, we aggregated the data to evaluate trends across a 12-state region in the Midwest encompassing much of the maize and soybean cultivated in the U.S.
Note that in cases when more fertiliser product was purchased than was used in a particular location over the course of the year, the unused mass was subtracted from the value in the following year.This resulted in reporting of some negative values year to year, reflecting a more accurate accounting of fertiliser sales over the entire period.
To estimate area-normalized S loads, we divided the total weight of each fertiliser (N, P, K, S) by the total crop acreage data reported for each state by year from the USDA Economics, Statistics, and Market Information System 20 .The total crop acreage includes area planted for maize, sorghum, oats, barley, rye, winter wheat, Durum wheat, other spring wheat, rice, soybean, peanuts, sunflower, cotton, dry edible beans, potatoes, sugar beets, canola, and proso millet.Hay, tobacco, and sugarcane are included as harvested acreage.These totals include double-cropped acres and unharvested small grains planted as cover crops.Again, crops other than maize and soybean were the minority (<25%) of the total acreage in the study region.To estimate S exported in maize tissues, we used tissue data reported by University of Nebraska-Lincoln Extension 31 and maize acres planted from the USDA NASS QuickStats database 21 .It is important to recognize that this estimate is conservative; we do not have a time series of S content of maize tissues, nor do we include soybean, as we could not parse areas single-or double-cropped with maize.Finally, we estimated yield trends in soybean, maize (grain), and maize (silage) for the Midwest region using a combination of yield (quantity per acre) and acres harvested reported in the USDA NASS QuickStats database 21 .
Information on estimating soil S pools is provided in Supplementary Note 2.
Estimation of atmospheric S deposition fluxes.We estimated total atmospheric S deposition for the study region using annual volume-weighted sulfate concentrations in wet-only deposition measurements from the National Atmospheric Deposition Program (NADP) 17 , estimates of dry S deposition from the U.S. Environmental Protection Agency Clean Air Status and Trends Network (CASTNET) 18 and precipitation quantity data from the PRISM spatial climate datasets 19 .We interpolated total deposition for unmonitored regions using estimates at point locations.This analysis was accomplished using a spatial model that incorporates precipitation quantity, annual volume-weighted mean S concentrations in precipitation, and the dry deposition data for particulate sulfate and sulfur dioxide.The model then uses a Kriging approach to determine the spatial pattern of S concentration in precipitation from the network of NADP stations within the 12 Midwestern U.S. states included in this analysis.Similarly, we used Kriging to generate spatial patterns of dry S deposition using point data obtained from 14 sites monitored as part of the CASTNET program.The annual total S deposition was generated for each of 12 states for all individual years between 1989-the year when adequate point data were available through the networks to develop the Kriging models-and 2017.The overall uncertainty in flux estimates is comprised of multiple sources of uncertainty associated with the components of annual total S deposition.For wet S deposition there is uncertainty in the weekly measurements of precipitation volume, as well as sulfate concentrations; weekly measurements are summed to give annual fluxes of wet sulfate deposition.The quality assurance procedures for the NADP are summarized by the laboratory 17 .For dry S deposition, there are uncertainties associated with the measurements of gaseous sulfur dioxide and particulate sulfate concentrations, as well as the modeled deposition velocity values to estimate deposition flux.Finally, there is uncertainty associated with the spatial extrapolation of point measurements from the deposition networks to the state scales.We estimated areal-normalized atmospheric S deposition using state areas for the Midwestern U.S. study region.
Statistical analyses.To evaluate the trends in atmospheric S deposition, maize and soybean yields, and fertiliser sales over time, we used the lm function in R 32 to determine the best-fit linear regression lines.We report slopes of the best-fit lines, as well as their R 2 and p-values.We also test the hypotheses that (1) trends in N, P, and K products are statistically different from the trend in S products (all S-containing products, and sulfur-S products only, p < 0.05) and ( 2) that the slopes of all regression lines are statistically different from zero (p < 0.05).

Fig. 1
Fig. 1 Mass balance estimate of S in agriculture from 1985-2017.a Conceptual model of S flows.b Trends in S inputs and estimated S output in maize.Areal loads are calculated as the total amount of S normalized by the total land area of the study region (atmospheric S deposition), the total cropland in the study region (all S fertiliser products), and area of maize planted (maize S).Best fit lines are slope = −0.12kg S ha −1 yr −1 , R 2 = 0.96 (atmospheric S deposition), 0.10 kg S ha −1 yr −1 , R 2 = 0.82 (fertiliser S inputs), and slope = 0.06 kg S ha −1 yr −1 , R 2 = 0.46 (maize S).The linear models are statistically significant (p < 0.001).The slopes of all regression lines are significantly different from zero (p < 0.0001).The mass balance estimate (i.e., atmospheric S deposition and fertiliser inputs less maize S export) is shown for reference.A negative value indicates that atmospheric S deposition plus fertiliser S application exceed losses from maize export.A positive value shows that atmospheric S deposition plus fertiliser S application does not compensate for losses from maize export.

Fig. 3
Fig. 3 Percent change in N, P, K, and S products relative to 1985.The total load of each element was calculated based on the stoichiometry of the fertiliser products and the mass of each product sold.Best fit lines are slope = 0.7% change in N yr −1 , R 2 = 0.52; slope = 1.2% change in P yr −1 , R 2 = 0.62; slope = 0.1% change in K yr −1 , R 2 = 0.01; slope = 7.7% change in all S products yr −1 , R 2 = 0.82; and slope = 6.6% change in sulfur-S yr −1 , R 2 = 0.71.For the period from 2009-2015, the slopes of the trendlines for all S-containing products and sulfur-S are 29.6% yr −1 (R 2 = 0.91 and p < 0.006) and 21.1% yr −1 (R 2 = 0.55 and p < 0.04), respectively.All linear models are statistically significant (p < 0.001) except K (p = 0.52).The slopes of N, P, and K trends are significantly different from those of all S-containing and sulfur-S products (p < 0.0001); the trend in sulfur-S is not significantly different from that of all S containing products (p = 0.093).Only the slopes of all S products and sulfur-S are significantly different from zero (p < 0.0001).