Articles | Volume 20, issue 3
Research article
 | Highlight paper
05 Feb 2020
Research article | Highlight paper |  | 05 Feb 2020

Mapping the drivers of uncertainty in atmospheric selenium deposition with global sensitivity analysis

Aryeh Feinberg, Moustapha Maliki, Andrea Stenke, Bruno Sudret, Thomas Peter, and Lenny H. E. Winkel

An estimated 0.5–1 billion people globally have inadequate intakes of selenium (Se), due to a lack of bioavailable Se in agricultural soils. Deposition from the atmosphere, especially through precipitation, is an important source of Se to soils. However, very little is known about the atmospheric cycling of Se. It has therefore been difficult to predict how far Se travels in the atmosphere and where it deposits. To answer these questions, we have built the first global atmospheric Se model by implementing Se chemistry in an aerosol–chemistry–climate model, SOCOL-AER (modeling tools for studies of SOlar Climate Ozone Links – aerosol). In the model, we include information from the literature about the emissions, speciation, and chemical transformation of atmospheric Se. Natural processes and anthropogenic activities emit volatile Se compounds, which oxidize quickly and partition to the particulate phase. Our model tracks the transport and deposition of Se in seven gas-phase species and 41 aerosol tracers. However, there are large uncertainties associated with many of the model's input parameters. In order to identify which model uncertainties are the most important for understanding the atmospheric Se cycle, we conducted a global sensitivity analysis with 34 input parameters related to Se chemistry, Se emissions, and the interaction of Se with aerosols. In the first bottom-up estimate of its kind, we have calculated a median global atmospheric lifetime of 4.4 d (days), ranging from 2.9 to 6.4 d (2nd–98th percentile range) given the uncertainties of the input parameters. The uncertainty in the Se lifetime is mainly driven by the uncertainty in the carbonyl selenide (OCSe) oxidation rate and the lack of tropospheric aerosol species other than sulfate aerosols in SOCOL-AER. In contrast to uncertainties in Se lifetime, the uncertainty in deposition flux maps are governed by Se emission factors, with all four Se sources (volcanic, marine biosphere, terrestrial biosphere, and anthropogenic emissions) contributing equally to the uncertainty in deposition over agricultural areas. We evaluated the simulated Se wet deposition fluxes from SOCOL-AER with a compiled database of rainwater Se measurements, since wet deposition contributes around 80 % of total Se deposition. Despite difficulties in comparing a global, coarse-resolution model with local measurements from a range of time periods, past Se wet deposition measurements are within the range of the model's 2nd–98th percentiles at 79 % of background sites. This agreement validates the application of the SOCOL-AER model to identifying regions which are at risk of low atmospheric Se inputs. In order to constrain the uncertainty in Se deposition fluxes over agricultural soils, we should prioritize field campaigns measuring Se emissions, rather than laboratory measurements of Se rate constants.

1 Introduction

Selenium (Se) is an essential dietary trace element for humans and animals, with the recommended intake ranging between 30 and 900 µg d−1 (Fairweather-Tait et al.2011). The amount of Se in crops depends on the amount of bio-available Se in the soils where the crops are grown (Winkel et al.2015). Levels of Se in soils, as well as Se dietary intakes, vary strongly around the world (Jones et al.2017). Selenium deficiency is considered a more widespread issue than Se toxicity, as around 0.5 to 1 billion people are estimated to have insufficient Se intake (Combs2001; Fordyce2013).

Atmospheric deposition is an important source of Se to soils. In several regions, Se concentrations in soils were found to correlate with precipitation (Låg and Steinnes1978; Blazina et al.2014). Further, several studies have attributed an increase in soil Se concentrations to regional anthropogenic Se emissions to the atmosphere (Haygarth et al.1993; Dinh et al.2018), suggesting a link between atmospheric Se inputs and soil concentrations. However, apart from some budget studies in the 1980s (Ross1985; Mosher and Duce1987), there has been a lack of research into global-scale atmospheric cycling of Se and the spatial patterns of Se deposition.

1.1 Atmospheric Se cycle

Since Se and sulfur (S) are in the same group on the periodic table, they share chemical properties and their biogeochemical cycles are similar. Like S, Se is emitted to the atmosphere by both natural and anthropogenic sources, with the total annual emissions estimated between 13 and 19 Gg Se yr−1. Natural sources of atmospheric Se include volatilization by the marine and terrestrial biospheres, volcanic degassing and eruptions, and minor contributions from sea salt and mineral dust. Anthropogenic Se is emitted during coal and oil combustion, metal smelting, and biomass burning. However, very few in situ measurements of Se emissions fluxes and speciation exist. Once in the atmosphere, volatile Se species are oxidized, eventually forming species like elemental Se and SeO2 (Wen and Carignan2007). These oxidized species are expected to partition to the particulate phase; previous measurements have found that 75 %–95 % of Se is in particulates (Mosher and Duce1983, 1987). The fate of atmospheric Se is dry and wet deposition, with wet deposition accounting for an estimated 80 % of total deposition globally (Wen and Carignan2007).

Atmospheric chemistry modeling studies have been applied to other trace elements to predict atmospheric lifetimes and spatial patterns of deposition. For example, atmospheric mercury models were developed more than 2 decades ago (Petersen and Munthe1995), and now there are around 16 global and regional atmospheric mercury models (Ariya et al.2015). A recent mechanistic modeling paper has advanced the understanding of atmospheric arsenic cycling (Wai et al.2016). To our knowledge, Se chemistry has never previously been included in an atmospheric chemistry–climate model (CCM), and thus many questions surrounding atmospheric Se transport remain unanswered. For example, the atmospheric lifetime of Se and thus the scales at which it can be transported (local, regional, hemispheric, global) are not well constrained.

1.2 Global sensitivity analysis

Since the atmospheric Se cycle has been investigated only by a limited number of studies, it is essential that we consider the relevant parametric uncertainties when building an atmospheric Se model. The reaction rate coefficients of Se compounds have either only been measured by one laboratory study or no laboratory measurement exists, and these rate coefficients need to be estimated. Selenium emission fluxes from certain sources have been measured; however it remains difficult to extrapolate these measurements to global fluxes due to the high degree of spatial and temporal variability. Atmospheric Se modeling can only be considered trustworthy when combined with full accounting of input parameter uncertainties and their propagation through the model. Through “global sensitivity analysis” (Saltelli et al.2008) we can identify which input uncertainties are the most important for the uncertainty in the model output. A sensitivity analysis is called global when the sensitivity is evaluated over the entire input parameter space, as opposed to local methods that test sensitivity only at a certain reference point in the space (i.e., based on the gradient of the output at this reference point). Sensitivity analysis provides a framework to prioritize which model inputs should be further constrained in order to reduce the uncertainty in the model output.

Until recently, most sensitivity analyses of atmospheric chemistry models consisted of local methods, principally the one-at-a-time (OAT) approach. In the OAT approach, the model is initially run with a set of default parameters to yield a “reference” simulation. Multiple sensitivity simulations are then conducted, so that for each simulation one parameter is perturbed from the reference set at a time. The influence of these perturbations on the model output of interest would then be analyzed. However, this approach may be flawed because it only considers the first-order response of the model to each parameter, ignoring interactions that might exist between parameters (Saltelli et al.2008; Lee et al.2011). Furthermore, the uncertainty ranges of the input parameters are rarely quantified and reported; much of the possible parameter space often remains unexplored. Global methods such as variance-based sensitivity analysis allow the uncertainty in model output to be apportioned to each input variable. Sobol' indices, which represent the fraction of model variance that one input variable explains, provide a ranking system for the importance of input variables (Sobol1993). The benefits of global sensitivity analysis include (1) identifying the most influential input variables, i.e., the ones that should be further constrained to yield the biggest reduction in model uncertainty; (2) identifying input variables that do not play any role in the model output variance; this could represent a route to simplify the model, since the process involving these input variables can be neglected; (3) understanding the behavior of the model, for example how the output depends on interactions between input variables; and (4) identifying possible model bugs or discontinuities, since the model will be tested with a wide range of input parameter values (Saltelli et al.2000; Ryan et al.2018).

There are several recent examples of atmospheric chemistry studies that included a global sensitivity analysis. The sensitivity of cloud condensation nuclei number density to input parameters in an aerosol model was investigated at the local (geographic) scale (Lee et al.2011) and at the global scale (Lee et al.2012, 2013), revealing regional differences in parameter rankings. Revell et al. (2018) investigated the sensitivity of the tropospheric ozone columns to emission and chemical parameters, to identify which processes are responsible for the bias in modeled tropospheric ozone. Marshall et al. (2019) employed global sensitivity analysis methods to identify how radiative forcing responds to volcanic emission parameters. In these examples, surrogate modeling techniques (also known as emulation) were employed to replace a process-oriented, computationally expensive model with an approximative statistical model. The statistical model has the advantage that it is quicker to evaluate; therefore, it can be used to calculate the model output throughout the parametric space (Lee et al.2011). The examples given above all used Gaussian process emulation; however other surrogate modeling techniques exist, including polynomial chaos expansions (PCEs) (Ghanem and Spanos2003). The PCE approach is well suited to sensitivity analysis, since the Sobol' sensitivity indices can be extracted analytically from the constructed PCE, with no need to evaluate the surrogate model through Monte Carlo sampling (Sudret2008). This can greatly reduce the computational time required to conduct the sensitivity analysis, especially when one is interested in conducting a separate sensitivity analysis for many model grid boxes.

1.3 Outline

This study focuses on the construction of the first global atmospheric Se model and the insights that this model reveals into which Se cycle uncertainties would be the most important to constrain. Section 2 describes the SOCOL-AER (modeling tools for studies of SOlar Climate Ozone Links – aerosol) model and the implementation of Se chemistry in the SOCOL-AER model. The SOCOL-AER model is a suitable tool to model the Se cycle, since it successfully describes the major properties of the atmospheric S cycle (Sheng et al.2015; Feinberg et al.2019b). The statistical methods that we use to conduct the sensitivity analysis are discussed in Sect. 3. Section 4 details the methods used to compile a database of measured wet Se deposition fluxes, which we use to evaluate the model. The results of the sensitivity analyses are presented in Sect. 5.1 for the atmospheric Se lifetime and Sect. 5.2 for the Se deposition patterns. Section 5.3 illustrates the comparison between the compiled Se deposition measurements and simulated results. A discussion of both sensitivity analyses follows in Sect. 6, and the paper is concluded in Sect. 7.

2 Model description


SOCOL-AERv2 is a global CCM that includes a sulfate aerosol microphysical scheme (Sheng et al.2015; Feinberg et al.2019b). The base CCM, SOCOLv3 (Stenke et al.2013), is a combination of the general circulation model ECHAM5 (Roeckner et al.2003) and the chemical model MEZON (Egorova et al.2003). The MEZON submodel comprises a comprehensive atmospheric chemistry scheme, with 89 gas-phase chemical species, 60 photolysis reactions, 239 gas-phase reactions, and 16 heterogeneous reactions. Chemical tracers are advected in the model using the Lin and Rood (1996) semi-Lagrangian method. Photolysis rates are calculated using a look-up table approach based on the simulated overhead ozone and oxygen columns. The MEZON model solves the system of differential equations representing chemical reactions with a Newton–Raphson iterative method for short-lived chemical species and an Euler method for long-lived species.

The sulfate aerosol model AER (Weisenstein et al.1997) was first coupled to SOCOL by Sheng et al. (2015). SOCOL-AER includes gas-phase S chemistry and 40 sulfate aerosol tracers, ranging in dry radius size from 0.39 nm to 3.2 µm. SOCOL-AER simulates microphysical processes that affect the aerosol size distribution, including binary homogeneous nucleation, condensation of H2SO4 and H2O, coagulation, evaporation, and sedimentation. SOCOL-AER was extended in Feinberg et al. (2019b) to include interactive wet and dry deposition schemes. The wet deposition scheme, based on Tost et al. (2006), calculates scavenging of gas-phase species depending on their Henry's law coefficients and aerosol species depending on the particle diameter. The wet removal of tracers is coupled to the grid cell simulated properties of clouds and precipitation. The dry deposition scheme is based on Kerkweg et al. (2006, 2009), who use the surface resistance approach of Wesely (1989). In addition to surface type and meteorology, the calculated dry deposition velocities depend on reactivity and solubility for gas-phase compounds and size for aerosol species. SOCOL-AER uses an operator splitting approach, wherein the model time step is 2 h for chemistry and radiation and 15 min for dynamics and deposition. Aerosol microphysical routines use a sub-time step of 6 min.

For the simulations in this study we use boundary conditions for the year 2000. Sea ice coverage and sea surface temperatures are prescribed from the Hadley Centre dataset (Rayner et al.2003). The year 2000 concentrations of the most relevant greenhouse gases (CO2, CH4, and N2O) derive from NOAA observations (Eyring et al.2008). Anthropogenic CO and NOx emissions are based on the RETRO dataset (Schultz and Rast2007), while natural emissions are taken from Horowitz et al. (2003). Sulfur dioxide (SO2) emissions from anthropogenic sources follow the year 2000 inventory from Lamarque et al. (2010) and Smith et al. (2011). Volcanic degassing SO2 emissions are assigned to surface grid boxes where volcanoes are located (Andres and Kasgnoc1998; Dentener et al.2006). Atmospheric emissions of dimethyl sulfide (DMS) are calculated using a wind-based parametrization (Nightingale et al.2000) and a marine DMS climatology (Lana et al.2011). To represent mean conditions for photolysis, the look-up table for photolysis rates is averaged over two solar cycles (1977–1998).

Figure 1Scheme of the atmospheric Se cycle in SOCOL-AER, based on information in Ross (1985) and Wen and Carignan (2007). For simplicity the two oxidized Se tracers in SOCOL-AER are represented with a single box. The impact on agriculture and human health is also shown, since it motivates the study of atmospheric Se.

