The Superstatistical Nature and Interoccurrence Time of Atmospheric Mercury
Concentration Fluctuations
F. Carbone1 , A. G. Bruno1,2, A. Naccarato1 , F. De Simone1 , C. N. Gencarelli1 ,
F. Sprovieri1 , I. M. Hedgecock1 , M. S. Landis3 , H. Skov4 , K. A. Pfaffhuber5 , K. A. Read6, L. Martin7 , H. Angot8,9, A. Dommergue8 , O. Magand8, and N. Pirrone10
1CNR-Institute of Atmospheric Pollution Research, Division of Rende, UNICAL-Polifunzionale, Rende, Italy,2Dipartimento di Fisica, Università della Calabria, Rende, Italy,3U. S. Environmental Protection Agency, Office of Research and
Development, Research Triangle Park, Durham, NC, USA,4Department of Environmental Science, Aarhus University, Roskilde, Denmark,5Norwegian Institute for Air Research (NILU), Kjeller, Norway,6NCAS, National Centre for Atmospheric Sciences, University of York, York, UK,7Cape Point GAW Station, Climate and Environment Research and Monitoring, South African Weather Service, Stellenbosch, South Africa,8University Grenoble Alpes, CNRS, IRD, IGE, Grenoble, France,
9Now at Institute for Data, Systems, and Society, Massachusetts Institute of Technology, Cambridge, MA, USA,
10CNR-Institute of Atmospheric Pollution Research, Area della Ricerca di Roma 1, Monterotondo, Italy
Abstract
The probability density function (PDF) of the time intervals between subsequent extreme events in atmospheric Hg0concentration data series from different latitudes has been investigated. The Hg0 dynamic possesses a long-term memory autocorrelation function. Above a fixed thresholdQin the data, the PDFs of the interoccurrence time of the Hg0data are well described by a Tsallis q-exponential function.This PDF behavior has been explained in the framework of superstatistics, where the competition between multiple mesoscopic processes affects the macroscopic dynamics. An extensive parameter𝜇, encompassing all possible fluctuations related to mesoscopic phenomena, has been identified. It follows a𝜒2distribution, indicative of the superstatistical nature of the overall process. Shuffling the data series destroys the long-term memory, the distributions become independent ofQ, and the PDFs collapse on to the same exponential distribution. The possible central role of atmospheric turbulence on extreme events in the Hg0data is highlighted.
1. Introduction
A number of studies using different methodologies have shown that long-term memory in atmospheric pol- lutant concentrations exists (Chelani, 2016; Lovejoy & Schertzer, 2013; Tuck, 2010), that is, up to a limit the concentrations maintain a certain correlation over time. Usually a long-term memory process is defined by a strong coupling between measured values at different time lags,𝓁, and the system’s dynamics are char- acterized by the presence of complex mesoscopic spatiotemporal patterns. These patterns are associated with the generation of high-amplitude fluctuations over a broad range of spatial and temporal scales giving rise to scale-free relationships for statistical quantities (Frisch, 1995; McComb, 1990; Monin & Yaglom, 2007).
These mesoscopic processes occur within macroscopic phenomena, and their behavior evolves into a power law decay of the autocorrelation function. Conversely, in a short-term memory process the autocorrelation function decreases exponentially or to zero after a certain time,𝓁. The dynamics of pollutant concentration variations depend on numerous processes, (for a review see Chelani, 2016); however, due to their complexity it is not possible to precisely describe their behavior and properties over space and time. One of the prin- cipal characteristics of complex dynamical systems is the intermittency (Briggs & Beck, 2007; Carbone &
Sorriso-Valvo, 2014; Carbone, Gencarelli, et al., 2016; Manshour et al., 2016; Warhaft, 2000). Intermittency rep- resents the strongly correlated fluctuations that lead to deviations from a normal probability distribution func- tion (PDF). In the atmospheric boundary layer, intermittency is an important part of a continuous spectrum of atmospheric motions (Katul et al., 2006; Vindel & Yagüe, 2011; Wyngaard, 1992). Within large-scale temporal variations of atmospheric motion, fluctuations in pollutant species concentrations result from interactions of a large ensemble of mesoscopic phenomena, occurring contemporaneously in the atmosphere: turbulence (Wyngaard, 1992), variation in anthropogenic and natural emission sources (Carbone, Landis, et al., 2016;
RESEARCH ARTICLE
10.1002/2017JD027384
Key Points:
• Superstatistics of atmospheric Hg0concentration extremes is demonstrated via the probability density functions of the interoccurrence time
• The statistics of interoccurrence times has been related to the long-term memory of Hg0concentrations via Tsallis q statistics
• Universal scaling in atmospheric Hg0concentration extremes was found, a possible link to atmospheric turbulence is suggested
Correspondence to:
F. Carbone, f.carbone@iia.cnr.it
Citation:
Carbone, F., Bruno, A. G., Naccarato, A., De Simone, F., Gencarelli, C. N., Sprovieri, F.,…Pirrone, N. (2018).
The superstatistical nature and interoccurrence time of atmospheric mercury concentration fluctuations.Journal of Geophysical Research: Atmospheres,123, 764–774.
https://doi.org/10.1002/2017JD027384
Received 30 JUN 2017 Accepted 30 NOV 2017
Accepted article online 6 DEC 2017 Published online 19 JAN 2018
©2017. American Geophysical Union.
All Rights Reserved.
Pirrone et al., 2010), variation in deposition velocity, loss through chemical reactions that is in turn deter- mined by fluctuating reactant/oxidant concentrations, and other extreme events. In the specific case of Hg0 these extreme events would include phenomena such as convective storms, forest fires, and atmospheric Hg0depletion events (De Simone et al., 2017; Dvonch et al., 2005; Holmes et al., 2016; Lindberg et al., 2002;
Schroeder et al., 1998). Understanding the dynamics of these emergent extreme events, meteorological, chemical, and anthropological, represents the key to understanding complex dynamical systems.
Due to their complexity, the analysis of these systems has focused on the understanding of where certain features are exhibited by a large class of phenomena, regardless of the details of their structure (Manshour et al., 2016). Beck and Cohen (2003) have shown that complex nonequilibrium systems, which possess a spa- tiotemporally fluctuating intensive quantity, can often be effectively described by a superposition of statistics, so-called “superstatistics.” The core idea is to characterize the system under consideration as a superposition of several statistics, which act on different time scales (Beck et al., 2005; Beck, 2010). Importantly for this study, it has been shown that a suitable intensive parameter𝜇of the complex system, (e.g., a local variance parameter extracted from an experimental data set) exists, which generates a Tsallis q-distribution if𝜇is a 𝜒2-distributed random variable (Abe & Okamoto, 2001; Briggs & Beck, 2007; Manshour et al., 2016; Tsallis, 1988; Tsallis et al., 1998). Therefore, an understanding of sudden (intermittent) atmospheric events could pro- vide useful information about climatological phenomena, teleconnections, and atmospheric transport pro- cesses (Bunde et al., 2005). Although the atmosphere is a relatively minor reservoir of Hg0compared to other environmental compartments, it is an important pathway by which Hg0is distributed globally over relatively short timescales, once thought to be roughly a year, but more recently closer to 6 months (Horowitz et al., 2017; Schroeder & Munthe, 1998). Efforts to measure Hg concentrations in monitoring networks worldwide and interpret these data with models have increased recently, seeking to elucidate the way Hg cycles in the environment (Pirrone et al., 2010). A number of important issues remain unexplained due to the com- plex interactions of Hg species, with and within, a multiphase atmospheric matrix that is in a constant state of flux. The focus here is on the long-term memory of atmospheric Hg0concentrations observed at dif- ferent latitudes and their sudden or extreme (intermittent) events through analysis of the interoccurrence times (IOTs) of Hg0concentration fluctuations. The superstatistical nature of these extreme events has been investigated and demonstrated.
1.1. Hg0Measurements Methods
The atmospheric Hg0data used in this analysis were all obtained using automated Tekran (Toronto, Canada) Model 2537A cold vapor atomic fluorescence spectrometer instruments. Tekran equipment provides a detec- tion limit below 0.1 ng/m3and a linear response over the range 1–200 ng/m3within 2%. The instruments are calibrated periodically using the internal permeation source in accordance with the Global Mercury Obser- vation System (GMOS) standard operating procedure (every 72 h with a permeation time of 120 s). Details of the instrument operating parameters can be found in Carbone, Landis, et al. (2016), Landis et al. (2002), Sprovieri et al. (2016), and Steffen et al. (2015). An exhaustive description of most of the sampling sites includ- ing their location, altitude, and climatology, as well as the Hg0data quality assurance/quality control (QA/QC) protocols, for data quality assessment, can be found in Sprovieri et al. (2016). The original data set had a temporal resolution ofΔt=300s; however, some data were excluded by the QA/QC procedure employed in order to ensure data quality. There, the temporal resolution was reduced toΔt=3,600s; this allowed a good compromise between the sample length and data set sampling at different thresholds,Q.
2. Superstatistics and Interoccurrence Times
Superstatistics seeks to represent a complex nonequilibrium system as a superposition of two (or several) statistics, described by an intensive parameter𝜇that fluctuates on a relatively large spatiotemporal scale (Beck & Cohen, 2003; Briggs & Beck, 2007). Thus,𝜇might itself be a stochastic variable that incorporates all the possible fluctuations (a local emergent dynamic or competition in the mesoscopic subparts of the system), which produces a sudden variation in the collective dynamic (Beck, 2010). If𝜇is distributed according to a probability densityf(𝜇), then the long-term associated marginal probability of the physical processP(𝜏) (𝜏 being the time lag between two subsequent events) may then defined as a mixture of exponential distributions within which𝜇fluctuates (Beck, 2010; Manshour et al., 2016; Tsallis & Souza, 2003):
P(𝜏) =
∫
∞
0 𝜇f(𝜇)e−𝜇𝜏d𝜇 (1)
As defined by Beck and Cohen (2003),f(𝜇)cannot be a generic function but must be a normalized proba- bility density, a physically relevant density from statistics (for instance, Gaussian, uniform, chi-squared, log normal but potentially others also), and must be normalizable (the integral∫0∞P(𝜏)d𝜏must exist). The super- statistical parameter𝜇does not necessarily have to be a variable such as inverse temperature, but it could be an effective parameter in a stochastic differential equation, a volatility measure in finance or simply a local variance parameter extracted from an experimental time series, as in the case described here (Beck, 2010).
In this last case,𝜇can be extracted from an experimental data series by partitioning the original data set in Mnonoverlapping windows of fixed size and taking the inverse variance of the data contained within each window (Kosun & Ozdemir, 2016; Rabassa & Beck, 2015).
When a process is random or uncorrelated, the PDF collapses on to an exponential distribution, and in general, the zeroth-order theoretical model for the distributionP(𝜏)can be written asP(𝜏) =𝜇e−𝜇𝜏(Briggs & Beck, 2007;
Manshour et al., 2016; van Kampen, 1981). In that case, equation (1) reduces to the exponential model if there are no fluctuations of the intensive parameter𝜇, and the distributionf(𝜇)is a delta function. Any deviation of f(𝜇)from a delta function yields a nonexponentialP(𝜏).
Iff(𝜇)follows a𝜒2distribution (equation (2))
f(𝜇) ∼𝜇k∕2−1exp [
−k𝜇 2𝜇0
]
, (2)
withkdegrees of freedom (𝜇0is a constant), the corresponding superstatistics, obtained by integrating over all𝜇, is described by aqstatistic (Abe & Okamoto, 2001; Tsallis, 1988; Tsallis et al., 1998).
Hence, the behavior of the PDF is described by a Tsallisq-exponential function (equation (3)) and possesses asymptotic power laws (Briggs & Beck, 2007; Manshour et al., 2016):
P(𝜏) = 𝛼
[1+𝛽(q−1)𝜏]1∕(q−1), (3)
where𝛼is a normalization factor (Beck, 2010; Wilk & Włodarczyk, 2000),qis a measure of the deviation from an exponential distribution, andq>1indicates a long-tailed distribution. A number of authors have shown that the limit of validity for the parameterqlies in the range 1≤q≤2. The upper limit arises from the normal- ization condition ofP(𝜏)to the unit area and the requirement that the normalization constant,𝛼, is positive (Briggs & Beck, 2007; Douglas et al., 2006; Tsallis, 1988; Wilk & Włodarczyk, 2000). The behavior of theqexpo- nential is principally related to a long-term memory process and also to the presence of strong ramp-cliff or extreme events in the data. The atmosphere is a complex system, described by a large number of variables, which are nonlinearly coupled by competing physical and chemical processes. The competition among pro- cesses becomes evident once the observations of a single variable are dominated by sudden and intermittent fluctuations. Variables can be either active or passive, transported by the atmospheric flow, or feeding back (and thus modifying) the flow itself (Celani et al., 2004; Mazzitelli & Lanotte, 2012).
This study focuses on fluctuations in Hg0concentration, which is hypothesized to be governed by a com- bination of large-scale atmospheric flow and mesoscopic (small-scale) processes. The aim is to investigate whether the statistics of the Hg0data series can provide insights into the similarities or differences between the measurement sites. Such an investigation requires that the data set is suitable for the proposed analysis.
Rabassa and Beck (2015) lay out a series of criteria that may be used in this context. These criteria concern the skewness and kurtosis of the data set, the existence of an appropriate time scale separation, and the rea- sonableness of a local Gaussian approximation. Applying the strategy proposed in Van der Straeten and Beck (2009), the short time scaletShas been found to be betweents∈ [6÷11]h, while the long time scale found is tL>50h. The ratio oftS∕tL≈0.2is comparable with the case presented in Rabassa and Beck (2015). As a further check on the assumption that during long-range transport, Hg0concentration may be considered a real pas- sive scalar quantity, and the Hurst exponent can be evaluated. This can then be compared to the value of H=5∕9predicted by generalized scale invariance, see Tuck (2010). The average Hurst exponent extracted from the data series used here is in good agreement with the valueH=5∕9(Istas & Lang, 1997; Tuck, 2010):
Troll (TRL 72∘S, 2∘E)H=0.58±0.02, Mauna Loa, Hawaii (MLO 19∘N, 155∘W)H=0.49±0.10, and Ny Alesund (NYA 78∘N, 11∘E)H=0.53±0.01.
Three sample Hg0concentration data sets are shown in Figure 1, for Villum Research Station at Station Nord, Greenland (81∘N, 17∘W), Mauna Loa, Hawaii, and Troll, Antarctica.
Figure 1.Hg0data measured at three different sites: (left) Villum Research Station at Station Nord (VRS81∘N,17∘W), (middle) Mauna Loa Observatory (MLO19∘N, 155∘W), and (right) Troll (TRL72∘S,2∘E).
The intermittent nature of the Hg0data is characterized by alternating periods of strong fluctuations and smoother periods characterized by smaller fluctuations. This behavior is determined by a large ensemble of mesoscopic phenomena occurring in the atmosphere. Despite this intermittent dynamic, the phenomenon preserves its non-Markovian nature, since a power law decay in the autocorrelation function is observed (Figure 2) (Alder & Wainwright, 1970; Bunde et al., 2005; Chelani, 2016; Schertzer & Lovejoy, 1985, 1987). Ifx(t) represents the instantaneous concentration of Hg0at timet, the associated autocorrelation function may be written as
C(𝓁) = 1 𝜎x2(T−𝓁)
T−𝓁
∑
t=1
(x(t) −⟨x⟩)(x(t+𝓁) −⟨x⟩) ∼𝓁−𝛾, (4)
where𝓁represents the time lag, and𝜎xis the standard deviation ofx(t). Due to the passive scalar nature of Hg0, the standard deviation of the concentration is related to atmospheric eddy diffusivity and represents a measure of the characteristic width of the plume (Hayley et al., 2002).
All the stations in a given hemisphere present the same scaling, and those close to the equator (tropics) present a scaling faster than those at higher latitudes,𝛾E≈0.4,𝛾N≈0.3, and𝛾S≈0.1, respectively, for tropics, Northern Hemisphere, and Southern Hemisphere. It appears thatC(𝓁)within either of the three zones is inde- pendent of latitude. A slight difference occurs in the tropical stations, the Cape Verde Observatory (CVO) and
Figure 2.Power law decay of Hg0data with slopes𝛾N=0.3,𝛾E=0.4, and𝛾S=0.1, respectively for the Northern Hemisphere, tropics, and Southern Hemisphere. The curves have been vertically shifted for clarity.
MLO, at small scale (small𝓁). At larger scales (longer𝓁) the scaling is in perfect agreement for both stations𝛾E≈0.4. The curvesC(𝓁)in Figure 2 have been vertically shifted for clarity.
The presence of long-range correlations suggest that there might be some other fundamental process (or processes) embedded in the temporal evolu- tion of the Hg0data series, as seen with other pollutants (Chelani, 2016).
The greater the irregularity of the mesoscale processes occurring in the atmo- sphere the faster the system “loses its memory.” In this case these processes would include sources, sinks such as dry deposition or wet scavenging, and chemical transformations. One possible hypothesis to explain the differences in the slopes in the three latitude zones may be the different characteristics of Hg0 emissions in each. Hg0 emission in the tropics includes biomass burn- ing that is irregular, another factor may be the different characteristics of Hg0 emissions in each. Hg0emission in the tropics includes biomass burning that is irregular, another factor may be that meteorological phenomena also tend to be of shorter duration (De Simone et al., 2015). The Northern Hemisphere is the home to most of the world’s anthropogenic emission sources and also has biomass burning events. The Southern Hemisphere in comparison has less variable and lower Hg0emissions, and synoptic rather than local weather conditions often dominate in large regions.
Figure 3.Sample illustration of IOT,𝜏i, obtained from a synthetic data set. Horizontal dashed lines represent the selected thresholdQin the data. By increasingQ, two sets of IOT𝜏1,2can be identified in the data, characterized by an increasing average⟨𝜏i⟩and standard deviation𝜎𝜏, respectively.
An optimal method used to obtain information concerning extreme events in physical or chemical processes is IOT,𝜏, series analysis. IOTs are a measure of the time between the occurrence of two or more subsequent events in the data that exceed a fixed, thresholdQi(Figure 3). The events exceedingQare defined as rare or more usually extreme. If a long-range correlation exists in the data, then the IOTs,𝜏i, are also long-range corre- lated (Bogachev & Bunde, 2008; Bogachev et al., 2007; Eichner et al., 2007; Ferri et al., 2010, 2012; Santhanam
& Kantz, 2005).
Before performing the IOT analysis, the standard normalization procedure, (subtracting the mean value of the data and dividing by the standard deviation) was applied to the Hg0data. Following this strategy, performing the analysis at a generic thresholdQimplies performing the analysis at a fixed standard deviation value of the data. This procedure is required in order to facilitate the comparison of concentration data measured at dif- ferent latitudes. For everyQ, an average⟨𝜏⟩and standard deviation𝜎𝜏are defined, and by increasingQ,⟨𝜏⟩
and𝜎𝜏become larger. The higher theQ, the rarer or more extreme are the events. Also, there is a one-to-one correspondence betweenQand the⟨𝜏⟩,𝜎𝜏 values (Figure 4, middle and right columns) (Chelani, 2016).
The differences between the three normalized data sets are shown in Figure 4, second row. At MLO the
Figure 4.(first row) Dependence on the thresholdQof the parametersNtotQ ,𝜎𝜏, and the average⟨𝜏⟩of the IOT, measured at three very different latitudes:
Mauna Loa (MLO19∘N,155∘W), Troll (TRL,72∘S,2∘E), and Ny Alesund (NYA78∘N,11∘E). (second row) Comparison of the three normalized data set.
dynamic is characterized by a large number of spikes whereQ>0, which can be related to volcanic (Landis et al., 2013), oceanic (Carbone, Landis, et al., 2016), Hg0emissions or to long-range atmospheric transport phenomena. TRL and NYA data show a very different dynamic, characterized by strong fluctuations at both positive and negativeQ. The negative fluctuation is related to Hg depletion events (Steffen et al., 2008).
However, in the rangeQ∈ [0÷1.5]the data set from each site demonstrates the same behavior, characterized by a wide distribution of IOT with the same⟨𝜏⟩and𝜎𝜏(Figure 4, first row). The strong departure of the value of 𝜎𝜏and⟨𝜏⟩for the NYA data set is related to the lower number of positive (Q>0) fluctuations (Figure 4, bottom right). In this case, by increasing the thresholdQ>1(dotted line in Figure 4, second row) the IOT duration becomes longer, but at the same time the number of IOTNtotQ decreases due to the distance between two subsequent peaks, resulting in a rapid increase in𝜎𝜏.
Investigating a real passive scalar such as Hg0is also interesting from a fundamental point of view, because its dynamics are mainly, if not completely, related to turbulent eddies in the atmosphere.
3. Discussion
The PDFs of the IOT,P(𝜏), from the normalized Hg0data were evaluated by using a number of bins,Nbins(Nbins∈ [10÷15]), depending on the specific data set; this choice being principally due to the temporal resolution of the data (Δt=3,600s). For every thresholdQ, the range of the PDF bins was set between the minimum and maximum values of𝜏in logarithmically spaced bins of lengthAbins. However, as it is the PDFs that are calculated, the number of bins and their width are irrelevant to the analysis. The histogram of the IOT,NQ(𝜏), has been used to calculate the PDFP(𝜏) =NQ(𝜏)∕(AbinsNtotQ ), where c is the total number of observations (IOT) at a fixed thresholdQ(Figure 4, top left). Figure 4 shows the standard deviation𝜎𝜏and the average⟨𝜏⟩of the IOT for three very different latitudes. All the experimental data sets demonstrate universal scaling for these two quantities over the rangeQ∈ [0÷1.5]. For this reason all the following analyses were performed using this range ofQ.
The sampling uncertainty on each bin ofP(𝜏)is calculated from the statistical (Poisson) error on the histogram Err[N(x)] =√
N(x), which, following the error on the probability density, becomes Err[P(𝜏)] =
√
P(𝜏)∕(AbinsNtotQ ).
Outside a given range, roughlyQ≥2.5÷3(related to the size of the Hg0data set), the statistics are too poor to permit accurate analysis, since the relative statistical error becomes too large.
Due to the strong fluctuations a thresholdQ>0.5(site dependent) is sufficient to observe an interesting dynamic. At this value a large number of IOTs can be identified, characterized by𝜎𝜏 =15±1days (Figure 6).
Figure 6 shows the PDFs, for three different sites, of the Hg0 IOT extremes (log-log plot). The plots have been shifted vertically for clarity. It is worth noting that for eachQ, the values of𝜏can be drastically differ- ent. Rather thanP(𝜏), in order to compare the different distributions of the IOT, the PDF of the normalized IOTP(𝜏∕𝜎𝜏) was evaluated. This procedure ensures that for everyQthe values of𝜏∕𝜎𝜏 always lie in the same interval. To check the universality of theqexponential, the analysis has been performed at different values ofQfor each normalized data set. Figure 5 shows the PDFP(𝜏∕𝜎𝜏)(symbols), obtained for differ- ent sites at various thresholdsQ. The same figure also shows the relativeqexponential fit obtained from equation (3). The fitting procedure was performed by minimizing the𝜒2statistic by varying the model param- eters (equation (3)):qwithin the closed intervalq∈ [1÷2], while the other parameters in the open interval 𝛼, 𝛽∈ (0÷ ∞). As stated in section 1.1, uncertainties should be minimized to obtain an accurate estimate of P(𝜏). A large number of gaps in the data set could potentially mask the real distribution. Long gaps give raise to a fat-tailed PDF, because the length of a certain number of IOTs can be overestimated, and the exponent qin equation 3 can exceed the theoretical valueq>2. In case of multiple short gaps, especially forQ≈0, the exponentqtends to unity. In this case the theoretical distribution, equation (3), is characterized by a fast decay, and the fitting procedure is unable to fully capture the tail of the PDF.
The value of𝜏ranges from 1.21×104s to 2.42×107s (approximately 3 h to 9 months, the latter being roughly compatible with the average residence time of Hg0in the atmosphere), which in terms of𝜏∕𝜎𝜏is translated into an interval𝜏∕𝜎𝜏∈ [10−2÷101]. The values of𝜎𝜏, used for the normalization, are shown in Figure 6. The behavior of the PDFs illustrates the good agreement between theqexponential (equation (3)) and the data over the whole range of normalized IOTs𝜏∕𝜎𝜏∈ [10−2÷101](Figure 5). The value ofqexponential reported in Figure 5 represent the average obtained over the different thresholdsQreported in the figure.
Figure 5.Log-log plot of PDFP(𝜏∕𝜎𝜏)for the normalized IOT extracted from Hg0data recorded at different latitudes and different thresholdQ. Dashed line is the averagedqexponential fit (equation (3)).
Figure 6.Log-log plot of PDFsP(𝜏∕𝜎𝜏)for the normalized IOT obtained at three different sites and at three different thresholdsQ: VRS81∘N, 17∘W), CVO, and TRL72∘S,2∘E). The plots have been arbitrarily shifted vertically for clarity. For VRS,Q=0.6,1.3,1.5; CVO,Q=0.7,0.9,1; and TRL, Q=0.9,1.1,1.3. For these thresholds all the data possess the same𝜎𝜏. Dashed lines are theqexponential fit (equation (3)).
From Figure 5 it is evident that the normalization process is robust since no difference can be observed in the distributionP(𝜏∕𝜎𝜏)when varying the thresholdQ.
In Figure 6 the PDFsP(𝜏∕𝜎𝜏)have been grouped according to the value of𝜎𝜏: 𝜎𝜏=15±1,𝜎𝜏=20±1, and𝜎𝜏=26±3days, respectively. In this way it may easily be observed that for three very different latitudes: VRS (Arctic), CVO (Tropics), and TRL (Antarctica) all the PDFs converges on the same universal q-exponential distribution 3. The curves have been vertically shifted for clarity.
The normalization factor𝛼(equation (3)) reflects the value𝛼th=𝛽(2−q)as found in Manshour et al. (2016): for MLO atQ=1the values obtained from the q–exponential fit for𝛼=18.7±1,q=1.55±0.1, and𝛽=40±3, using the relation above, give𝛼th=18±3. The same result is obtained for the other stations:
VRSQ=1.2,𝛼=47±1,q=1.53±0.07,𝛽=81±6, and𝛼th=38.1±6.1; CVO Q=1.2,𝛼=17±2,q=1.54±0.12,𝛽=37±7, and𝛼th=17.02±7. Similar agree- ment was found for all the other stations. These results assure the goodness of theqstatistics in describing extreme events in the Hg0time series at dif- ferent latitudes. The PDFs clearly display the same behavior and could be represented by the universal law, equation (3), with the sameqparameter.
The values ofqobtained over all latitudes lie in the rangeq∈ [1.3÷1.7]and are distributed around an average value⟨q⟩≈1.59±0.05(Figure 7). This result makes evident the universality of the mechanism that lies behind the IOTs within the Hg0concentration time series. The results shown in Figure 2, on the other hand, demonstrate the variation with latitude of the local, that is mesoscale, phenomena. Figure 7 shows that the formation of IOTs is a
Figure 7.Average⟨q⟩entQthresholds. All values are distributed around the dashed line⟨q⟩=1.59±0.05. Dash-dotted lines represent the standard deviation of the experimentalqvalues. Hg0data from the following sites were included: TRL72∘S,2∘E), Dumont d’Urville (DDU66∘S,140∘E), Cape Point (CPT)34∘S,18∘E), CVO16∘N,24∘W), Mace Head (MH53∘N,9∘W), Waldhorf (WAL52∘N,10∘E), Ny Alesund (NYA78∘N,11∘E), Andoya (ADY69∘N,16∘E), VRS (81∘N,17∘W), and Alert (ALT82∘N,62∘W).
large-scale characteristic of the time series, and that it is universal and inde- pendent of latitude.
Thisqvalue is significant since a similar value has been reported in a dif- ferent context, namely velocity fluctuations in fully developed turbulence experiments (Manshour et al., 2016), whereq≈1.6was obtained from the IOT statistics. In light of this, it becomes clear that atmospheric turbulence, although it is never fully developed in the atmosphere, potentially plays a central role in the Hg0IOT dynamic.
For comparison Figure 8 shows the PDFs of the normalized IOT𝜏W, obtained from 1 min (top row) and 1 h (bottom row), wind velocity data recorded at Mauna Loa Observatory. The analysis was performed on two different data sets sampled during the years 2015 and 2016, while a third data set covers the period 2012–2016. In these cases theqvalues obtained are in perfect agree- ment with the results from laboratory experiments (Manshour et al., 2016), demonstrating the universality of the phenomenon. Theqvalues found were q=1.58,1.59and 1.61 for the 2015, 2016, and 2012–2016 1 min data sets, and q=1.6,1.55and 1.6 for the 1 h data sets. It should be pointed out that the full data series were used and not selected to use only daytime or nighttime data even though the wind regimes change significantly between day and night on MLO, see Ryan (1997) and Sharma and Barnes (2016).
In light of this result, it becomes evident that the variation fromq≈1.6could be related to other factors that influence the dynamics. For example, at MLO
⟨q⟩=1.67, however, volcanic emissions near MLO are known to emit Hg0and may act as a local perturbation to the large-scale dynamics. Continuous perturbations of this nature could shift the exponentqaway from the dominant dynamic, where the exponent should beq≈1.6.
Figure 8.(top row) A log-log plot of PDFsP(𝜏W∕𝜎W𝜏)for the normalized IOT evaluated from 1 min wind speed data at MLO, for three different data sets.
(bottom row) Log-log plot of PDFs relative to the 1 h wind speed data, the same data sets has been used. Dashed line isq-exponential fit 3.
Figure 9.Log linear plot of PDFsP(𝜏∕𝜎𝜏), for four stations, after the random permutation of the original data. The PDFs are independent of the thresholdQand collapse on the same exponential distribution, exposing their uncorrelated nature. Numbers next to the station label in legend represent theQvalue. The stations reported in the plot are Mace Head (MH53∘N,9∘W), MLO, Cape Point (CPT34∘S,18∘E), TRL, and VRS.
All the other stations show the same behavior but have not been included for figure clarity.
In the case of a nonexponentially decaying PDF, the shape indicates the presence of correlated structures in the process (long-term memory in the phenomenon) (Bogachev & Bunde, 2008; Bogachev et al., 2007; Chelani, 2016;
Eichner et al., 2007; Ferri et al., 2010, 2012). To check the dependence of the correlated structure in the data set, the Hg0data were randomly rearranged over a large number of trials, keeping the random seed fixed. Figure 9 shows the PDF of the shuffled data set in a log linear plot. The random permutation destroys the correlation in the data, and moreover, the distribution becomes independent of the thresholdQand collapses onto the same exponential shape:P(𝜏) =𝜎𝜏−1exp−𝜏∕𝜎𝜏.
In this framework a single measurement can be related to a local emergent dynamic in the system or can encompass all the possible mechanisms capa- ble of producing a sudden variation in the collective dynamic of the system.
Recalling equation (1),𝜇=𝜎𝜏−1is the parameter that incorporates all the pos- sible fluctuations due to competing processes in the mesoscopic dynamics of the system. Since most IOTs are very sharp for large values of𝜇(𝜎𝜏), whereas small values of𝜇correspond to the frequent occurrence of long IOTs, the exis- tence of correlations in the dynamics makes𝜇a fluctuating random variable with a probability densityf(𝜇)(Briggs & Beck, 2007; Manshour et al., 2016).
To evaluate the distributionf(𝜇), the IOTs obtained from each data set were collected in a large number of windowsW𝜇wider than 48 h and𝜇=𝜎−1𝜏 was measured in each window (Rabassa & Beck, 2015; Van der Straeten & Beck, 2009). Their distribution was evaluated for both the original and shuffled data sets. Figure 10 (left and middle) shows the distribution off(𝜇)obtained for differentW𝜇and differentQ. The thresholdsQhave been selected in order to compare different sites possessing the same value of𝜎𝜏. The good agree- ment with the theoretical𝜒2distribution is clearly observable. The number of degrees of freedomkcan be related to the exponentqof equation (3) with the relationqth=1+2∕(k+2) (Beck, 2010; Briggs & Beck, 2007; Kosun & Ozdemir, 2016). For MLO atQ=1,q=1.67±0.10, andk=2.797±0.40, using the relation,qth=1.4±0.4was obtained that, within the error bars, is in good agreement with the exper- imental estimation (for CVOQ=1,q=1.690±0.07andk=2.7±0.5givingqth=1.42±0.50; for VRSQ=1.5, q=1.582±0.1,k=2.823±0.5, andqth=1.41±0.5was obtained; and for MHQ=1.5,q=1.61±0.2,k=2.88±0.2, andqth=1.41±0.2).
Figure 10 (right) shows the distributionf(𝜇)for the shuffled data set. In this case the distribution collapses onto the same narrow distribution, completely differently from the other cases, and independent of the
Figure 10.(left) The probability densityf(𝜇)for the normalized IOT sequences obtained for different window sizesW𝜇at a fixed thresholdQ, dashed line𝜒2 distribution withk≈3∘of freedom. (middle) Distributionf(𝜇)for different sites and different thresholdsQobtained for a window length,W𝜇=12days.
The dashed line represents a theoretical𝜒2distribution withk=3∘of freedom. (right) The probability densityf(𝜇)for shuffled data for different sites and different thresholdsQ.
thresholdQ. This profile is in good agreement with the theoretical assumption of a delta function; however, the low data sampling rate and gaps in the data set can affect this shape and introduce a small broadening of the peak.
4. Conclusions
Through an extensive analysis of Hg0concentration data, measured at different latitudes and in different climatological regions, the universal behavior of the IOT extremes has been identified. The Hg0dynamic is characterized by a long-term memory (power law) autocorrelation function, and above a fixed thresholdQ the PDFs of the IOT are described by the Tsallisq-exponential function. The PDF behavior can be explained in terms of the superstatistical nature of the IOTs, where the competition between multiple mesoscopic pro- cesses affects the macroscopic dynamics. Additionally, it seems possible that atmospheric turbulence plays a central role in the dynamics of extreme Hg0concentration events, since the average⟨q⟩=1.59is comparable with the value obtained in fully developed turbulence experiments. The small differences seen in theqvalues obtained from theq-exponential function may be attributable to differences in the evolution of the Hg0con- centration at each site, which would depend on the specific mesoscopic processes that occur locally and their interactions. All processes occurring in atmosphere (e.g., emission, scavenging, chemical reactions, etc.) act as a perturbation on the macroscopic dynamics of the predominant phenomenon, in this case, atmospheric transport. The extensive parameter𝜇encompassing all the possible fluctuations related to the mesoscopic phenomena follows a𝜒2distribution, which is effectively a “fingerprint,” identifying the superstatistical nature of the overall process. By destroying the long-term memory in the data (by shuffling), the PDFs become inde- pendent from the thresholdQand the distributions all collapse on to an exponential distributionP(𝜏∕𝜎Q) = 𝜎−1exp[−𝜏∕𝜎−1]exposing the uncorrelated nature of the shuffled data. However, it is worth mentioning that theqvalues may change slightly depending on the data resolution.
References
Abe, S., & Okamoto, Y. (Eds.) (2001).Nonextensive statistical mechanics and its applications. Berlin Heidelberg: Springer-Verlag.
Alder, B. J., & Wainwright, T. E. (1970). Decay of the velocity autocorrelation function.Physical Review A,1, 18–21.
https://doi.org/10.1103/PhysRevA.1.18
Beck, C. (2010). Generalized statistical mechanics for superstatistical systems.Philosophical Transactions of the Royal Society of London A:
Mathematical, Physical and Engineering Sciences,369(1935), 453–465.
Beck, C., & Cohen, E. (2003). Superstatistics.Physica A: Statistical mechanics and its applications,322, 267–275.
Beck, C., Cohen, E. G. D., & Swinney, H. L. (2005). From time series to superstatistics.Physical Review E,72, 56133.
https://doi.org/10.1103/PhysRevE.72.056133
Bogachev, M. I., & Bunde, A. (2008). Memory effects in the statistics of interoccurrence times between large returns in financial records.
Physical Review E,78, 36114.
Bogachev, M. I., Eichner, J. F., & Bunde, A. (2007). Effect of nonlinear correlations on the statistics of return intervals in multifractal data sets.
Physical Review Letters,99, 240601.
Briggs, K., & Beck, C. (2007). Modelling train delays with q-exponential functions.Physica A: Statistical Mechanics and its Applications,378(2), 498–504.
Bunde, A., Eichner, J. F., Kantelhardt, J. W., & Havlin, S. (2005). Long-term memory: A natural mechanism for the clustering of extreme events and anomalous residual times in climate records.Physical Review Letters,94, 48701.
Carbone, F., & Sorriso-Valvo, L. (2014). Experimental analysis of intermittency in electrohydrodynamic instability.The European Physical Journal E,37, 1–11.
Carbone, F., Gencarelli, C. N., & Hedgecock, I. M. (2016). Lagrangian statistics of mesoscale turbulence in a natural environment: The Agulhas return current.Physical Review E,94, 63101.
Carbone, F., Landis, M. S., Gencarelli, C. N., Naccarato, A., Sprovieri, F., De Simone, F.,…Pirrone, N. (2016). Sea surface temperature variation linked to elemental mercury concentrations measured on Mauna Loa.Geophysical Research Letters,43, 7751–7757.
https://doi.org/10.1002/2016GL069252
Celani, A., Cencini, M., Mazzino, A., & Vergassola, M. (2004). Active and passive fields face to face.New Journal of Physics,6(1), 72.
Chelani, A. (2016). Long-memory property in air pollutant concentrations.Atmospheric Research,171, 1–4.
De Simone, F., Cinnirella, S., Gencarelli, C. N., Yang, X., Hedgecock, I. M., & Pirrone, N. (2015). Model study of global mercury deposition from biomass burning.Environmental Science and Technology,49(11), 6712–6721.
De Simone, F., Artaxo, P., Bencardino, M., Cinnirella, S., Carbone, F.,…Pirrone, N. (2017). Particulate-phase mercury emissions from biomass burning and impact on resulting deposition: A modelling assessment.Atmospheric Chemistry and Physics,17(3), 1881–1899.
Douglas, P., Bergamini, S., & Renzoni, F. (2006). Tunable Tsallis distributions in dissipative optical lattices.Physical Review Letters,96, 110601.
https://doi.org/10.1103/PhysRevLett.96.110601
Dvonch, J. T., Keeler, G. J., & Marsik, F. J. (2005). The influence of meteorological conditions on the wet deposition of mercury in southern Florida.Journal of Applied Meteorology,44(9), 1421–1435.
Eichner, J. F., Kantelhardt, J. W., Bunde, A., & Havlin, S. (2007). Statistics of return intervals in long-term correlated records.Physical Review E, 75, 11128.
Ferri, G., Savio, M. R., & Plastino, A. (2010). Tsallis’ triplet and the ozone layer.Physica A: Statistical Mechanics and its Applications,389(9), 1829–1833.
Ferri, G. L., Figliola, A., & Rosso, O. A. (2012). Tsallis’ statistics in the variability of El Niño/Southern Oscillation during the Holocene epoch.
Physica A: Statistical Mechanics and its Applications,391(5), 2154–2162.
Acknowledgments
F. C., A. G. B., and I. M. H. gratefully acknowledge the two anonymous referees and all the coauthors for all the fruitful discussions and suggestions and for providing the Hg0data used in the analysis.
In particular, A. Steffen from Air Quality Research Division, Environment and Climate Change Canada, Toronto, Ontario M3H 5T4, Canada for the Alert station data. Part of this work was funded through the EU GMOS project (FP7-265113); we also acknowledge the UNEP-GEF (Global Environment Facility) project for funding support to this work.
The EPA through its Office of Research and Development and Office of International and Tribal Affairs partially funded and contributed to this research. The views expressed in this paper are those of the authors and do not necessarily reflect the views or policies of EPA. It has been subjected to Agency review and approved for publication. Mention of trade names or commercial products does not constitute an endorsement or rec- ommendation for use. We thank Ram Vedantham, Sania Tong-Argao, and Carry Croghan (EPA ORD) for their assistance with data Quality Assurance;
and Russell Schell, Aidan Colton, Daryl Kuniyuki, Alan Yoshinaga, and Poai Suganuma (NOAA) for MLO site logistical support. The Hg0data sets can be requested on http://www.
gmos.eu/, with exception for data for Waldhof (WAL) and Mace Head (MH), which are available from the European Monitoring and Evaluation Programme (EMEP) http://www.
nilu.no/projects/ccc/emepdata.html.
The wind speed data for Mauna Loa Observatory can be obtained at https://www.esrl.noaa.gov/gmd/obop/
mlo/summary.html.
Frisch, U. (Ed.) (1995).Turbulence: The legacy of A. N. Kolmogorov. Cambridge UK: Cambridge University Press.
Hayley, H. S., Alexander, H. D. C., Keh-Han, W., Michelle, H. T., & Clark, C. K. L. (Eds.) (2002).Environmental fluid mechanics: Theories and applications. Reston, VA: American Society of Civil Engineer.
Holmes, C. D., Krishnamurthy, N. P., Caffrey, J. M., Landing, W. M., Edgerton, E. S., Knapp, K. R., & Nair, U. S. (2016). Thunderstorms increase mercury wet deposition.Environmental Science & Technology,50(17), 9343–9350.
Horowitz, H. M., Jacob, D. J., Zhang, Y., Dibble, T. S., Slemr, F., Amos, H. M.,…Sunderland, E. M. (2017). A new mechanism for atmospheric mercury redox chemistry: Implications for the global mercury budget.Atmospheric Chemistry and Physics Discussions,17, 6353–6371.
Istas, J., & Lang, G. (1997). Quadratic variations and estimation of the local Hölder index of a Gaussian process.Annales de l’Institut Henri Poincare (B) Probability and Statistics,33(4), 407–436.
Katul, G., Porporato, A., Cava, D., & Siqueira, M. (2006). An analysis of intermittency, scaling, and surface renewal in atmospheric surface layer turbulence.Physica D: Nonlinear Phenomena,215(2), 117–126.
Kosun, C., & Ozdemir, S. (2016). A superstatistical model of vehicular traffic flow.Physica A: Statistical Mechanics and its Applications, 444(Supplement C), 466–475.
Landis, M. S., Stevens, R. K., Schaedlich, F., & Prestbo, E. M. (2002). Development and characterization of an annular denuder methodology for the measurement of divalent inorganic reactive gaseous mercury in ambient air.Environmental Science & Technology,36(13), 3000–3009.
Landis, M. S., Vedantham, R., & Alvarez-Aviles, L. (2013). Measurements of ambient mercury and related species at the Mauna Loa observatory 2002–2010. In11th International Conference of Mercury as a Global Pollutant. Edinburgh, Scotland. Retrieved from https://www.esrl.noaa.gov/gmd/publications/annual_meetings/2010/pdfs/3-Krnavek.pdf
Lindberg, S. E., Brooks, S., Lin, C.-J., Scott, K. J., Landis, M. S., Stevens, R. K.,…Richter, A. (2002). Dynamic oxidation of gaseous mercury in the arctic troposphere at polar sunrise.Environmental Science & Technology,36(6), 1245–1256.
Lovejoy, S., & Schertzer, D. (2013).The weather and climate: Emergent laws and multifractal cascades. Cambridge: Cambridge University Press.
Manshour, P., Anvari, M., Reinke, N., Sahimi, M., & Tabar, M. R. R. (2016). Interoccurrence time statistics in fully-developed turbulence.
Scientific Reports,6, 27452.
Mazzitelli, I., & Lanotte, A. S. (2012). Active and passive scalar intermittent statistics in turbulent atmospheric convection.Physica D:
Nonlinear Phenomena,241(3), 251–259. https://doi.org/10.1016/j.physd.2011.07.009 McComb, W. D. (Ed.) (1990).The physics of fluid turbulence. New York: Oxford University Press.
Monin, A. S., & Yaglom, A. M. (Eds.) (2007).Statistical fluid mechanics: Mechanics of turbulence. New York: Dover.
Pirrone, N., Cinnirella, S., Feng, X., Finkelman, R. B., Friedli, H. R., Leaner, J.,…Telmer, K. (2010). Global mercury emissions to the atmosphere from anthropogenic and natural sources.Atmospheric Chemistry and Physics,10(13), 5951–5964.
Rabassa, P., & Beck, C. (2015). Superstatistical analysis of sea-level fluctuations.Physica A: Statistical Mechanics and its Applications,417, 18–28.
Ryan, S. (1997). The wind field around Mauna Loa derived from surface and balloon observations.Journal of Geophysical Research,102(D9), 10,711–10,725.
Santhanam, M., & Kantz, H. (2005). Long-range correlations and rare events in boundary layer wind fields.Physica A: Statistical Mechanics and its Applications,345(3–4), 713–721. https://doi.org/10.1016/j.physa.2004.07.012
Schertzer, D., & Lovejoy, S. (1985).The dimension and intermittency of atmospheric dynamics(pp. 7–33). Berlin, Heidelberg: Springer Berlin Heidelberg.
Schertzer, D., & Lovejoy, S. (1987). Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes.Journal of Geophysical Research,92(D8), 9693–9714. https://doi.org/10.1029/JD092iD08p09693
Schroeder, W. H., & Munthe, J. (1998). Atmospheric mercury—An overview.Atmospheric Environment,32(5), 809–822.
Schroeder, W. H., Anlauf, K. G., Barrie, L. A., Lu, J. Y., Steffen, A., Schneeberger, D. R., & Berg, T. (1998). Arctic springtime depletion of mercury.
Nature,394(6691), 331–332. https://doi.org/10.1038/28530
Sharma, N. C., & Barnes, J. E. (2016). Boundary layer characteristics over a high altitude station, Mauna Loa Observatory.Aerosol and Air Quality Research,16(3), 729–737.
Sprovieri, F., Pirrone, N., Bencardino, M., D’Amore, F., Carbone, F., Cinnirella, S.,…Norstrom, C. (2016). Atmospheric mercury concentrations observed at ground-based monitoring sites globally distributed in the framework of the GMOS network.Atmospheric Chemistry and Physics,16(18), 11,915–11,935. https://doi.org/10.5194/acp-16-11915-2016
Steffen, A., Douglas, T., Amyot, M., Ariya, P., Aspmo, K., Berg, T.,…Temme, C. (2008). A synthesis of atmospheric mercury depletion event chemistry in the atmosphere and snow.Atmospheric Chemistry and Physics,8(6), 1445–1482. https://doi.org/10.5194/acp-8-1445-2008 Steffen, A., Lehnherr, I., Cole, A., Ariya, P., Dastoor, A., Durnford, D.,…Pilote, M. (2015). Atmospheric mercury in the Canadian Arctic. Part I: A
review of recent field measurements.Science of The Total Environment,509–510, 3–15. Special Issue: Mercury in Canada’s North.
Tsallis, C. (1988). Possible generalization of Boltzmann-Gibbs statistics.Journal of Statistical Physics,52(1), 479–487.
Tsallis, C., & Souza, A. M. C. (2003). Constructing a statistical mechanics for Beck-Cohen superstatistics.Physical Review E,67(2), 26106.
https://doi.org/10.1103/PhysRevE.67.026106
Tsallis, C., Mendes, R., & Plastino, A. (1998). The role of constraints within generalized nonextensive statistics.Physica A: Statistical Mechanics and its Applications,261(3–4), 534–554.
Tuck, A. F. (2010). From molecules to meteorology via turbulent scale invariance.Quarterly Journal of the Royal Meteorological Society, 136(650), 1125–1144.
Van der Straeten, E., & Beck, C. (2009). Superstatistical fluctuations in time series: Applications to share-price dynamics and turbulence.
Physical Review E,80, 36108. https://doi.org/10.1103/PhysRevE.80.036108
van Kampen, N. (Ed.) (1981).Stochastic processes in physics and chemistry. Amsterdam: North-Holland.
Vindel, J. M., & Yagüe, C. (2011). Intermittency of turbulence in the atmospheric boundary layer: Scaling exponents and stratification influence.Boundary-Layer Meteorology,140(1), 73–85.
Warhaft, Z. (2000). Passive scalars in turbulent flows.Annual Review of Fluid Mechanics,32(1), 203–240.
Wilk, G., & Włodarczyk, Z. (2000). Interpretation of the nonextensivity parameterqin some applications of Tsallis statistics and Levy distributions.Physical Review Letters,84, 2770–2773. https://doi.org/10.1103/PhysRevLett. 84.2770
Wyngaard, J. C. (1992). Atmospheric turbulence.Annual Review of Fluid Mechanics,24, 205.