2.2 Implementing Se emissions and chemistry in SOCOL-AER

2.2.1 Selenium species overview

We included the Se cycle in SOCOL-AER (Fig. 1) based on the existing literature on atmospheric Se (Ross1985; Wen and Carignan2007). Seven Se gas-phase tracers have been added to SOCOL-AER (Table 1): carbonyl selenide (OCSe), thiocarbonyl selenide (CSSe), carbon diselenide (CSe2), dimethyl selenide (DMSe), hydrogen selenide (H2Se), oxidized inorganic Se (OX_Se_IN), and oxidized organic Se (OX_Se_OR). The oxidized inorganic Se tracer represents species such as elemental Se, selenium dioxide (SeO2), selenous acid (H2SeO3), and selenic acid (H2SeO4). Very little is known about the kinetics of interconversion between the oxidized inorganic Se species (Wen and Carignan2007), and therefore in our model they are treated as one tracer. However, these species all have very low vapor pressures under atmospheric conditions (Rumble2017) and likely partition to the particulate phase. Oxidized organic Se species include dimethyl selenoxide (DMSeO) and methylseleninic acid (CH3SeO2H), which form after oxidation of DMSe (Atkinson et al.1990; Rael et al.1996). Similar to the oxidized inorganic Se compounds, oxidized organic Se species also partition to the particulate phase due to their low volatilities (Rael et al.1996).

2.2.2 Selenium emissions

To determine which Se compounds are emitted by the different sources, we have reviewed studies that investigated the speciation of Se emissions. Thermodynamic modeling and in situ measurements of combustion exhaust gases have detected the following Se species in anthropogenic emissions: oxidized inorganic Se, H2Se, OCSe, CSe2, and CSSe (Yan et al.2000, 2004; Pavageau et al.2002, 2004). Oxidized inorganic Se species and minor amounts of H2Se are expected to be emitted from volcanic degassing (Symonds and Reed1993; Wen and Carignan2007). A variety of methylated Se species have been observed from biogenic emissions, including DMSe, dimethyl diselenide (DMDSe), dimethyl selenylsulfide (DMSeS), and methane selenol (MeSeH) (Amouroux and Donard1996, 1997; Amouroux et al.2001; Wen and Carignan2007). Since DMSe is usually the dominant emitted compound, and little is known about the oxidation kinetics of the other methylated species, DMSe is the only species emitted by marine and terrestrial biogenic emissions in SOCOL-AER.

Table 1Description of Se tracers included in SOCOL-AER.

Download Print Version | Download XLSX

Atmospheric S emissions have been measured more extensively than Se emissions, so we scale inventories of S emissions to yield the spatial distribution of emitted Se. For the sensitivity analysis we assume that the Se emissions have the same distribution as S emissions, and we focus on the uncertainties in the global scaling factors for each source. The range in scaling factors will be discussed in Sect. 3.1.4. The spatial distribution of anthropogenic and biomass burning Se emissions comes from the SO2 inventory for the year 2000 (Lamarque et al.2010; Smith et al.2011). For sea surface DMSe concentrations, we scale the DMS climatology calculated by Lana et al. (2011). We determine DMSe emissions using a wind-driven parametrization (Nightingale et al.2000). The locations and strength of background degassing volcanic emissions are taken from the Global Emissions InitiAtive (GEIA) inventory (Andres and Kasgnoc1998; Dentener et al.2006). Since little is known about both terrestrial biogenic Se emissions and terrestrial S emissions (Pham et al.1995), we assume that terrestrial Se is emitted in all land surface grid boxes, excluding glaciated locations.

2.2.3 Chemistry of Se species

We conducted a literature review to develop the model's chemical scheme of the Se cycle. Reactions of Se species that have been measured by laboratory studies are compiled in Table 2. We neglect any temperature dependency in the Se reaction rates, since the Se reactions have only been studied at around 298 K. For all compiled reactions, atmospheric Se compounds react much quicker than the analogous S compounds, due to the stronger bonds that S forms with carbon and hydrogen than Se (Rumble2017). In addition, there are reactions that are known to occur for the analogous S compound, but have never been studied for the Se compound (OCSe+OH, CSSe+OH, CSe2+OH, and H2Se+OH). These reaction rate constants were estimated in Fig. 2, which shows the ratio of the Se compounds’ rates to analogous S rates (i.e., an enhancement factor for replacing an S atom with Se) plotted versus the S reaction rate. For S reactions which have a fast reaction rate (e.g., DMS+Cl, k=1.8×10-10cm3molec.-1s-1), replacing S with Se does not yield a large difference in measured rates (DMSe+Cl, k=5.0×10-10cm3molec.-1s-1). This is because these reactions are already close to the collision-controlled limit, and thus lowering the activation energy by substituting a Se atom for S has little impact on the overall rate. On the other hand, slow reactions like DMS+O3 are sped up by more than 4 orders of magnitude when Se is substituted for S. We used the log-log relationship in Fig. 2 to predict the reaction rates for OCSe, CSSe, and H2Se with OH (blue stars). The CSe2 reaction with OH is calculated from the CSSe reaction, assuming an enhancement for the substitution of a second Se atom similar to that between the measured CSe2+O and CSSe+O reaction rates (Li et al.2005). The branching ratio for the CSSe+OH reaction products was assumed to be 30 % OCSe and 70 % OX_Se_IN, the same as the measured CSSe+O branching ratio (Li et al.2005). We recognize that these estimates are inherently uncertain, and therefore we address these uncertainties in our sensitivity analysis (Sect. 3.1.2).

Li et al. (2005)Li et al. (2005)Li et al. (2005)Atkinson et al. (1990)Atkinson et al. (1990)Atkinson et al. (1990)Thompson et al. (2002)Agrawalla and Setser (1987)Agrawalla and Setser (1986)Belyaev et al. (2012)

Table 2Rate constants of Se compound gas-phase reactions at around 298 K and the corresponding rate constant of the analogous S compound. All S reaction rates are from Burkholder et al. (2015), except the DMS+O3 reaction rate which is taken from Wang et al. (2007). No corresponding rate constants for CSe2 reactions are listed, since CSe2 is obtained from doubly substituting Se in CS2.

Download Print Version | Download XLSX

The photolysis of gas-phase Se compounds was included using absorption cross sections of H2Se (Goodeve and Stein1931) and OCSe (Finn and King1975). The absorption cross section of CSe2 (King and Srikameswaran1969) has been measured, however in too low resolution to be incorporated into the model. Therefore, we assume that CSe2 and CSSe have the same cross section as CS2 (Burkholder et al.2015). Given the lack of available information, quantum yields for all Se photolysis reactions were assumed to be 1.

2.2.4 Condensation of Se on preexisting aerosol particles

As nonvolatile species, oxidized inorganic and organic Se would condense on available atmospheric surfaces. In the SOCOL-AER model, the uptake of these oxidized Se species by sulfate aerosols is calculated similarly to the existing scheme of gas-phase H2SO4 uptake on sulfate particles (Sheng et al.2015). We track the size distribution of Se in the aerosol phase with 40 tracers, one for each sulfate aerosol size bin. The sulfate aerosol size distribution changes through processes like growth, evaporation, and coagulation. We track how these microphysical processes change the size distribution of condensed Se through mass-conserving schemes. Evaporation of condensed Se only occurs when the smallest sulfate aerosol bin evaporates, releasing the Se stored in that bin as gas-phase inorganic oxidized Se. Sedimentation and deposition of the host sulfate particles are concurrent with sedimentation and deposition of the condensed Se tracers. Gas-phase oxidized Se tracers are also removed by dry and wet deposition, with the assumption that they have the same Henry's law constant as gas-phase H2SO4.

One limitation of using SOCOL-AER for the Se cycle is that the model only includes online sulfate aerosols. This means that the transport of Se on other aerosols, including dust, sea salt, and organic aerosols, would be neglected. This may not be a poor assumption, since Se and S are often co-emitted and have been found to be highly correlated in atmospheric aerosol measurements (e.g., Eldred1997; Weller et al.2008). Nevertheless, we included a “dummy” aerosol tracer to test the effect of missing aerosol species in SOCOL-AER. The dummy aerosol tracer represents monodisperse particles that are emitted in a latitudinal band in the model and undergo Se uptake, sedimentation, and wet and dry deposition. This dummy aerosol tracer is clearly a simplification of true atmospheric processes, as in reality other aerosols are distributed in size and can coagulate with sulfate aerosol particles. However, by varying the radius, location, and emission magnitude of these particles (Sect. 3.1.7), we can determine whether missing aerosols affect atmospheric Se cycling.

Figure 2Estimation of unknown Se reaction rates from the analogous S reaction rate. A power regression is performed, with statistics shown in the upper right corner of the plot. For the H2S+O3 reaction only an upper limit estimate was available, and therefore it was not included in the analysis.


3 Statistical methods

To conduct the sensitivity analysis of our Se model, we first need to select the input parameters that would be included in the sensitivity analysis. The probability distributions of these input parameters' uncertainties were determined by reviewing literature sources and using our best judgment. Variance-based sensitivity analysis methods usually require 104 to 106 model runs, which would be prohibitively expensive for the full SOCOL-AER model. We therefore replace the SOCOL-AER model with a surrogate PCE model. The SOCOL-AER model is run at carefully selected points within the parameter space, creating a set of “training” runs. The training runs are used to produce a surrogate PCE model, which approximates the outputs of the full SOCOL-AER model throughout the input parameter space. Sensitivity indices can then be derived from the surrogate model. All statistical methods presented in this section are available in UQLab, an open-source MATLAB-based framework for uncertainty quantification (Marelli and Sudret2014).

3.1 Uncertainty ranges of input parameters

We restricted the scope of our sensitivity analysis to parameters that have been implemented in the model as part of the Se cycle. We neglect other model parameters, including those related to climate, deposition parametrizations, or the emissions of sulfate aerosol precursors. The focus of our sensitivity analysis is to prioritize which Se-related parameters should be further constrained. Since we do not vary all other model parameters, the uncertainties of output quantities may be underestimated. However, given the large dimension of our parameter space with 34 Se-related input parameters, including additional non-Se-related parameters would be challenging. In the following section, we will discuss the uncertainty distributions for each of the 34 input parameters included in our study. Due to the lack of detailed information available in literature about the parameter distributions, we chose loguniform or uniform distributions for all but one of the parameters. This follows the conservative approach recommended by the maximum entropy principle, as the uniform and loguniform distributions maximize entropy while fulfilling the data constraints (Kapur1989). The uncertainty distributions of all input parameters are listed in Table 3.

Table 3Probability distributions of the model input parameters selected for the sensitivity analysis.

Download Print Version | Download XLSX

3.1.1 Measured rate constants (k1k12)

The Se reactions studied in the literature each have only been measured by one laboratory group (Table 2). Since only one measurement technique has been applied, the reported measurement uncertainties may underestimate the true uncertainties of these rate constants. To approximate an uncertainty distribution for these rate constants, we reviewed S compound reactions that have been studied by multiple research groups. The reaction that had the largest spread in reported rate constants at ∼298 K was OCS+OH, which has been measured in six studies (Atkinson et al.1978; Kurylo1978; Cox and Sheppard1980; Leu and Smith1981; Cheng and Lee1986; Wahner and Ravishankara1987; Burkholder et al.2015). The measured reaction rate constant varied over multiple orders of magnitude; therefore, we calculated its variability on a logarithmic scale. The coefficient of variation (standard deviation divided by mean) of this reaction rate in logarithmic space was around 6 %. We assumed that this maximum S coefficient of variation would apply to the measured Se reaction rates. The bounds were calculated as 88 % and 112 % (±2 standard deviations) of the available measured rate constant in logarithmic space, i.e.,

(1) bounds = k 1 ± 0.12 ,

where k is the measured rate constant and “bounds” are its upper and lower bounds, all expressed in cubic centimeters per molecule per second. The maximum upper bound was set to 1.0×10-9cm3molec.-1s-1, since at this order of magnitude the Se reaction rates are collision-limited (Seinfeld and Pandis2016).

3.1.2 Estimated rate constants (k13k17)

Five Se reaction rate constants have not been measured previously in the laboratory and were estimated based on their relationship with analogous S rate constants (Fig. 2). We calculated the uncertainty bounds of estimated rate constants using the 95 % error interval of prediction with a linear regression (Wackerly et al.2014):

(2) bounds = 10 x + Y ^ ± t 0.025 SSE n - 2 1 + 1 n + ( x - x ) 2 S x x ,

where x is the logarithm (to the base 10) of the corresponding S rate constant, Y^ is the logarithm of the predicted ratio of the Se rate to the S rate, n is the number of data points in the regression, t0.025 is the 2.5th percentile value of the Student t distribution for n−2 degrees of freedom, x is the mean of the logarithm of S rate constants in Fig. 2, SSE is the sum of squares of the residuals, and Sxx is the variance of the S rate constants in Fig. 2. All rate constants are in units of cubic centimeters per molecule per second.

3.1.3 Photolysis rate scaling factors

Uncertainties in our calculated Se photolysis rates arise from uncertainties in the measured OCSe and H2Se cross sections, the assumption that CSe2 and CSSe have the same cross section as CS2, the quantum yields of photolysis reactions, and the look-up table approach that SOCOL-AER applies to calculate photolysis rates. Given the lack of specific information about these uncertainties, we apply a scaling factor on the calculated photolysis rates ranging from 0 to 2 (Table 3).

3.1.4 Global emissions from source categories

For the sensitivity analysis, we do not alter the spatial distribution of Se emissions from anthropogenic, marine biogenic, terrestrial biogenic, and volcanic sources (Fig. 3). The parameters that we varied are the scaling factors for each map, i.e., the global total emissions from each source. We reviewed atmospheric Se budget estimates to determine the range in total emissions for different sources (Table S1). The estimates for global DMSe emissions ranged from the lower limit value of 0.4 Gg Se yr−1 in Nriagu (1989) to 35 Gg Se yr−1 in Amouroux et al. (2001). DMSe emissions are calculated online in the model from a marine DMSe concentration map. Based on the results of a previous simulation, we normalize the marine DMSe concentration map in Fig. 3a so that it leads to 1 Gg Se yr−1 emissions globally. All other maps are also normalized to 1 Gg Se yr−1 emissions, so that we can directly apply the total global emissions as a scaling factor. The widest uncertainty range of terrestrial emissions was given by Nriagu (1989), from 0.15 to 5.25 Gg Se yr−1. Total global anthropogenic Se emissions in 1983 were estimated between 3.0 and 9.6 Gg Se yr−1 (Nriagu and Pacyna1988). We applied the same uncertainty range to the total anthropogenic emissions in the year 2000, because it is unclear how global Se emissions have changed during this period. To estimate global Se emissions from degassing volcanoes, we reviewed measurements of Se-to-S ratios in volcanic emissions, extending the studies reviewed in Floor and Román-Ross (2012) (Table S2 in the Supplement). There is a high degree of variability in the emitted Se-to-S ratios between different volcanoes, and even temporally and spatially for the same volcano (Floor and Román-Ross2012). The Se-to-S ratio in volcanic emissions ranges from 6×10-6 for White Island, New Zealand (Wardell et al.2008), to 3.9×10-3 for Merapi Volcano, Indonesia (Symonds et al.1987). By multiplying this range in ratios with the global mean total degassing SO2 emissions from Andres and Kasgnoc (1998) and Dentener et al. (2006), 12.6 Tg S yr−1, we calculate an uncertainty range for total volcanic Se emissions: 0.076–49.1 Gg Se yr−1. Loguniform distributions were used for the source types whose total emission ranges vary by more than 2 orders of magnitude (volcanic and marine biogenic), whereas uniform distributions were used for the terrestrial and anthropogenic Se emissions.

Figure 3Spatial distribution of Se sources. Each source map is normalized so that 1 Gg Se yr−1 would be emitted globally. This is scaled during the sensitivity analysis.

3.1.5 Speciation of Se emission sources

The available speciation information for Se emissions is largely qualitative; possible emission species have been identified but not quantified (Sect. 2.2.2). To estimate the uncertainty range for Se emission speciation, we use estimates of speciation from an atmospheric S budget (Watts2000). The second most important anthropogenic S species after SO2 is H2S. Watts (2000) estimates the anthropogenic emission of H2S is 3.1±0.3Tg S yr−1, compared to an anthropogenic SO2 emission total of 53.2 Tg S yr−1 (Lamarque et al.2010; Smith et al.2011) in the year 2000. H2S therefore contributes at most 6 % of total S emissions. Since OCSe, CSe2, CSSe, and H2Se are in general less stable than the analogous S species (Table 2), we consider 6 % to be a maximum value for the mass fraction of anthropogenic Se emissions that come from each of these species. The anthropogenic speciation fractions for OCSe, CSe2, CSSe, and H2Se are varied between 0 % and 6 %. The rest of the anthropogenic emissions, 76 %–100 %, are attributed to OX_Se_IN, representing species such as Se and SeO2.

Regarding the speciation of volcanic S emissions, Watts (2000) estimates that 0.99±0.88Tg S yr−1 is in the form of H2S. Comparing this to the estimate for volcanic SO2 emissions, 12.6 Tg S yr−1 (Andres and Kasgnoc1998; Dentener et al.2006), H2S contributes at most 13 % of volcanic S emissions. Therefore, in our sensitivity analysis the percentage of volcanic Se emissions in the form of H2Se ranges from 0 % to 13 %. Conversely, the percentage of OX_Se_IN in volcanic emissions ranges from 87 % to 100 %.

3.1.6 Accommodation coefficient of oxidized Se

The accommodation coefficient represents the probability that a gas-phase oxidized Se molecule will stick to an aerosol particle upon collision. No laboratory studies have investigated the accommodation coefficient of oxidized Se on aerosol surfaces. However, review papers suggest that Se efficiently partitions to the aerosol phase upon oxidation (Mosher and Duce1987; Wen and Carignan2007), indicating a high accommodation coefficient. Due to the lack of further information, we assume an uncertainty range of 0.02–1 for the accommodation coefficient, as selected by Lee et al. (2011) for H2SO4. This accommodation coefficient applies to uptake of Se on sulfate and dummy aerosols.

3.1.7 Dummy aerosol tracers

SOCOL-AER only includes sulfate aerosol, lacking other aerosol types (e.g., dust, sea salt, organic aerosol) that may also transport Se. To test how these other aerosol types might affect Se transport and deposition, in the sensitivity analysis we vary the emission location, particle radius, and emission magnitude of the dummy aerosol tracer (introduced in Sect. 2.2.4). These input parameters are different from other uncertainties, in that they are intended to investigate a lack of model completeness rather than uncertainties in measurable quantities. In our experiments the aerosols are emitted in one of 18 latitude bands, ranging from 90–80 S to 80–90 N. The latitude parameter can demonstrate whether Se deposition is affected by these missing aerosol types only in certain latitudinal bands. The emission radius affects the transport of Se on these missing particles, since the atmospheric lifetime of these particles depends on radius (Seinfeld and Pandis2016). For example, additional coarse particles (radius >1µm) could inhibit far-range Se transport, since they sediment and are removed quicker than the accumulation mode sulfate particles. We vary the emission radius between 0.01 and 3 µm, as this covers the range of atmospherically relevant particle sizes.

To determine a reasonable range for the emission magnitude of additional aerosol particles, we analyzed the particle emission inventories from the AEROCOM I project (Dentener et al.2006). Aerosol types are segregated into size classes based on their effective radius, and the emission magnitude for 10 latitude bands is calculated for each size class. The mass emission of particles correlates with the particle size, with generally larger mass emissions for larger particle sizes (Fig. S1 in the Supplement). Therefore, the value of emission magnitude in our experiments depends lognormally on the particle radius (r), with mean μ=2.54r+19.0 and standard deviation σ=2.88. We acknowledge that in the real atmosphere particles are emitted globally as size distributions, not monodisperse particles in one latitudinal band. Nevertheless, by including these input parameters as part of the sensitivity analysis, we can identify whether the lack of tropospheric aerosols other than sulfate impacts the simulated Se deposition maps in SOCOL-AER.

3.2 Experimental setup and model outputs

The creation of surrogate models requires a set of training runs with the full SOCOL-AER model. The values of the input parameters are varied simultaneously between training runs, so that interactions between parameters can also be detected. A Latin hypercube design is used to draw N samples from an M-dimensional input parameter space. In Latin hypercube sampling, each parameter range is divided into N equally probable subintervals. N values for each parameter are drawn randomly from within each subinterval. The selected values for all input parameters are then matched randomly with each other, to yield N points that cover the parameter space better than purely random Monte Carlo sampling (McKay et al.1979). The general rule of thumb is to select around 10M training simulations, to adequately cover the sample space (Loeppky et al.2009). We ran N=400 training simulations in our sensitivity analysis with 34 input parameters. Producing training simulations from the full SOCOL-AER model is the most computationally expensive step in the uncertainty and sensitivity analysis. With 48 cores, each training run takes around 14 h, corresponding to around 109 core seconds for 400 training simulations.

The initial conditions file for the training runs was created from a previous 10-year spinup of the model under year 2000 conditions. The atmospheric mixing ratios of Se tracers are initially set to 0 in each simulation. The simulations are 18 months long, with the first 6 months considered to be an equilibration period for the Se tracers. Therefore, we analyze only the last 12 months of the 18-month simulation. The model is run with T42 horizontal resolution (approximately 2.8×2.8 or 300 km×300 km in latitude and longitude) with 39 vertical levels up to 0.01 hPa (∼80km). The interactions between chemical species (i.e., greenhouse gases and aerosols) and radiation are decoupled in our simulations. The decoupling between chemistry and climate prevents meteorological differences between training runs with different input parameters, eliminating the influence of precipitation variability on deposition. Deposition flux differences in these relatively short simulations can then be more easily attributed to changes in the input parameters.

All relevant Se cycle fluxes and reservoir burdens are output by SOCOL-AER as monthly averages. For the sensitivity analysis, the target outputs are annual mean values of the global atmospheric Se lifetime and Se surface deposition fluxes. The global Se lifetime is calculated as the total atmospheric Se burden divided by total deposition (Seinfeld and Pandis2016). Deposition fluxes of Se were calculated by summing up the dry and wet deposition fluxes of aerosol and gas-phase Se species. Since the model is run in T42 horizontal resolution, there are 8192 surface grid boxes representing geographical coordinates on the globe. We therefore have 8192 model outputs for Se deposition and we construct a PCE model and conduct a sensitivity analysis for each one. It can be argued that the computational cost can further be reduced by dimensionality reduction techniques (e.g., Blatman and Sudret2011b, 2013; Ryan et al.2018), like building the PCE models on a reduced output set coming from a principal component analysis (PCA). We did not consider such an approach in this work because the cost of building the 8192 PCE models is marginal compared to that of evaluating the full SOCOL-AER model. In summary, the 400 SOCOL-AER training runs yield 400 values for atmospheric Se lifetimes and 400×8192 values for deposition fluxes.

3.3 Surrogate modeling with PCE

Sudret (2008) provided a detailed description of using PCE to build surrogate models, which was updated by Le Gratiet et al. (2017). We will summarize the most important features of the approach here. For this study, a certain output of the SOCOL-AER model (Y) can be thought of as a finite variance random variable which is a function of the M=34 input parameters (X=[X1,X2,,X34]):

(3) Y = M ( X ) .

The input X∈ℝM is modeled by a joint probability density function (PDF) fX whose marginals are assumed independent, i.e., fX=i=1MfXi(xi).

In polynomial chaos decomposition, the output variable Y is decomposed into the following infinite series (Ghanem and Spanos2003):

(4) Y = α = N M y α Ψ α ( X ) ,

where yα are coefficients to be determined and α is a multi-index set that defines the degree of the multivariate orthonormal polynomial Ψα(X)=i=1MΨαi(Xi). The latter belongs to a family of polynomials that are orthogonal with respect to the marginal PDF fXi. For example, univariate uniform probability distributions correspond to the family of Legendre polynomials, and Gaussian probability distributions correspond to the family of Hermite polynomials (Xiu and Karniadakis2002). Multivariate orthogonal polynomials can be constructed through multiplying the relevant univariate polynomials. The order of a polynomial term is defined as the total number of variables included in the polynomial term. The degree of a polynomial term represents the sum of the exponents of all variables appearing in the term. The number of coefficients to estimate grows exponentially with both the dimension and the degree. To allow calculation of the coefficients, the terms in Eq. (4) are truncated by restricting the maximum degree of the polynomials. Advanced truncation schemes can also be used to reduce the number of terms and thus the computational budget. Here we consider hyperbolic truncation which removes high-order interaction terms from the PCE, while retaining high-degree polynomials of a single variable (Blatman and Sudret2011a). Hyperbolic truncation is based on selecting only the terms that satisfy

(5) A = α : i = 1 M α i q 1 / q p ,

where q is a value selected between 0 and 1, and p is a value selected as the maximum degree of the PCE. Setting q=1 corresponds to the standard truncation scheme where all terms with a degree below p are selected. The lower the value of q, the more higher-order interaction terms are removed from the PCE. Based on previous experience, we selected a q value of 0.75.

PCE coefficients are generally calculated by least-square regression (Berveiller et al.2006) using an experimental design which consists of uniformly sampled realizations of the input variables 𝒳={x(1)x(N)} and corresponding model evaluations {ℳ(x(1))…ℳ(x(N))}, i.e.,

(6) y α = arg min y α R card A 1 N i = 1 N M x ( i ) - α A y α Ψ α x ( i ) 2 .

When the dimension is large, regression techniques that allow for sparsity, i.e., by forcing some coefficients to be zero, are favored. In this work, we consider least-angle regression as proposed in Blatman and Sudret (2011a) and follow the implementation in the UQLab PCE module (Marelli and Sudret2019).

The accuracy of the PCE in representing the full SOCOL-AER model is evaluated with a cross-validation metric named the leave-one-out (LOO) error, ϵLOO (Blatman and Sudret2010a). The leave-one-out procedure consists of creating a PCE using all but one of the training simulations, PC∕i. This PCE is then used to predict the output value of the model at the removed training simulation point, x(i). The process is repeated for all N training points so that a residual sum of squares can be calculated. This residual is then divided by the output variance in the training dataset, yielding the LOO error:

(7) ϵ LOO = i = 1 N ( M ( x ( i ) ) - M PC / i ( x ( i ) ) ) 2 i = 1 N ( M ( x ( i ) ) - μ ^ Y ) 2 ,

where μ^Y is the sample mean of the model output in the training dataset. In practice, the LOO error is estimated using a single PCE model considering the entire experimental design (Marelli and Sudret2014; Blatman and Sudret2010b), to avoid the procedure of creating N PCE models. Equation (7) is evaluated as

(8) ϵ LOO = i = 1 N M ( x ( i ) ) - M PC ( x ( i ) ) 1 - h i 2 / i = 1 N M ( x ( i ) ) - μ ^ Y 2 ,

where hi is the ith component of the vector defined by

(9) h = diag A ( A T A ) - 1 A T ,

and A is the experimental matrix whose components read

(10) A i j = Ψ j x ( i ) i = 1 , , N ; j = 0 , , P - 1 ,

where P≡card𝒜 is the number of terms in the PCE. In our study, we use a degree-adaptive scheme to construct the PCE models sequentially from degree 1 to a maximum of degree 13. If the LOO error does not decrease over two steps in degree, the algorithm is stopped and the maximum degree is selected as the one with the lowest LOO error. This method reduces the risk of overfitting the training set. The maximum degree of 13 was selected because the PCE becomes too computationally expensive to calculate above this degree. In any case, almost all PCEs calculated in this study are below degree 10.

To improve the accuracy of the approximation, we applied post-processing steps to the construction of PCE models (Table 4). The global atmospheric Se lifetime is a function of the atmospheric Se burden and the total Se deposition. We found improved overall accuracy when separate PCE models were created for the atmospheric burden and total deposition, rather than first calculating the lifetime for each simulation and constructing a PCE model of the lifetime. To calculate a surrogate model of the Se lifetime, we divided the PCE model of the atmospheric Se burden by the PCE model of total deposition. When calculating the PCE models of Se deposition fluxes, we first used the deposition fluxes directly from the training simulations. However, 873 of the 8192 grid boxes showed LOO errors higher than 0.05, which was selected as the acceptable threshold for this study. Improved LOO errors were achieved when the simulated deposition fluxes were log-transformed before constructing the PCE model. However, after log-transformation the sensitivity indices for deposition cannot be extracted analytically from the PCE (Sect. 3.4), increasing computational expense. Therefore, we only log-transformed the data from these 873 grid boxes before creating their PCE models. Surrogate models for deposition fluxes in these 873 grid boxes are calculated by taking the exponential of the PCE model of log-transformed data.

The cross-validation approach would usually remove the need for an independent validation dataset, saving computational expense. However, to evaluate the post-processing steps applied to the surrogate models, we also produced an independent validation dataset of 50 SOCOL-AER runs. The parameters for these runs were chosen by enriching the training experimental design to create a pseudo-Latin hypercube of 450 runs, ensuring that the distance between the validation runs and existing training runs is maximal.

Table 4Summary of methods to construct surrogate models and calculate Sobol' sensitivity indices of the SOCOL-AER output parameters.

Download Print Version | Download XLSX

3.4 Sensitivity analysis

A Sobol' sensitivity index represents the fraction of model variance caused by the parametric uncertainty of a certain input variable or the interaction between multiple variables (Sobol1993; Marelli et al.2019). It is a global measure, meaning that the sensitivity index considers the entire parametric uncertainty space. Let us consider a model whose inputs, defined over a domain 𝒟X, are assumed independent. It can be shown that it admits the following so-called Sobol' decomposition:

(11) M ( x ) = M 0 + i = 0 M M i ( x i ) + 1 i j M M M i j ( x i , x j ) + + M 1 , 2 , , M ( x 1 , , x M ) ,

where 0 is a constant and the other summands are defined such that their integrals with respect to any of their arguments are 0, i.e.,

(12) D X M i 1 , , i s ( x i 1 , , x i s ) f x i k d x i k = 0 , 1 k s .

Following this decomposition, it can be shown that the variance of the random variable Y=ℳ(X) can be cast as (Saltelli et al.2008)

(13) D = Var [ M ( X ) ] = i = 1 M D i + 1 i j M M D i j + + D 12 M ,

where a partial variance with respect to several variables Xi1,,Xis can be calculated as

(14) D i 1 , , i s = D X M i 1 , , i s 2 x i 1 , , x i s f X i 1 x i 1 f X i s x i s d x i 1 d x i s .

The Sobol' indices (Si1,,is) for a subset of input parameters are eventually obtained by normalizing the corresponding variance, i.e.,

(15) S i 1 , , i s = D i 1 , , i s D .

The number of variables involved in a Sobol' index determines its order. A first-order Sobol' index (Si=DiD) apportions the variance in the model output to the effect of a single variable, Xi. Second-order indices (Sij=DijD) represent the impact of the interaction of two variables (e.g., Xi and Xj, ij) on the model output variance, not already accounted for by Si and Sj. Higher-order indices can be calculated as well. The summation of all individual Sobol' indices yields a value of 1, i.e., accounting for all of the variance in the output. However, the number of higher-order indices to calculate can become very large:

(16) number of n th - order indices = n choose M = n ! M ! ( n - M ) ! .

It therefore becomes computationally demanding in our case with M=34 input parameters to calculate sensitivity indices with an order higher than 2. Furthermore, by the sparsity-of-effect principle (Goupy and Creighton2006), it is expected that the model is primarily driven by first-order effects and low-order interactions. We therefore only individually calculate the first- and second-order Sobol' indices, which is common practice in global sensitivity analysis.

The total Sobol' index (SiT) summarizes all sensitivity indices which include the effect of a given input variable:

(17) S i T = S i + j i S i j + j i k i k j S i j k +

In practice, it is possible to calculate the total Sobol' index without individually calculating all of the higher-order indices (Marelli et al.2019). The total Sobol' indices can be used to rank the influence of input variables on the variability in the model output. It should be noted that the sum of all total Sobol' indices will be greater than 1 if the model is nonadditive, since interaction terms would be counted multiple times.

Other studies have emphasized the computational expense of conducting a sensitivity analysis for all grid boxes in a chemistry–climate model (Ryan et al.2018). By using PCE as the surrogate model, we can reduce the computational expense since it is possible to compute the Sobol' sensitivity indices analytically. As shown by Sudret (2008), Sobol' sensitivity indices are a function of the calculated coefficients in the PCE model (Eq. 4), due to the similarity of the PCE decomposition with variance decomposition. First-order Sobol' indices can be calculated as Di=αAiyα2, where 𝒜i is the set of polynomial terms involving only variable i. Similarly, higher-order Sobol' indices can be calculated by summing the squares of the coefficients of the polynomial terms that include the variables of interest. The analytical solution of Sobol' indices is no longer possible when post-processing steps are applied to the original PCE models to create surrogate models. For example, (1) the surrogate model for the atmospheric lifetime is calculated from two PCE models, one for total Se burden and one for total Se deposition; (2) Se deposition in 873 grid boxes is calculated using the exponential of a PCE model (see Table 4). In these cases, it is still possible to calculate the Sobol' indices through Monte Carlo estimation (Marelli et al.2019). To accomplish this, the surrogate models are sampled around 106 times, which remains tractable.

To summarize categories of input variables, we aggregate the Sobol' indices of several input parameters to yield a total Sobol' index for that category. For example, we summarize the total dummy aerosol effect by summing the total Sobol' indices of the dummy aerosol radius, emission magnitude, and latitude. It would anyways be difficult to separate the effects of the dummy aerosol input parameters since they are correlated inputs in the experimental design. The second-order indices involving two dummy aerosol input parameters may be double-counted with this method. However, since these indices are small (<0.05) we do not expect a large error in the aggregated total Sobol' index (Sect. 5.1).

3.5 Resampling of surrogate models

In order to estimate distribution statistics (mean, standard deviation, quantiles) of the output variable, we resample each surrogate model 40 000 times. We also use these 40 000 samples of the parameter space to calculate relationships between input parameters and output variables. To visualize marginal relationships between a certain input parameter and the output, we replace the value of the input parameter in the 40 000 samples by a fixed value and calculate the mean and variance of the surrogate model output. Repeating this step with other evenly spaced values in the input parameter range, we can produce the univariate relationship between the model output and the input parameter.

4 Compilation of Se wet deposition flux data

We decided to compare SOCOL-AER results with measurements of Se wet deposition, since wet deposition contributes an estimated 80 % of total deposition (Wen and Carignan2007). We conducted a systematic literature review to assemble a dataset of measured wet Se deposition fluxes, extending an earlier review from Conde and Sanz Alaejos (1997). We searched in Web of Science (Clarivate Analytics) for combinations of the following criteria: “Selenium”, “Se”, “rain”, “precipitation”, “wet deposition”, and “trace element”. The last search was completed in July 2019, yielding a total of 672 search results. We screened these search results for studies that measured Se concentrations in rainwater for at least 1 month, neglecting studies that measured bulk (wet and dry) deposition or that extrapolated wet deposition fluxes from aerosol measurements. The compiled dataset, which is available in the Supplement, consists of 29 papers and data from two measurement networks, the European Monitoring and Evaluation Programme (EMEP) and Canadian National Atmospheric Chemistry (NAtChem) database, for a total of 73 measurement sites (Table 5). From these studies, we extracted the annual mean Se deposition fluxes, if available, or we used rainwater volume-mean weighted Se concentrations combined with the annual precipitation depths to calculate the deposition flux. If the paper did not provide the annual precipitation depth, we calculated the mean annual precipitation depth for the time period and location of the study from the historical Multi-Source Weighted-Ensemble Precipitation (MSWEP) dataset (Beck et al.2019). The calculated annual deposition fluxes for nine sites are extrapolated from fewer than 12 months of rainwater Se measurements; the majority of sites (64) were measured for longer than a year. For studies that spanned multiple years, a multi-annual mean deposition flux was calculated to compare with the model. Additional metadata were extracted from the papers, including geographic coordinates, the time period, collection methods, sampling frequency, and analytical methods. Despite the fact that the model is based on year 2000 conditions for emission maps and meteorology, we compare the model with data from all measurement years, since the dataset is relatively small. We created surrogate models of the Se wet deposition flux in the grid boxes where measurements were made (as in Sect. 3.3 for the total deposition fluxes). We can then compare the resampled distribution of Se wet deposition in these model grid boxes with the corresponding observation.

Suzuki et al. (1981)Arimoto et al. (1985)Cutter and Church (1986)Arimoto et al. (1987)Dasch and Wolff (1989)Heaton et al. (1990)Cutter (1993)Scudlark et al. (1994)Al-Momani et al. (1997)Cutter and Cutter (1998)Lawson and Mason (2001)De Gregori et al. (2002)Scudlark et al. (2005)Sakata et al. (2006)Shimamura et al. (2007)Landing et al. (2010)Zhou et al. (2012)Liu et al. (2012)Gratz et al. (2013)Pan and Wang (2015)Lynam et al. (2015)Xing et al. (2017)Nie et al. (2017)Blazina et al. (2017)Savage et al. (2017)Uchiyama et al. (2019)Suess et al. (2019)Pearson et al. (2019)Savage et al. (2019)

Table 5Previous studies measuring Se wet deposition fluxes.

Download Print Version | Download XLSX

5 Results

5.1 Atmospheric Se lifetime

From the 400 training runs of SOCOL-AER, we created PCE models of the global and annual mean total Se burden and deposition flux. The LOO error of the global burden is around 0.02 and the LOO error of the deposition flux is on the order of 10−6. The LOO error is low for total deposition since for a mass-conserving model total deposition should equal the sum of the emission input parameters. Indeed, all 400 training runs showed very good Se mass conservation: total annual Se deposition fluxes were 98 %–102 % of total emission fluxes. We derive a surrogate model for the atmospheric Se lifetime by dividing the Se burden PCE model by the Se deposition PCE model (Sect. 3.3). This surrogate model shows a higher error (0.35) than the burden and deposition flux PCE models, which would only be reduced by running more training runs (Fig. S2). The surrogate model for the atmospheric lifetime is resampled to calculate the probability distribution of this output, given our assumptions for the uncertainty ranges of all 34 parameters (Sect. 3.5). The histogram of the atmospheric lifetime (Fig. 4) shows a narrow range for the atmospheric Se lifetime, with a median of 4.4 d (days), 2nd percentile value of 2.9 d, and 98th percentile value of 6.4 d. Despite the large uncertainty ranges for the reaction rate constants of Se, spanning multiple orders of magnitude (Table 3), the uncertainty of the simulated atmospheric lifetime is less than 1 order of magnitude.

Figure 4Distribution of the atmospheric Se lifetime, resampled from the surrogate model. Summary statistics (median, 2nd percentile, and 98th percentile values) are listed on the plot.


In order to identify the input parameters that drive the variability in the simulated Se lifetime, we apportion the variance into contributions from the most important parameters (Fig. 5). The most important variable is k13, which is the rate constant for the OCSe+OH reaction, followed by the dummy aerosol input parameters. Nonlinearities are also clearly important for the global Se lifetime, since all first-order terms only account for 53 % of the variance in the Se lifetime, with interaction terms accounting for the other 47 %. However, the interaction contribution is made up of many small individual interaction terms. Only two interaction terms account for more than 5 % of the total variance in the Se lifetime: the interaction between k1 and k13 (5.3 % of variance) and the interaction between the dummy aerosol radius and dummy aerosol emissions (5.0 % of variance). Through resampling the surrogate model for the Se lifetime, we investigate the relationships between the Se lifetime and input parameters in Figs. 6 and 7.

Figure 5Sensitivity indices of the most important parameters (SiT>0.05) for the atmospheric Se lifetime. Total Sobol' indices are partitioned into first-order and interaction effects. The total Sobol' indices for the dummy aerosol parameters are aggregated, since it is difficult to separate the effects of the correlated input parameters (aerosol radius and emission magnitude).


Several of the most influential input parameters for the Se lifetime are related to OCSe processes (Fig. 6a–c). Given the median estimates for the reaction rate constants in Table 2, OCSe has the longest lifetime of any gas-phase Se species. Therefore, chemical reaction rates and emissions of OCSe have a strong effect on the overall Se burden and Se lifetime. Slower reaction rates of OCSe with OH (k13) and O (k1) lead to longer Se lifetimes. Since OH is more prevalent in the lower atmosphere than O, the dependence of the Se lifetime on k13 is stronger than the dependence of the lifetime on k1. The influence of k13 on the Se lifetime is mostly saturated above 10−12cm3molec.-1s-1. Above this threshold in k13, the OCSe burden becomes a minor part of the total Se burden (accounting for less than 2 % of total Se), and therefore the Se lifetime is not affected by OCSe processes. The second-order interaction between k13 and k1 is also shown in Fig. 6b. When the reaction of OCSe with OH is fast (k13>10-12.5cm3molec.-1s-1), the value of k1 is not important for the Se lifetime since almost all OCSe will react with OH. In cases when the OCSe + OH reaction is slow (k13<10-12.5cm3molec.-1s-1), the value of k1 has a stronger effect on the Se lifetime since not all OCSe has reacted with OH. The Se lifetime increases for higher fractions of anthropogenic Se emitted as OCSe, again showing that higher OCSe burdens prolong the atmospheric Se lifetime. OCSe has been associated with anthropogenic emissions by only one study, which inferred its presence from a similarity in boiling point with a detected Se species (Pavageau et al.2002). OCSe has never been identified in the ambient atmosphere; on the other hand, SOCOL-AER predicts maximum OCSe concentrations of sub-parts-per-trillion levels, which would be difficult to measure analytically. Still, processes related to OCSe, a highly speculative atmospheric Se compound, contribute to the model's variability in the Se lifetime.

Figure 6Relationships between the atmospheric Se lifetime and the most important input parameters from Fig. 5. Using resampling techniques (Sect. 3.5), we calculate the mean and standard deviation Se lifetime over the range of input parameters. Interaction effects are illustrated in the k1 plot (b), by grouping the samples into cases where k13 is high and low.


After OCSe, the most important input parameters for the atmospheric Se lifetime are related to dummy aerosols. It is difficult to interpret the univariate effects of aerosol emissions and radius (Fig. 7a, b), since by design the two inputs are correlated. We hypothesize that the influence of dummy aerosols on the Se lifetime is related to the uptake rate on these dummy aerosols, which would be proportional to the surface area density. The mass emissions of the aerosols are proportional to volume (r3), which when divided by dummy aerosol radius yields a metric that is proportional to the emitted surface area (r2). The Se lifetime increases monotonically with the surface area of emitted dummy aerosols, although the response is not linear (Fig. 7c). We compare this response with the range of available surface area of sulfate aerosol for 10 latitude bands, shaded in yellow. For lower dummy aerosol surface areas, the available sulfate aerosol dominates the uptake of gas-phase Se, and additional dummy aerosols do not play an important role. The dummy aerosols have a stronger effect on Se lifetime at the upper limit of the sulfate aerosol surface area range. The absence of aerosols other than sulfate leads to a lower Se lifetime in SOCOL-AER. However, the effect of dummy aerosols is not drastic, only increasing the mean Se lifetime from 4.3 to 4.8 d (Fig. 7c).

Figure 7Relationship between the atmospheric Se lifetime and the dummy aerosol input parameters (a, b). A quantity related to surface area is shown in (c), which also includes the comparable range for sulfate aerosol as yellow shading. For sulfate, this is calculated by dividing the deposited mass of sulfate (equal to emitted mass in steady state) in each 10 latitudinal band by the effective sulfate radius in that latitudinal band. Since aerosol mass emissions are a lognormally distributed quantity, the bounds are infinite and thus we include only the 5th to 95th percentile range in (a, c).


The other inputs have minor impacts on the global Se lifetime. Stronger emissions of volcanic Se lead to shorter overall global Se lifetimes, while emissions of marine and terrestrial biogenic Se lead to longer Se lifetimes (Fig. 6d–f). In SOCOL-AER, biogenic sources emit DMSe, which is not removed by wet and dry deposition, while volcanic sources emit mainly oxidized Se species, which can be removed by wet and dry deposition. Biogenic emissions must first be oxidized before deposition can occur, which can prolong the Se lifetime. The influence of the two other terms with SiT>0.05, the reaction rate constant of CSSe+OH (k15) and the fraction of volcanic emissions emitted as H2Se, is mainly through interaction terms, and therefore the Se lifetime responds only weakly to univariate variations in these variables (Fig. 6g, h).

5.2 Spatial patterns of Se deposition

Surrogate models for Se deposition in all surface grid boxes were calculated according to Sect. 3.3. After log-transforming Se deposition in 873 grid boxes, the LOO error is below 0.1 in all grid boxes, but still above 0.05 in 354 grid boxes, mainly in polar regions (Fig. S2). However, due to limitations in computational time to run more training runs, we proceeded with the sensitivity analyses of all grid boxes. By resampling the surrogate models, we calculate maps of mean Se deposition and its coefficient of variation in all input uncertainties (Fig. 8). Selenium deposition is highest in areas close to anthropogenic and volcanic emissions (Fig. 3c, d), which are point sources. The deposition pattern is clearly affected by precipitation: dry areas (e.g., eastern portion of ocean basins, polar regions, and the Sahara) show the lowest Se deposition fluxes globally. Other interesting features of the Se deposition patterns, for example identifying the regions influenced by long-range transport, will be investigated in upcoming studies. In this study, we focus on determining the most important sources of uncertainty to the simulated deposition maps. The relative uncertainty in simulated deposition, illustrated by the coefficient of variation (Fig. 8b), is highest in areas affected by marine and volcanic emissions because these emission sources have wider uncertainty ranges compared to anthropogenic and terrestrial emissions (Table 3). In the following global sensitivity analysis we can identify the input factors that contribute to the variation in Se deposition in each grid box.

Figure 8Map of the mean Se deposition (a) and associated coefficient of variation (b), calculated by resampling surrogate models of deposition in each grid box.

Figure 9Maps of the total Sobol' indices of emission parameters and the accommodation coefficient for total Se deposition (left column). A blue circle indicates the grid box where the total Sobol' index is at a maximum. The relationship between total Se deposition and the input parameters in that grid box, calculated by resampling the surrogate model for deposition, are shown (right column). Note that the magnitude of deposition (y axis) varies in each plot, depending on the grid box shown.

Figures 9 and 10 illustrate the spatial variation in the importance of input parameters. We chose example grid boxes (indicated by blue circles) to illustrate the marginal relationships between Se deposition and the input parameters. Overall, the most important input parameters are the total emissions from each source. In the example grid boxes shown in Fig. 9a–h, Se deposition increases linearly with increasing emissions from the different source categories. The linear relationship is logical since deposition balances emission in the steady state. In areas that are more remote from emission regions (the Sahara, Antarctic, and Arctic), other factors become more important but are still minor compared to the emission inputs. The aerosol accommodation coefficient affects areas where the precipitation is very low, for example in the Sahara (Fig. 9i). In these dry regions, total Se deposition is dominated by dry deposition. When the model is run with low accommodation coefficients, less oxidized Se partitions to the particulate phase and more remains in the gas phase. Dry deposition of particles in the 0.1–1 µm diameter size range, within the range of sulfate and dummy aerosols, is slower than gas compounds due to the slower Brownian motion of particles (Seinfeld and Pandis2016). Chemical reaction rate constants, specifically the reaction rates of DMSe, impact Se deposition in polar regions. Slower reaction rates of DMSe with OH (k6) and O3 (k8) enhance deposition over the example Antarctic grid box (Fig. 10a–c). Longer DMSe lifetimes allow more marine Se emissions to reach polar regions, which have few local Se emission sources. The chemical rate coefficients are more important in the Antarctic than the Arctic since there is more O3 and OH in the Northern Hemisphere than the Southern Hemisphere, meaning that the DMSe lifetime is longer in the Southern Hemisphere than the Northern Hemisphere. Dummy aerosol parameters are only important for Se deposition in the Arctic (Fig. 10d–f). We calculate the relationship between Se deposition and dummy aerosols using the same quantity for emitted surface area of the dummy aerosols as in Fig. 7c. With increasing dummy aerosol surface area, the Se deposition in the Arctic increases, but only after surpassing the threshold of the available sulfate aerosol surface area in that latitude band (Fig. 10e). The dummy aerosols also have a stronger effect on Se deposition when they are emitted in a latitude band closer to the example grid box at 76.7 N (Fig. 10f). Attachment of oxidized Se to dummy aerosols increases the overall lifetime of Se (Sect. 5.1), leading to enhanced transport of Se to the Arctic region. The transport of Se on dummy aerosols does not lead to higher deposition in the Antarctic, perhaps because wet deposition in the Antarctic circumpolar storm track impedes the transport of aerosol poleward.

It is also important to note which input parameters do not influence Se deposition in any of the grid boxes. Variations in the speciation of emissions, photolysis rates, and 15 of the Se reaction rate constants have a negligible influence on deposition.

Although other parameters may play a role in certain grid boxes, the emission parameters are most important on the global scale, evidenced by their higher mean total Sobol' index (Fig. 11). Figure 9 illustrates which regions are affected by different emission sources. Variations in marine emissions impact the most grid boxes; however, their influence is mainly confined to the oceans, coastal areas, and Southern Hemispheric continents. Since the motivation of studying Se deposition is to understand its impact on agricultural soils, we also calculated the mean influence of parameters in pasture and cropland areas (Fig. 11), using maps from Ramankutty et al. (2008). The importance of anthropogenic emissions increases when looking only at pasture and cropland areas, since agricultural areas coincide with human settlement. All four emission parameters show similar levels of importance for agricultural regions, i.e., showing similar mean Sobol' indices for pasture and cropland areas. Therefore, further work in understanding any of these emission processes would be valuable to reducing the uncertainty in deposition fluxes.

Figure 10Maps of the total Sobol' indices of reaction rate constants and dummy aerosol parameters for total Se deposition (a, d). A blue circle indicates the grid box where the aggregated total Sobol' index is at a maximum. The relationship between total Se deposition and the multiple relevant input parameters within the aggregated index for that grid box is shown (b, e). The emitted dummy aerosol surface area is compared to the corresponding sulfate quantity in the latitude band of the grid box, shown as a dashed brown line in (e).

Figure 11Bar plot summarizing the importance of the input parameters to total Se deposition globally and in agricultural areas. For the cropland and pasture means, Sobol' indices are averaged over grid boxes that are covered by more than 25 % cropland or pasture area in the Ramankutty et al. (2008) database.


5.3 Comparison with deposition flux measurements

With surrogate models of wet Se deposition, we can estimate the modeled wet Se deposition throughout the parametric uncertainty space. By comparing these modeled distributions of wet deposition with observed values, one could constrain input parameters to which deposition is sensitive. However, we do not attempt at this stage to calibrate the parameters to existing measurements because of several challenges in comparing the compiled measurement dataset with the simulations in this study. Firstly, the emissions and meteorology in this study are representative of the year 2000, whereas the measurements were made between 1975 and 2017. Secondly, Se is a difficult element to measure at environmental concentrations, which might lead to inaccurate reported deposition fluxes. The most popular analytical method is inductively coupled plasma mass spectrometry (ICP-MS), which was used for 58 of the 73 sites in the database. However, it is difficult to measure Se with ICP-MS due to the low ionization of Se, the Se signal being split on the five stable isotopes, and especially mass interferences (Winkel et al.2012). Several studies reported that Se concentrations in rainwater samples were often below the detection limit of the analytical method (e.g., Arimoto et al.1987; Gratz et al.2013). Unfortunately, other studies often do not explicitly report the detection limit, the fraction of samples under the detection limit, and how these samples are treated statistically. Thirdly, many of the measurement sites were located in urban locations close to point-source emissions. Due its coarse resolution (2.8×2.8), the model would have difficulty reproducing point values for Se concentrations and deposition fluxes.

Figure 12 compares the measured Se wet deposition fluxes with the resampled median deposition fluxes from the surrogate models, also showing the likely range predicted by the models (defined with bounds of the 2nd and 98th percentile values). With the results from the deposition sensitivity analysis (Fig. 9), we categorize each measurement location by the input parameter that predominates the uncertainty in modeled deposition, indicated by the symbol in Fig. 12. The compiled data and measurements show good agreement at the lower end of deposition values, where the measurement sites are more remote from point-source emissions. The agreement worsens at higher values of observed Se deposition, which correspond to more urban measurement sites. As discussed before, it is not surprising that the model has difficulty matching the measurements for higher observed values of deposition, since the model has coarse resolution and the simulation year may be mismatched from the measurement year. Indeed, the model underestimates several Se deposition measurements at urban sites from East Asia after 2005. Anthropogenic SO2 emissions, an analogue of anthropogenic Se emissions, have increased in East Asia since 2000 (Smith et al.2011). Natural emission factors dominate the variability in several of these locations in the year 2000 simulations, likely because the input anthropogenic emission maps do not correspond to the measurement time period. Nevertheless, we find overall that 53 % of existing measurements are within the likely range of the model's prediction. The agreement improves to 79 % when comparing the model range to measurements in background locations, defined as having observed deposition values below 150 µg Se m−2. These results provide confidence to SOCOL-AER's predictions of Se deposition fluxes in nonurban locations.

6 Discussion

Through our consideration of model uncertainties related to Se cycling, we derived a median atmospheric Se lifetime of 4.4 d. This is the first bottom-up estimate for Se made with a mechanistic global atmospheric model. Our estimate for the Se lifetime matches the global arsenic lifetime calculated in the global model of Wai et al. (2016), namely 4.5 d. This agreement is due to both elements attaching to submicrometer aerosol particles; therefore, their lifetime is determined by the lifetime of these particles. Since Se and arsenic have similar atmospheric sources as well (e.g., metal smelting, coal combustion, volcanoes), it may be possible to draw analogies between their atmospheric cycles. The range of previously estimated Se lifetimes from global atmospheric budgets is between 0.8 and 6 d, similar to our result (Ross1985; Mosher and Duce1987). The recent value from Mason et al. (2018) of a 0.15-year (55 d) Se lifetime seems overestimated compared to our results and past budgets, especially since Mason et al. (2018) only consider gas-phase Se in their model, which tends to be shorter lived in the atmosphere than aerosol-bound Se. According to our sensitivity analysis results, the atmospheric Se lifetime could be further constrained by measuring the OCSe+OH reaction rates, and in general knowing more about whether OCSe is present in the atmosphere. Since dummy aerosols also impact the Se lifetime in our model, implementing a more complex tropospheric aerosol parametrization in SOCOL-AER would also further constrain the atmospheric lifetime of Se. However, since the main interest in Se is its atmospheric input to agricultural soils, it may be a higher priority to constrain the input parameters that affect the deposition of Se in agricultural regions rather than the Se lifetime.

The results of the sensitivity analyses raise an obvious question: why do the input parameters that influence the atmospheric Se lifetime not appear as important for the Se deposition fluxes? One would expect that Se deposition fluxes close to areas of high emissions would be dominated by the magnitude of these emissions. One would also expect that variation in the Se lifetime would play a role, if anywhere, over remote regions, where the amount of locally emitted Se is low, and thus the amount that can be transported from emission regions has a larger effect on deposition. However, the range in the atmospheric Se lifetime in our simulations is relatively narrow, between 2.9 and 6.4 d if we consider the 2nd percentile and 98th percentile bounds (Fig. 4). On the other hand, emissions of various Se species can vary by orders of magnitude (Table 3). These larger variations in the amount of emitted Se have a larger impact on deposition than smaller variations in the Se lifetime, even in many remote places. Only in extremely remote areas, for example in the Arctic, do some of the parameters that affect the Se lifetime show up as important, like the dummy aerosols. Parameters with regional rather than global importance for the Se lifetime, like the DMSe reaction parameters, impact deposition of Se in the Antarctic by controlling the amount of transported Se. It is not surprising that the parameter that has the largest impact on lifetime, the OCSe+OH reaction rate constant, has little impact on deposition fluxes, since emissions of OCSe are assumed to be a minor flux of Se (maximum 6 % of the anthropogenic emissions flux). Like all sensitivity analyses, the results are dependent on the choice of uncertainty ranges for the different parameters; if we had selected narrower uncertainties for the Se emission sources, the uncertainties of parameters that affect Se lifetime (e.g., chemical reaction rates, dummy aerosols) may have been more important in remote regions. However, the choice of wide uncertainty ranges for the Se emissions is justified, given the variability in natural emission processes and a lack of field campaigns assessing Se emission fluxes (Sect. 3.1.4). The different results for the two types of sensitivity analyses (lifetime and deposition fluxes) highlight that the “important” parameters to constrain depend on the choice of research question.

It must be noted that our simulations were performed only for the year 2000 and focused on uncertainties in the Se-related input parameters, neglecting variations in the Se cycle due to interannual variability in meteorology and sulfate aerosol properties. In future studies we intend to investigate how the Se cycle varies under different climate conditions.

The global sensitivity analyses in this paper provide clear next steps for atmospheric Se research. The magnitude and spatial distribution of Se emissions remain the most important uncertainty to constrain, in order to improve the predictions of Se deposition patterns. Further investigations of chemical reactivity of Se species or the speciation of emissions are a lower priority, although measuring the speciation of emissions can give mechanistic insights into emission processes. The emission uncertainties could be constrained by conducting field campaigns that measure either emission fluxes of Se close to sources (e.g., Amouroux et al.2001) or separate Se source contributions at an ambient measurement site through trajectory modeling and/or speciation measurements (e.g., Suess et al.2019). Our model results can help identify ambient locations that would be interesting to study for field campaigns, by mapping the contribution of the Se emission sources to deposition in different regions (Fig. 9). In addition to new field measurements, we can also compile and reanalyze previously collected data from the literature to evaluate estimates of emission fluxes. Bayesian inverse modeling techniques (e.g., Stohl et al.2009) could be employed in conjunction with the SOCOL-AER model to provide posterior estimates for Se emission fluxes. Global sensitivity analysis is an invaluable first step before such model calibration techniques, since the parameter dimensionality can be reduced by neglecting uninfluential parameters. As shown in Sect. 5.3, the heterogeneity of compiled literature data represents a challenge to comparing models and measurements. Therefore, standardized measurement techniques and adequate reporting of sampling, analytical, and post-processing methods are required so that the model is not calibrated to an errant measurement.

Figure 12Comparison of wet deposition flux measurements (Table  5) with modeled fluxes. Medians for the modeled values are shown, along with vertical bars representing the 2nd and 98th percentile values. The symbols correspond to the model input parameter that is most important for deposition at the measurement location. The color of the symbols represents the year when the measurement was taken; multi-year measurements show the middle year.


The ultimate motivation for studying biogeochemical Se cycling is to better understand how this essential element enters food chains and ecosystems. It has been hypothesized that atmospheric Se could be an important source of Se to soils and thus terrestrial food chains (Winkel et al.2015); however, until now this could not be proven due to missing knowledge of global atmospheric Se cycling. Using SOCOL-AER, the first model of its kind for Se, we can quantify the atmospheric inputs of Se to terrestrial systems and investigate how these inputs are impacted by different climate conditions and anthropogenic activities. Although human health effects cannot be directly inferred from atmospheric deposition maps, SOCOL-AER will be coupled to soil, plant, and health system models in future work to accurately predict Se deficiency risks.

7 Conclusions

Now that it includes Se cycling, the SOCOL-AER model can be used to predict atmospheric Se transport and deposition globally. We created surrogate PCE-based models that are able to predict the output of the model throughout the uncertainty space of the input parameters. With these surrogate models, we determined that the atmospheric Se lifetime is around 4.4 d, similar to the lifetime of submicron aerosol particles in the atmosphere. Assuming that longitudinal wind speeds are around 10 m s−1 (Jacob1999), the likely Se lifetime range of 2.9–6.4 d corresponds to a distance of 2500–5000 km that Se is transported in the atmosphere. The global sensitivity analysis of Se deposition fluxes shows that reducing uncertainties in Se emissions would lead to the biggest reductions in the uncertainty of deposition maps. Field measurements that elucidate and quantify Se emission processes should be prioritized, so that model predictions of Se deposition maps can be improved. Available measurements of Se in rainwater are within the likely range of model results at 79 % of background sites; remaining discrepancies may be due to the time period of the simulations in this study, the coarse resolution of the model, and analytical challenges leading to measurement inaccuracies. In future work, deposition maps from SOCOL-AER can be linked to soil, plant, and health system models to identify the regions at risk of Se deficiency in current times and the future.

Code and data availability

The SOCOL-AER code is available upon request from the authors, after users have signed the ECHAM5 license agreement: (last access: 30 January 2020). The relevant simulation data, along with the experimental design of the training and validation runs, is available at (Feinberg et al.2019a). The compiled Se precipitation database is available in the Supplement. The NAtChem precipitation database is available online at (Environment and Climate Change Canada, 2020). Selenium in precipitation measured by EMEP was extracted from annual reports on heavy metals and POP measurements available at (EMEP/CCC, 2020) UQLab is freely available for academic and degree-granting institutions by registering at (last access: 30 January 2020). All of the scientific source code is available under the BSD three-clause license and can be downloaded from (last access: 30 January 2020).


The supplement related to this article is available online at:

Author contributions

AS, TP, and LHEW initiated the project of studying the atmospheric Se cycle using a chemistry–climate model. AF implemented Se into SOCOL-AER with assistance from AS, conducted all simulations, analyzed the results, and wrote the paper with contributions from all authors during the revision process. AF was directly supervised by AS, TP, and LHEW during the project. MM and BS guided the statistical analysis during the project. MM developed the scripts to conduct the global sensitivity analysis on the Se model and aided with writing the Statistical methods section.

Competing interests

The authors declare that they have no conflict of interest.


Thanks are due to Severin Gysin for his work on sensitivity analysis in a Se box model during his bachelor thesis at ETH. Thanks are due to David Amouroux and Emmanuel Tessier from CNRS/University of Pau for helpful discussions about marine DMSe emissions. Thanks are due to Franziska Genter for her help in compiling past studies for the Se precipitation database. We acknowledge Environment and Climate Change Canada's Monitoring and Surveillance in the Great Lakes Basin (GLB) under the Chemicals Management Plan (CMP) for the provision of Se in precipitation data and specifically Helena Dryfhout-Clark for providing additional information about this dataset. We acknowledge Barbara Trost and Anna Breuninger for providing data from the Alaska Department of Environmental Conservation. Thanks are due to Andrew Meharg for providing additional data of Se in precipitation. We acknowledge all scientists who contributed to measurements that were compiled into the Se precipitation datasets. Thanks are due to the developers of the UQLab for providing the software used in this study. Thanks are due to This Meschini for the design of the Se cycle schematic. Color tables from ColorBrewer 2.0 (, last access: 30 January 2020) were used for the figures in this paper.

Financial support

This research has been supported by the ETH Zurich (grant no. ETH-39 15-2).

Review statement

This paper was edited by Frank Dentener and reviewed by four anonymous referees.


Agrawalla, B. and Setser, D.: Infrared chemiluminescence and laser-induced fluorescence studies of energy disposal by reactions of F and Cl atoms with H2S (D2S), H2Se, H2O (D2O), and CH3OH, J. Phys. Chem., 90, 2450–2462,, 1986. a

Agrawalla, B. and Setser, D.: Vibrational distributions and rate constants from reactions of oxygen atoms with HI, GeH4, SiH4, H2Se, and H2S, J. Chem. Phys., 86, 5421–5432,, 1987. a

Al-Momani, I., Güllü, G., Ölmez, I., Eler, Ü., Örtel, E., Sirin, G., and Tuncel, G.: Chemical composition of eastern Mediterranean aerosol and precipitation: Indications of long-range transport, Pure Appl. Chem., 69, 41–46,, 1997. a

Amouroux, D. and Donard, O. F.: Maritime emission of selenium to the atmosphere in eastern Mediterranean seas, Geophys. Res. Lett., 23, 1777–1780,, 1996. a

Amouroux, D. and Donard, O. F.: Evasion of selenium to the atmosphere via biomethylation processes in the Gironde estuary, France, Mar. Chem., 58, 173–188,, 1997. a

Amouroux, D., Liss, P. S., Tessier, E., Hamren-Larsson, M., and Donard, O. F.: Role of oceans as biogenic sources of selenium, Earth Planet. Sc. Lett., 189, 277–283,, 2001. a, b, c

Andres, R. and Kasgnoc, A.: A time‐averaged inventory of subaerial volcanic sulfur emissions, J. Geophys. Res.-Atmos., 103, 25251–25261,, 1998. a, b, c, d

Arimoto, R., Duce, R., Ray, B., and Unni, C.: Atmospheric trace elements at Enewetak Atoll: 2. Transport to the ocean by wet and dry deposition, J. Geophys. Res.-Atmos., 90, 2391–2408,, 1985. a

Arimoto, R., Duce, R. A., Ray, B. J., Hewitt, A. D., and Williams, J.: Trace elements in the atmosphere of American Samoa: Concentrations and deposition to the tropical South Pacific, J. Geophys. Res.-Atmos., 92, 8465–8479,, 1987. a, b

Ariya, P. A., Amyot, M., Dastoor, A., Deeds, D., Feinberg, A., Kos, G., Poulain, A., Ryjkov, A., Semeniuk, K., and Subir, M.: Mercury physicochemical and biogeochemical transformation in the atmosphere and at atmospheric interfaces: A review and future directions, Chem. Rev., 115, 3760–3802,, 2015. a

Atkinson, R., Perry, R., and Pitts Jr., J.: Rate constants for the reaction of OH radicals with COS, CS2 and CH3SCH3 over the temperature range 299–430 K, Chem. Phys. Lett., 54, 14–18,, 1978. a

Atkinson, R., Aschmann, S. M., Hasegawa, D., Thompson-Eagle, E. T., and Frankenberger Jr., W. T.: Kinetics of the atmospherically important reactions of dimethyl selenide, Environ. Sci. Technol., 24, 1326–1332,, 1990. a, b, c, d

Beck, H. E., Wood, E. F., Pan, M., Fisher, C. K., Miralles, D. G., van Dijk, A. I., McVicar, T. R., and Adler, R. F.: MSWEP V2 global 3-hourly 0.1 precipitation: methodology and quantitative assessment, B. Am. Meteorol. Soc., 100, 473–500,, 2019. a

Belyaev, A., Sozin, A., Aleksandrov, Y., and Churbanov, M.: Homogeneous oxidation of hydrogen selenide with ozone in an argon-oxygen medium, Russian J. Appl. Chem., 85, 1919–1922,, 2012. a

Berveiller, M., Sudret, B., and Lemaire, M.: Stochastic finite elements: a non intrusive approach by regression, Eur. J. Comput. Mech., 15, 81–92,, 2006. a

Blatman, G. and Sudret, B.: An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probabilist. Eng. Mech., 25, 183–197,, 2010a. a

Blatman, G. and Sudret, B.: Efficient computation of global sensitivity indices using sparse polynomial chaos expansions, Reliab. Eng. Syst. Safe., 95, 1216–1229,, 2010b. a

Blatman, G. and Sudret, B.: Adaptive sparse polynomial chaos expansion based on least angle regression, J. Comput. Phys., 230, 2345–2367,, 2011a. a, b

Blatman, G. and Sudret, B.: Principal component analysis and Least Angle Regression in spectral stochastic finite element analysis, in: Proc. 11th Int. Conf. on Applications of Stat. and Prob. in Civil Engineering (ICASP11), Zurich, Switzerland, edited by: Faber, M., Köhler, J., and Nishijima, K., 2011b. a

Blatman, G. and Sudret, B.: Sparse polynomial chaos expansions of vector-valued response quantities, in: Proc. 11th Int. Conf. Struct. Safety and Reliability (ICOSSAR'2013), New York, USA, edited by: Deodatis, G., 3245–3252,, 2013. a

Blazina, T., Sun, Y., Voegelin, A., Lenz, M., Berg, M., and Winkel, L. H.: Terrestrial selenium distribution in China is potentially linked to monsoonal climate, Nat. Commun., 5, 4717,, 2014. a

Blazina, T., Läderach, A., Jones, G. D., Sodemann, H., Wernli, H., Kirchner, J. W., and Winkel, L. H.: Marine primary productivity as a potential indirect source of selenium and other trace elements in atmospheric deposition, Environ. Sci. Technol., 51, 108–118,, 2017. a

Burkholder, J., Sander, S., Abbatt, J., Barker, J., Huie, R., Kolb, C., Kurylo, M., Orkin, V., Wilmouth, D., and Wine, P.: Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies, Evaluation No. 18, Tech. Rep. 15-10, Jet Propulsion Laboratory, Pasadena, (last access: 30 August 2019), 2015. a, b

Cheng, B. and Lee, Y.: Rate constant of OH+ OCS reaction over the temperature range 255–483 K, Int. J. Chem. Kinet., 18, 1303–1314,, 1986. a

Combs, G. F.: Selenium in global food systems, Brit. J. Nutr., 85, 517–547,, 2001. a

Conde, J. and Sanz Alaejos, M.: Selenium concentrations in natural and environmental waters, Chem. Rev., 97, 1979–2004,, 1997. a

Cox, R. and Sheppard, D.: Reactions of OH radicals with gaseous sulphur compounds, Nature, 284, 330–331,, 1980. a

Cutter, G. A.: Metalloids in wet deposition on Bermuda: concentrations, sources, and fluxes, J. Geophys. Res.-Atmos., 98, 16777–16786,, 1993. a

Cutter, G. A. and Church, T. M.: Selenium in western Atlantic precipitation, Nature, 322, 720–722,, 1986. a

Cutter, G. A. and Cutter, L. S.: Metalloids in the high latitude North Atlantic Ocean: sources and internal cycling, Mar. Chem., 61, 25–36,, 1998. a

Dasch, J. M. and Wolff, G. T.: Trace inorganic species in precipitation and their potential use in source apportionment studies, Water Air Soil Pollut., 43, 401–412, 1989. a

De Gregori, I., Lobos, M. G., and Pinochet, H.: Selenium and its redox speciation in rainwater from sites of Valparaıso region in Chile, impacted by mining activities of copper ores, Water Res., 36, 115–122,, 2002. a

Dentener, F., Kinne, S., Bond, T., Boucher, O., Cofala, J., Generoso, S., Ginoux, P., Gong, S., Hoelzemann, J., and Ito, A.: Emissions of primary aerosol and precursor gases in the years 2000 and 1750 prescribed data-sets for AeroCom, Atmos. Chem. Phys., 6, 4321–4344,, 2006. a, b, c, d, e

Dinh, Q. T., Cui, Z., Huang, J., Tran, T. A. T., Wang, D., Yang, W., Zhou, F., Wang, M., Yu, D., and Liang, D.: Selenium distribution in the Chinese environment and its relationship with human health: a review, Environ. Int., 112, 294–309,, 2018. a

Egorova, T., Rozanov, E., Zubov, V., and Karol, I.: Model for investigating ozone trends (MEZON), Izv. Atmos. Oceanic Phy., 39, 277–292, 2003. a

Eldred, R. A.: Comparison of Selenium and Sulfur at Remote Sites, J. Air Waste Manage., 47, 204–211,, 1997. a

EMEP/CCC: Co-operative Programme for Monitoring and Evaluation of the Long Range Transmission of Air Pollutants in Europe/Chemical Co-ordinating Centre: Heavy metals and POP measurements reports, available at:, last access: 30 January 2020. 

Environment and Climate Change Canada, Metals in Precipitation, available at:, last access: 30 January 2020. 

Eyring, V., Chipperfield, M. P., Giorgetta, M. A., Kinnison, D. E., Manzini, E., Matthes, K., Newman, P. A., Pawson, S., Shepherd, T. G., and Waugh, D. W.: Overview of the new CCMVal reference and sensitivity simulations in support of upcoming ozone and climate assessments and the planned SPARC CCMVal report, SPARC Newsletter, 30, 20–26, (last access: 9 July 2019), 2008. a

Fairweather-Tait, S. J., Bao, Y., Broadley, M. R., Collings, R., Ford, D., Hesketh, J. E., and Hurst, R.: Selenium in human health and disease, Antiox. Redox Sign., 14, 1337–1383,, 2011. a

Feinberg, A., Moustapha, M., Stenke, A., Sudret, B., Peter, T., and Winkel, L. H. E.: Simulation data used for sensitivity analysis of atmospheric selenium in SOCOL-AER,, 2019a. a

Feinberg, A., Sukhodolov, T., Luo, B.-P., Rozanov, E., Winkel, L. H. E., Peter, T., and Stenke, A.: Improved tropospheric and stratospheric sulfur cycle in the aerosol–chemistry–climate model SOCOL-AERv2, Geosci. Model Dev., 12, 3863–3887,, 2019b. a, b, c

Finn, E. J. and King, G.: The absorption spectrum of carbonyl selenide: Sub-Rydberg transitions, J. Mol. Spectrosc., 56, 39–51,, 1975. a

Floor, G. H. and Román-Ross, G.: Selenium in volcanic environments: a review, Appl. Geochem., 27, 517–531,, 2012. a, b

Fordyce, F. M.: Selenium deficiency and toxicity in the environment, in: Essentials of Medical Geology, Springer, 375–416, 2013. a

Ghanem, R. and Spanos, P.: Stochastic Finite Elements: A Spectral Approach, Courier Dover Publications, Mineola, 2nd Edn., 45–65, 2003. a, b

Goodeve, C. and Stein, N.: The absorption spectra and the optical dissociation of the hydrides of the oxygen group, Trans. Faraday Soc., 27, 393–402,, 1931. a

Goupy, J. and Creighton, L.: Introduction aux plans d'expériences, vol. 3, Dunod Paris, 2006. a

Gratz, L. E., Keeler, G. J., Morishita, M., Barres, J. A., and Dvonch, J. T.: Assessing the emission sources of atmospheric mercury in wet deposition across Illinois, Sci. Total Environ., 448, 120–131,, 2013. a, b

Haygarth, P., Cooke, A., Jones, K., Harrison, A., and Johnston, A.: Long‐term change in the biogeochemical cycling of atmospheric selenium: Deposition to plants and soil, J. Geophys. Res.-Atmos., 98, 16 769–16 776,, 1993. a

Heaton, R. W., Rahn, K. A., and Lowenthal, D. H.: Determination of trace elements, including regional tracers, in Rhode Island precipitation, Atmos. Environ. Pt. A, 24, 147–153,, 1990. a

Horowitz, L. W., Walters, S., Mauzerall, D. L., Emmons, L. K., Rasch, P. J., Granier, C., Tie, X., Lamarque, J., Schultz, M. G., and Tyndall, G. S.: A global simulation of tropospheric ozone and related tracers: Description and evaluation of MOZART, version 2, J. Geophys. Res.-Atmos., 108, D24,, 2003. a

Jacob, D. J.: Introduction to Atmospheric Chemistry, Princeton University Press, 49–50, 1999. a

Jones, G. D., Droz, B., Greve, P., Gottschalk, P., Poffet, D., McGrath, S. P., Seneviratne, S. I., Smith, P., and Winkel, L. H.: Selenium deficiency risk predicted to increase under future climate change, P. Natl. Acad. Sci. USA, 114, 2848–2853,, 2017. a

Kapur, J.: Maximum-entropy Models in Science and Engineering, Wiley, 44–87, 1989. a

Kerkweg, A., Buchholz, J., Ganzeveld, L., Pozzer, A., Tost, H., and Jöckel, P.: An implementation of the dry removal processes DRY DEPosition and SEDImentation in the Modular Earth Submodel System (MESSy), Atmos. Chem. Phys., 6, 4617–4632,, 2006. a

Kerkweg, A., Buchholz, J., Ganzeveld, L., Pozzer, A., Tost, H., and Jöckel, P.: Corrigendum to “Technical Note: An implementation of the dry removal processes DRY DEPosition and SEDImentation in the Modular Earth Submodel System (MESSy)” published in Atmos. Chem. Phys., 6, 4617–4632, 2006, Atmos. Chem. Phys., 9, 9569–9569,, 2009.+ a

King, G. and Srikameswaran, K.: Analysis of the 4500 Å absorption system (R system) of carbon diselenide, J. Mol. Spectrosc., 31, 269–291,, 1969. a

Kurylo, M. J.: Flash photolysis resonance fluorescence investigation of the reactions of OH radicals with OCS and CS2, Chem. Phys. Lett., 58, 238–242,, 1978. a

Låg, J. and Steinnes, E.: Regional distribution of selenium and arsenic in humus layers of Norwegian forest soils, Geoderma, 20, 3–14,, 1978. a

Lamarque, J.-F., Bond, T. C., Eyring, V., Granier, C., Heil, A., Klimont, Z., Lee, D., Liousse, C., Mieville, A., and Owen, B.: Historical (1850–2000) gridded anthropogenic and biomass burning emissions of reactive gases and aerosols: methodology and application, Atmos. Chem. Phys., 10, 7017–7039,, 2010. a, b, c

Lana, A., Bell, T., Simó, R., Vallina, S., Ballabrera‐Poy, J., Kettle, A., Dachs, J., Bopp, L., Saltzman, E., and Stefels, J.: An updated climatology of surface dimethlysulfide concentrations and emission fluxes in the global ocean, Global Biogeochem. Cy., 25, GB1004,, 2011. a, b

Landing, W., Caffrey, J., Nolek, S., Gosnell, K., and Parker, W.: Atmospheric wet deposition of mercury and other trace elements in Pensacola, Florida, Atmos. Chem. Phys., 10, 4867–4877,, 2010. a

Lawson, N. M. and Mason, R. P.: Concentration of mercury, methylmercury, cadmium, lead, arsenic, and selenium in the rain and stream water of two contrasting watersheds in western Maryland, Water Res., 35, 4039–4052,, 2001. a

Le Gratiet, L., Marelli, S., and Sudret, B.: Metamodel-based sensitivity analysis: polynomial chaos expansions and Gaussian processes, chap. 38, Springer International Publishing, cham, Switzerland, 1289–1325, 2017. a

Lee, L., Carslaw, K., Pringle, K., Mann, G., and Spracklen, D.: Emulation of a complex global aerosol model to quantify sensitivity to uncertain parameters, Atmos. Chem. Phys., 11, 12253–12273,, 2011. a, b, c, d

Lee, L., Carslaw, K., Pringle, K., and Mann, G.: Mapping the uncertainty in global CCN using emulation, Atmos. Chem. Phys., 12, 9739–9751,, 2012. a

Lee, L., Pringle, K., Reddington, C., Mann, G., Stier, P., Spracklen, D., Pierce, J., and Carslaw, K.: The magnitude and causes of uncertainty in global model simulations of cloud condensation nuclei, Atmos. Chem. Phys., 13, 8879–8914,, 2013. a

Leu, M.-T. and Smith, R. H.: Kinetics of the gas-phase reaction between hydroxyl and carbonyl sulfide over the temperature range 300–517 K, The J. Phys. Chem., 85, 2570–2575,, 1981. a

Li, S., Chwee, T. S., and Fan, W. Y.: FTIR studies of O(3P) atom reactions with CSe2, SCSe, and OCSe, J. Phys. Chem. A, 109, 11815–11822,, 2005. a, b, c, d, e

Lin, S.-J. and Rood, R. B.: Multidimensional flux-form semi-Lagrangian transport schemes, Mon. Weather Rev., 124, 2046–2070,<2046:MFFSLT>2.0.CO;2, 1996. a

Liu, J. Q., Yang, Y. J., Di, Y., Yang, J., Huang, W. W., Wen, T. X., and Li, Y. W.: Trace Elements of Precipitation in the Shigatse Region, Southern Tibetan Plateau, Adv. Mat. Res., 518–523, 1652–1656, 2012. a

Loeppky, J. L., Sacks, J., and Welch, W. J.: Choosing the sample size of a computer experiment: A practical guide, Technometrics, 51, 366–376,, 2009. a

Lynam, M. M., Dvonch, J. T., Barres, J. A., Morishita, M., Legge, A., and Percy, K.: Oil sands development and its impact on atmospheric wet deposition of air pollutants to the Athabasca Oil Sands Region, Alberta, Canada, Environ. Pollut., 206, 469–478,, 2015. a

Marelli, S. and Sudret, B.: UQLab: A framework for uncertainty quantification in Matlab, in: Vulnerability, uncertainty, and risk: quantification, mitigation, and management, Second International Conference on Vulnerability and Risk Analysis and Management (ICVRAM) and the Sixth International Symposium on Uncertainty, Modeling, and Analysis (ISUMA), 13–16 July, 2554–2563,, 2014. a, b

Marelli, S. and Sudret, B.: UQLab user manual – Polynomial chaos expansions, Tech. Rep., Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, report #UQLab-V1.3-10, 2019. a

Marelli, S., Lamas, C., Konakli, K., Mylonas, C., Wiederkehr, P., and Sudret, B.: UQLab user manual – Sensitivity analysis, Tech. Rep., Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, report #UQLab-V1.2-106, 2019. a, b, c

Marshall, L., Johnson, J. S., Mann, G. W., Lee, L., Dhomse, S. S., Regayre, L., Yoshioka, M., Carslaw, K. S., and Schmidt, A.: Exploring How Eruption Source Parameters Affect Volcanic Radiative Forcing Using Statistical Emulation, J. Geophys. Res.-Atmos., 124, 964–985,, 2019. a

Mason, R. P., Soerensen, A. L., DiMento, B. P., and Balcom, P. H.: The Global Marine Selenium Cycle: Insights From Measurements and Modeling, Global Biogeochem. Cy., 32, 1720–1737,, 2018. a, b

McKay, M. D., Beckman, R. J., and Conover, W. J.: Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, 21, 239–245,, 1979. a

Mosher, B. W. and Duce, R. A.: Vapor phase and particulate selenium in the marine atmosphere, J. Geophys. Res.-Ocean., 88, 6761–6768,, 1983. a

Mosher, B. W. and Duce, R. A.: A global atmospheric selenium budget, J. Geophys. Res.-Atmos., 92, 13289–13298,, 1987. a, b, c, d

Nie, X., Wang, Y., Li, Y., Sun, L., Li, T., Yang, M., Yang, X., and Wang, W.: Characteristics and impacts of trace elements in atmospheric deposition at a high-elevation site, southern China, Environ. Sci. Pollut. Res., 24, 22839–22851,, 2017. a

Nightingale, P. D., Malin, G., Law, C. S., Watson, A. J., Liss, P. S., Liddicoat, M. I., Boutin, J., and Upstill‐Goddard, R. C.: In situ evaluation of air‐sea gas exchange parameterizations using novel conservative and volatile tracers, Global Biogeochem. Cy., 14, 373–387,, 2000. a, b

Nriagu, J. O.: A global assessment of natural sources of atmospheric trace metals, Nature, 338, 47–49,, 1989. a, b

Nriagu, J. O. and Pacyna, J. M.: Quantitative assessment of worldwide contamination of air, water and soils by trace metals, Nature, 333, 134–139,, 1988. a

Pan, Y. and Wang, Y.: Atmospheric wet and dry deposition of trace elements at 10 sites in Northern China, Atmos. Chem. Phys., 15, 951–972,, 2015. a

Pavageau, M.-P., Pecheyran, C., Krupp, E. M., Morin, A., and Donard, O. F.: Volatile metal species in coal combustion flue gas, Environ. Sci. Technol., 36, 1561–1573,, 2002. a, b

Pavageau, M.-P., Morin, A., Séby, F., Guimon, C., Krupp, E., Pécheyran, C., Poulleau, J., and Donard, O. F.: Partitioning of metal species during an enriched fuel combustion experiment. Speciation in the gaseous and particulate phases, Environ. Sci. Technol., 38, 2252–2263,, 2004. a

Pearson, C., Howard, D., Moore, C., and Obrist, D.: Mercury and trace metal wet deposition across five stations in Alaska: controlling factors, spatial patterns, and source regions, Atmos. Chem. Phys., 19, 6913–6929,, 2019. a

Petersen, G. and Munthe, J.: Atmospheric mercury species over central and northern Europe, Model calculations and comparison with observations from the Nordic air and precipitation network for 1987 and 1988, Atmos. Environ., 29, 47–67,, 1995. a

Pham, M., Müller, J.-F., Brasseur, G. P., Granier, C., and Mégie, G.: A three-dimensional study of the tropospheric sulfur cycle, J. Geophys. Res.-Atmos., 100, 26061–26092,, 1995. a

Rael, R. M., Tuzaon, E. C., and Frankenberger Jr., W. T.: Gas-phase reactions of dimethyl selenide with ozone and the hydroxyl and nitrate radicals, Atmos. Environ., 30, 1221–1232,, 1996. a, b

Ramankutty, N., Evan, A. T., Monfreda, C., and Foley, J. A.: Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000, Global Biogeochem. Cy., 22, GB1003,, 2008. a, b

Rayner, N., Parker, D. E., Horton, E., Folland, C., Alexander, L., Rowell, D., Kent, E., and Kaplan, A.: Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century, J. Geophys. Res.-Atmos., 108, ACL 2-1 to ACL 2-29,, 2003. a

Revell, L. E., Stenke, A., Tummon, F., Feinberg, A., Rozanov, E., Peter, T., Abraham, N. L., Akiyoshi, H., Archibald, A. T., and Butchart, N.: Tropospheric ozone in CCMI models and Gaussian process emulation to understand biases in the SOCOLv3 chemistry–climate model, Atmos. Chem. Phys., 18, 16155–16172,, 2018. a

Roeckner, E., Bäuml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Rhodin, A., Schlese, U., Schulzweida, U., and Tompkins, A.: The atmospheric general circulation model ECHAM 5, Part I: Model description, Tech. Rep. 349, Max-Planck-Institut für Meteorologie, Hamburg, (last access: 9 July 2019), 2003. a

Ross, H. B.: An atmospheric selenium budget for the region 30 N to 90 N, Tellus B, 37, 78–90,, 1985. a, b, c, d

Rumble, J.: CRC Handbook of Chemistry and Physics, CRC Press, chap. 6: Fluid Properties, chap. 9: Molecular Structure and Spectroscopy, 2017. a, b

Ryan, E., Wild, O., Voulgarakis, A., and Lee, L.: Fast sensitivity analysis methods for computationally expensive models with multi-dimensional output, Geosci. Model Dev., 11, 3131–3146,, 2018. a, b, c

Sakata, M., Marumoto, K., Narukawa, M., and Asakura, K.: Regional variations in wet and dry deposition fluxes of trace elements in Japan, Atmos. Environ., 40, 521–531,, 2006. a

Saltelli, A., Tarantola, S., and Campolongo, F.: Sensitivity anaysis as an ingredient of modeling, Stat. Sci., 15, 377–395, 2000. a

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S.: Global Sensitivity Analysis: the Primer, John Wiley & Sons, 11–42, 159–160, 2008. a, b, c

Savage, L., Carey, M., Hossain, M., Islam, M. R., de Silva, P. M. C., Williams, P. N., and Meharg, A. A.: Elevated trimethylarsine oxide and inorganic arsenic in northern hemisphere summer monsoonal wet deposition, Environ. Sci. Technol., 51, 12210–12218,, 2017. a

Savage, L., Carey, M., Williams, P. N., and Meharg, A. A.: Maritime Deposition of Organic and Inorganic Arsenic, Environ. Sci. Technol., 53, 7288–7295,, 2019. a

Schultz, M. and Rast, S.: Emission data sets and methodologies for estimating emissions, REanalysis of the TROpospheric chemical composition over the past 40 years, A long-term global modeling study of tropospheric chemistry funded under the 5th EU framework programme, Tech. Rep., EVK2-CT-2002-00170, 2007. a

Scudlark, J. R., Conko, K. M., and Church, T. M.: Atmospheric wet deposition of trace elements to Chesapeake Bay: CBAD study year 1 results, Atmos. Environ., 28, 1487–1498,, 1994. a

Scudlark, J. R., Rice, K. C., Conko, K. M., Bricker, O. P., and Church, T. M.: Transmission of atmospherically derived trace elements through an undeveloped, forested Maryland watershed, Water Air Soil Pollut., 163, 53–79,, 2005. a

Seinfeld, J. H. and Pandis, S. N.: Atmospheric Chemistry and Physics: from Air Pollution to Climate Change, John Wiley & Sons, 19–23, 71–74, 829–851, 856–873, 2016. a, b, c, d

Sheng, J., Weisenstein, D. K., Luo, B., Rozanov, E., Stenke, A., Anet, J., Bingemer, H., and Peter, T.: Global atmospheric sulfur budget under volcanically quiescent conditions: Aerosol‐chemistry‐climate model predictions and validation, J. Geophys. Res.-Atmos., 120, 256–276,, 2015. a, b, c, d

Shimamura, T., Iwashita, M., Iijima, S., Shintani, M., and Takaku, Y.: Major to ultra trace elements in rainfall collected in suburban Tokyo, Atmos. Environ., 41, 6999–7010,, 2007. a

Smith, S. J., Aardenne, J. v., Klimont, Z., Andres, R. J., Volke, A., and Delgado Arias, S.: Anthropogenic sulfur dioxide emissions: 1850–2005, Atmos. Chem. Phys., 11, 1101–1116,, 2011. a, b, c, d

Sobol, I. M.: Sensitivity estimates for nonlinear mathematical models, Mathematical Modelling and Computational Experiments, 1, 407–414, 1993. a, b

Stenke, A., Schraner, M., Rozanov, E., Egorova, T., Luo, B., and Peter, T.: The SOCOL version 3.0 chemistry–climate model: description, evaluation, and implications from an advanced transport algorithm, Geosci. Model Dev., 6, 1407–1427,, 2013. a

Stohl, A., Seibert, P., Arduini, J., Eckhardt, S., Fraser, P., Greally, B. R., Lunder, C., Maione, M., Mühle, J., O'Doherty, S., Prinn, R. G., Reimann, S., Saito, T., Schmidbauer, N., Simmonds, P. G., Vollmer, M. K., Weiss, R. F., and Yokouchi, Y.: An analytical inversion method for determining regional and global emissions of greenhouse gases: Sensitivity studies and application to halocarbons, Atmos. Chem. Phys., 9, 1597–1620,, 2009. a

Sudret, B.: Global sensitivity analysis using polynomial chaos expansions, Reliab. Eng. Syst. Safety, 93, 964–979,, 2008. a, b, c

Suess, E., Aemisegger, F., Sonke, J. E., Sprenger, M., Wernli, H., and Winkel, L. H. E.: Marine versus Continental Sources of Iodine and Selenium in Rainfall at Two European High-Altitude Locations, Environ. Sci. Technol., 53, 1905–1917,, 2019. a, b

Suzuki, Y., Sugimura, Y., and Miyake, Y.: The content of selenium and its chemical form in rain water and aerosol in Tokyo, J. Meteorol. Soc. Jpn. Ser. II, 59, 405–409,, 1981. a

Symonds, R. B. and Reed, M. H.: Calculation of multicomponent chemical equilibria in gas-solid-liquid systems: Calculation methods, thermochemical data, and applications to studies of high-temperature volcanic gases with examples from Mt. St. Helens, Am. J. Sci., 293, 758–864,, 1993. a

Symonds, R. B., Rose, W. I., Reed, M. H., Lichte, F. E., and Finnegan, D. L.: Volatilization, transport and sublimation of metallic and non-metallic elements in high temperature gases at Merapi Volcano, Indonesia, Geochim. Cosmochim. Ac., 51, 2083–2101,, 1987. a

Thompson, K., Canosa-Mas, C., and Wayne, R.: Kinetics and mechanism of the reaction between atomic chlorine and dimethyl selenide; comparison with the reaction between atomic chlorine and dimethyl sulfide, Phys. Chem. Chem. Phys., 4, 4133–4139,, 2002. a

Tost, H., Jöckel, P., Kerkweg, A., Sander, R., and Lelieveld, J.: A new comprehensive SCAVenging submodel for global atmospheric chemistry modelling, Atmos. Chem. Phys., 6, 565–574,, 2006. a

Uchiyama, R., Okochi, H., Ogata, H., Katsumi, N., and Nakano, T.: Characteristics of trace metal concentration and stable isotopic composition of hydrogen and oxygen in “urban-induced heavy rainfall” in downtown Tokyo, Japan; The implication of mineral/dust particles on the formation of summer heavy rainfall, Atmos. Res., 217, 73–80,, 2019. a

Wackerly, D., Mendenhall, W., and Scheaffer, R. L.: Mathematical Statistics with Applications, Cengage Learning, 593–596, 2014. a

Wahner, A. and Ravishankara, A.: The kinetics of the reaction of OH with COS, J. Geophys. Res.-Atmos., 92, 2189–2194,, 1987. a

Wai, K.-M., Wu, S., Li, X., Jaffe, D. A., and Perry, K. D.: Global atmospheric transport and source-receptor relationships for arsenic, Environ. Sci. Technol., 50, 3714–3720,, 2016. a, b

Wardell, L., Kyle, P., and Counce, D.: Volcanic emissions of metals and halogens from White Island (New Zealand) and Erebus volcano (Antarctica) determined with chemical traps, J. Volcanol. Geoth. Res., 177, 734–742,, 2008. a

Watts, S. F.: The mass budgets of carbonyl sulfide, dimethyl sulfide, carbon disulfide and hydrogen sulfide, Atmos. Environ., 34, 761–779,, 2000. a, b, c

Weisenstein, D. K., Yue, G. K., Ko, M. K., Sze, N.-D., Rodriguez, J. M., and Scott, C. J.: A two-dimensional model of sulfur species and aerosols, J. Geophys. Res.-Atmos., 102, 13019–13035,, 1997. a

Weller, R., Wöltjen, J., Piel, C., Resenberg, R., Wagenbach, D., König-Langlo, G., and Kriews, M.: Seasonal variability of crustal and marine trace elements in the aerosol at Neumayer station, Antarctica, Tellus B, 60, 742–752,, 2008. a

Wen, H. and Carignan, J.: Reviews on atmospheric selenium: emissions, speciation and fate, Atmos. Environ., 41, 7151–7165,, 2007. a, b, c, d, e, f, g, h, i

Wesely, M.: Parameterization of surface resistances to gaseous dry deposition in regional-scale numerical models, Atmos. Environ., 23, 1293–1304,, 1989.  a

Winkel, L., Vriens, B., Jones, G., Schneider, L., Pilon-Smits, E., and Bañuelos, G.: Selenium cycling across soil-plant-atmosphere interfaces: a critical review, Nutrients, 7, 4199–4239,, 2015. a, b

Winkel, L. H. E., Johnson, C. A., Lenz, M., Grundl, T., Leupin, O. X., Amini, M., and Charlet, L.: Environmental Selenium Research: From Microscopic Processes to Global Understanding, Environ. Sci. Technol., 46, 571–579,, 2012. a

Xing, J., Song, J., Yuan, H., Wang, Q., Li, X., Li, N., Duan, L., and Qu, B.: Atmospheric wet deposition of dissolved trace elements to Jiaozhou Bay, North China: Fluxes, sources and potential effects on aquatic environments, Chemosphere, 174, 428–436,, 2017. a

Xiu, D. and Karniadakis, G. E.: The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing, 24, 619–644,, 2002. a

Yan, R., Gauthier, D., and Flamant, G.: Possible interactions between As, Se, and Hg during coal combustion, Combust. Flame, 120, 49–60,, 2000. a

Yan, R., Gauthier, D., Flamant, G., and Wang, Y.: Behavior of selenium in the combustion of coal or coke spiked with Se, Combust. Flame, 138, 20–29,, 2004. a

Zhou, J., Wang, Y., Yue, T., Li, Y., Wai, K.-M., and Wang, W.: Origin and distribution of trace elements in high-elevation precipitation in southern China, Environ. Sci. Pollut. Res., 19, 3389–3399,, 2012. a

Short summary
The amount of the micronutrient selenium in food largely depends on the amount and form of selenium in soil. The atmosphere acts as a source of selenium to soils through deposition, yet little information is available about atmospheric selenium cycling. Therefore, we built the first global atmospheric selenium model. Through sensitivity and uncertainty analysis we determine that selenium can be transported thousands of kilometers and that measurements of selenium emissions should be prioritized.
Final-revised paper