Atmospheric CO2 and δ13C Measurements from 2012 to 2014 at the Environmental Research Station Schneefernerhaus, Germany: Technical Corrections, Temporal Variations and Trajectory Clustering

Received: January 10, 2018
Revised: July 6, 2018
Accepted: July 10, 2018

Download Citation: RIS|BibTeX|https://doi.org/10.4209/aaqr.2018.01.0010  

Cite this article:

Ghasemifard, H., Yuan, Y., Luepke, M., Schunk, C., Chen, J., Ries, L., Leuchner, M. and Menzel, A. (2019). Atmospheric CO2 and δ13C Measurements from 2012 to 2014 at the Environmental Research Station Schneefernerhaus, Germany: Technical Corrections, Temporal Variations and Trajectory Clustering. Aerosol Air Qual. Res. 19: 657-670. https://doi.org/10.4209/aaqr.2018.01.0010

HIGHLIGHTS

Keywords: CO2 mixing ratio; δ13C; Mountain station; Trajectory HYSPLIT; PSCF.

INTRODUCTION

The sources and sinks of atmospheric carbon dioxide (CO₂) play a critical role in governing global climate. Human activities, especially fossil fuel combustion, cement manufacturing and land use changes trigger anthropogenic emissions and thus contribute to a steady increase of the atmospheric CO₂ mixing ratio (IPCC, 2014; WMO, 2016; Xu et al., 2017). Nearly half of the emitted CO₂ is taken up by both the terrestrial biosphere and absorption in the ocean reservoirs with a similar share each (Le Quéré et al., 2016). The atmosphere and these two major sinks are linked via the balance between photosynthesis, respiration fluxes and CO₂ dissolving in ocean water (Ciais et al., 2013). Considering the risk of the ocean sink to become saturated, the risk of the biosphere sink to turn into a source and of a consequential rapid increase in atmospheric CO₂ growth rates, it is crucial to monitor atmospheric CO₂ mixing ratio continuously with high temporal resolution in order to identify the source and sink contributions (Cramer et al., 2001; IPCC, 2013).

The stable carbon isotope composition of CO₂ (δ¹³C) allows distinguishing the exchange of fossil carbon from the atmosphere and surface reservoir fluxes (Keeling et al., 2011). During photosynthesis, most plants of the terrestrial biosphere prefer to take up the light isotopologue ¹²CO₂ and thus discriminate against ¹³CO₂ (~18‰) (Farquhar et al., 1989). In contrast, during oceanic uptake of CO₂, there is almost no discrimination (~2‰) (Mook et al., 1974). Since the discrimination by oceans is small, the carbon isotope signature (δ¹³C) can be applied to study the contribution of the biosphere on carbon cycles in the atmosphere (Keeling et al., 1989; Miller et al., 2003).

Atmospheric CO₂ and δ¹³C have been studied on both regional and global scales at different sites with various techniques (Levin et al., 1995; Tuzson et al., 2011; Sturm et al., 2013; Moore and Jacobson, 2015; Xia et al., 2015; Pang et al., 2016a). For instance, Xia et al. (2015) analyzed measurements of atmospheric CO₂ concentration and its stable isotope ratios (δ¹³C) at the regional background station Lin’an (LAN) in China to identify the isotopic signature of CO₂ sinks and sources. They concluded that during the winter season (Dec–Feb) coal combustion is the major CO₂ source due to domestic heating. The relatively high isotopic signature (~21.32‰) of sources and sinks during the vegetation season (Mar–Nov) was attributed to the significant contribution of biological activities at LAN. Sturm et al. (2013) presented continuous measurements of atmospheric CO₂ isotopes (δ¹³C and δ¹⁸O) at the High Altitude Research Station Jungfraujoch (JFJ), Switzerland. Based on the high temporal resolution of their measurements (e.g., in comparison to flask samples), diurnal and hourly variations could be analyzed. They determined that diurnal cycles of atmospheric CO₂ and its isotopic compositions were small; however, the day-to-day variability depending on the origin of the air masses arriving at JFJ was larger. Based on this, the footprint clustering by the backward Lagrangian particle dispersion model FLEXPART revealed different CO₂, δ¹³C, and δ¹⁸O values depending on the origin of air masses and surface residence time. Moore and Jacobson (2015) differentiated sources of atmospheric CO₂ by records of CO₂ concentration and carbon isotope composition (δ¹³C) at Evanston, an urban site north of Chicago, the third largest city in the United States.

Measurements of CO₂ mixing ratios in combination with stable CO₂ isotopes in the atmosphere at high altitude (background) sites may provide valuable information on carbon source and sink mechanisms (Zhou et al., 2006; Sturm et al., 2013; Yuan et al., 2018). The variability of the atmospheric CO₂ mixing ratio and δ¹³C is related to the history of air masses arriving at a measurement site, such as whether the air traveled through the free troposphere or has been influenced by the planetary boundary layer (PBL). Therefore, many atmospheric studies have applied trajectory statistical methods to identify the sources and their contribution to mixing ratios. Fleming et al. (2012) reviewed more than 150 studies dealing with back trajectories, cluster analyses and residence time. However, most of them applied air-mass history to identify transport pathways and potential sources but not sink contributions. It is well known that the effect of biospherical activity as well as strong influences of thermal convection (local transport processes from the boundary layer) may complicate the study of atmospheric CO₂ and δ¹³C (Zellweger et al., 2003; Sturm et al., 2013).

In the present study, we (1) describe the specific corrections and selections of the data recorded by a non-upgraded G1101-i analyzer, (2) analyze the diurnal and seasonal patterns of atmospheric CO₂ mixing ratios and δ¹³C and (3) identify potential CO₂ source and sink locations with the help of backward trajectories and the potential source contribution function.

MATERIAL AND METHODS

Measurement Site

The Environmental Research Station Schneefernerhaus (UFS) is part of the Global Atmosphere Watch (GAW) Global Observatory Zugspitze/Hohenpeissenberg in Germany. It is located in the northern limestone Alps at the border of Germany and Austria (47°25ʹ00ʺN, 10°58ʹ46ʺE) about 90 km southwest of Munich. The UFS is situated at an altitude of 2650 m above sea level (a.s.l.) on the southern slope of Zugspitze mountain massif. It can receive representative free tropospheric air but is also influenced by polluted air masses from the planetary boundary layer (PBL), especially during the day and in summer (Zellweger et al., 2003; Pandey Deolal et al., 2014; Leuchner et al., 2016). Due to its location on a southern slope, northerly winds are mostly blocked (see Fig. 1) while there is wind channeling of western and eastern winds by the west side of the mountain ridge and the Rein Valley, respectively (Gantner et al., 2003; Risius et al., 2015). Due to diverse and intensive human activities in the vicinity of the UFS, from skiing and hiking areas as well as associated tourist services, the GAW site is categorized as a “weakly influenced, constant deposition” site (Henne et al., 2010; Ferrarese et al., 2015). Further detailed information can be found in the UFS station information system (http://www.schneefernerhaus. de/en/home.html).

Fig. 1. Windrose at UFS, obtained from hourly wind data over the entire study period (May 2012–November 2014). Wind speed and direction are classified into colored segments; gray circles show the cumulated percentage of occurrence.

Instrumentation

A wavelength-scanned cavity ring down spectrometer (G1101-i, Picarro Inc. USA) was installed in the laboratory on the 4^th floor of the UFS building. It was operated from May 2012 to November 2014 (except for a downtime period in May and June 2013 due to a pump defect and its subsequent replacement). CO₂ mixing ratios and δ¹³C were measured at a rate of 0.1 Hz. The air inlet for the instrument was located on a terrace above the laboratory roof at an altitude of 2670 m a.s.l. The inlet cap was constantly heated to prevent ice formation. The inner part of the inlet tube was made of borosilicate glass and was constantly regulated to a low positive temperature of ca. 5°C to avoid water condensation. The inlet tube (length: 3.5 m) was connected to a common manifold (glass, 4.2 m in length, 8 cm inner diameter) that was used for all measuring devices and species at the GAW station (Zellweger et al., 2011). A 2 m 1/8ʺ stainless steel pipe connected the instrument to the inlet line, via a VICI (Valco Instruments Company Inc.) rotary valve that switched between ambient air and standard gases.

Isotopic ratio measurements of ¹³C/¹²C are expressed in per mil (‰), defined by the following Eq. (1) as:

Isotopic values were given relative to the international standard VPDB-CO₂ (Vienna Pee Dee Belemnite) (Brand et al., 2010). Measurements of CO₂ and δ¹³C were calibrated using two working standard gases with high and low concentrations (DEUSTE Steininger GmbH) of the certified CO₂ mixing ratios and the isotopic compositions (Standard 1: 350.1 ± 0.50 ppm CO₂ mixing ratio and –3.28 ± 0.164‰ δ¹³C; Standard 2: 503.4 ± 0.50 ppm CO₂ mixing ratio and –20.03 ± 1.002‰ δ¹³C) in synthetic air. Each standard gas was fed into the analyzer for one hour every 7 days using an open-split configuration. After the data collection period, both gas standards were reanalyzed at the laboratory of the Max Planck Institute for Biogeochemistry (MPI-BGC), Jena, Germany, in order to confirm the quality of the calibration and the stability of the standards. The measured δ¹³C values (Standard 1: –3.19 ± 0.009‰ δ¹³C; Standard 2: –20.26 ± 0.01‰ δ¹³C) agreed well with the original values within the given range of GAW data quality objectives (±0.01‰). The standard gases were also analyzed for other trace gas concentrations. Except for N₂O (Standard 1 and 2: 0.1 ppm) and CO (Standard 1: 6 ppb; Standard 2: 5.47 ppb), no other species (CH₄, H₂O, SF₆) could be detected.

The German Meteorological Service (Deutscher Wetterdienst, DWD) provided wind speed and wind direction data from the UFS in a one-minute time interval, and for this study, hourly averaged data are presented. In addition, we used CO₂ and CH₄ data measured with a Picarro Envirosense 3000i instrument to correct our measured data (see the section “Correction of CO₂”). This instrument, using cavity ring down spectroscopy, was operated by the German Environment Agency (Umweltbundesamt, UBA) connected to the same inlet as our Picarro G1101-i, therefore measured identical air sample, and was operated in the same air conditioned laboratory. For this device, water vapor in the ambient air was removed by cold traps. Gas mixtures for working standards of this device were equally delivered by DEUSTE Steininger GmbH, Germany. The station standards for interconnection with the international standard reference were reported on WMO X-2007 scale for CO₂ and WMO X2004a scale for CH₄ by NOAA, Boulder, Colorado, USA. Calibration and quality assurance for these atmospheric compounds followed the standard operating procedures of UBA in accordance with GAW quality standards. Since CO₂ measurements from UBA are reported to the WMO World Data Centre for Greenhouse Gases, the respective CO₂ data in the section “Correction of CO₂” will be referred to as CO₂ data of Global Atmosphere Watch (CO_{2, GAW}).

Optimal Integration Time and Precision

The Allan variance method (Allan, 1966; Werle et al., 1993; Chen et al., 2016) was applied to determine the optimal integration time for the measurements and to determine the best precision, using Standard 1 (CO₂: 350.1 ± 0.50 ppm; δ¹³C: –3.28 ± 0.164‰; Fig. 2). For this purpose, a long-term measurement of 24 h duration in February 2014 was analyzed. When integrating less than the optimum integration time, the Allan deviation follows a slope of –1/2 in the double logarithmic scale, indicating that white noise is dominating. When integrating beyond the optimum integration time, the Allan deviation rises and follows a slope of 1/2, indicating an instrument drift.

Fig. 2. Allan deviation plots for CO₂ (left) and δ¹³C (right) as a function of the integrating time τ, based on a long-term measurement of Standard 1 (CO₂: 350.1 ± 0.50 ppm; δ¹³C: 3.28 ± 0.164‰). f stands for frequency and the black dashed lines represent slopes of –1/2 and 1/2, which correspond to power spectral densities S(f) = f⁰ and S(f)= f², respectively. The Allan deviation follows a slope of –1/2 up to an integration time of 1 hour for CO₂ and of 2 hours for δ¹³C and then turns over to a slope of 1/2 which defines a drift.

The optimum integration time minimizing the Allan deviation was around 1 hour for CO₂ and around 2 hours for δ¹³C. The best achievable precision (1 sigma) was 4 ppb for CO₂ and 0.2‰ for δ¹³C, respectively (see Fig. 2). However, we used a shorter integration time of 30 minutes for the ambient measurements, which gave a precision of 0.4‰ for δ¹³C and 5–6 ppb for CO₂. For the very same analyzer, Wen et al. (2013) found the best precision of 0.08‰ at 2000 s, Vogel et al. (2013) had a precision of 0.2‰ at 5 min averaging, and Pang et al. (2016b) achieved optimum values of 0.08, 0.15, and 0.10‰ at 7600, 1900, and 1900 s for three reference gases. The precision of all studies on Picarro G1101-i shows better precision than the specification provided by the manufacturer (0.3‰).

Data Calibration and Correction

The G1101-i analyzer in this study was manufactured in 2010 and not upgraded during the measurement period to minimize downtime. As an alternative to upgrading, Picarro Inc. recommended later removal of spectral interference caused by CH₄ that can bias δ¹³C by 0.4‰ ppm^–1 (Vogel et al., 2013) as well as water interferences. The latter include water vapor dilution, water vapor pressure broadening, and HDO spectral interference effects (Wen et al., 2013). As we did not use any drying system and therefore measured the humid ambient gas, corrections for the dominating water vapor dilution effect were done according to Hoffnagle (2013); however, the smaller effects due to the water vapor pressure broadening and the HDO spectral interference effects were not corrected and their effects are about 2 ppm %v^–1 water at 400 ppm of carbon dioxide, and up to 5‰ at ambient humidity (Nara et al., 2012; Rella et al., 2013; Wen et al., 2013). According to Wen et al. (2013), even an upgrade of the analyzer would have led to overcorrections. The dependency of δ¹³C on the CO₂ mixing ratio is of no concern to this study, as Vogel et al. (2013) showed that no dependency could be detected in the range of 303–437 ppm. Background gas concentrations and ratios (e.g., N₂/O₂ ratio and argon content) are known to influence the CRDS technology in general. While their natural variations usually generate negligible effects, large differences between ambient air and synthetic standards, such as the ones used in our study, may be problematic (Nara et al., 2012; Rella et al., 2013). However, no specific information or corrections are documented for the G1101-i analyzer.

Correction of δ¹³C

The isotopic analyzer reports the peak height of the near-infrared absorption spectrum for the rovibronic transition of ¹²C¹⁶O₂ and ¹³C¹⁶O₂ with arbitrary labeled C12_peak-BookAve and C13_peak-BookAve, respectively. The ratio of the peak absorption values and the ratio of isotopic abundances are linear.

A and B are instrument-specific constants and were provided by the manufacturer; they are the linear slope and intercept terms for computing the delta value from the ratio between C13_peak-BookAve and C12_peak-BookAve (Hoffnagle, 2013). The peak height can be affected by the absorption peaks for water, methane and other gases (Nara et al., 2012; Rella et al., 2013).

δ¹³C ratios were corrected for methane and water vapor according to Hoffnagle (2013) using Eq. (3):

H₂O is the water vapor concentration (in percent) and
(–2.98648648 × H₂O) is the correction for the water interference. CH₄ is the methane concentration (in ppm), measured by the UBA device (Picarro Envirosense 3000i instrument).

Correction of CO₂

To quantify and correct water vapor effects on the CO₂ measurements, we compared the CO₂ mixing ratio for the wet gas stream measured by Picarro G1101-i analyzer (CO_2,meas) with CO₂ data of the accompanying instrument (CO_2,GAW). Although similar seasonal cycles were observed in Fig. 3, CO_{2, meas} mixing ratios only showed a good agreement with CO_{2, GAW} in wintertime when water vapor concentrations were low (0.39 ± 0.004%), while clear differences were observed in the summertime when water vapor concentrations were higher (1.81 ± 0.0048%). Wen et al. (2013) also measured a wet gas stream and compared the same analyzer as in this study with the Los Gatos DLT-100, observing that the mixing ratio measured by Picarro G1101-i was 2.2 ± 1.0 ppm lower than Los Gatos. They suggested that the water vapor dilution effect was partly responsible for this difference.

Fig. 3. Comparison of CO₂ mixing ratio measured by Picarro G1101-i (in red) without drying system to CO₂ measured by Picarro Envirosense 3000i (in black) dried by cold traps.

We corrected the dilution effect using ordinary least squares regression of the CO₂ mixing ratio difference between CO_2,GAW and CO_{2, meas}. Concurrent with a pump failure and subsequent replacement (end of May/beginning of June 2013), the extent of the dilution effect changed and thus OLS regressions were carried out separately for the period before (Eq. (4)) and after (Eq. (5)) the pump failure, yielding

[CO₂]_calculated = 6.5[H₂O] – 2.7 + [CO₂]_meas                  (4)

[CO₂]_calculated = 5.6[H₂O] – 0.85 + [CO₂]_meas                 (5)

where [H₂O] is water vapor concentration in percent, [CO₂]_meas is CO₂ mixing ratio measured by our analyzer and [CO₂]_calculated is the calculated and corrected mixing ratio. Fig. 4 illustrates the regressions before and after the pump failure. The physical reasons for the change in G1101-i behavior are not known and it seems unlikely that they were caused by the pump failure and replacement alone.

Fig. 4. Scatterplot of water vapor concentrations and difference of GAW and Picarro CO₂ measurements in the two phases before (a) and after (b) the pump replacement. The dashed line is an ordinary least squares regression line (coefficient of determination of linear regressions in (a): R² = 0.97, p-value < 0.001 and (b): R² = 0.92, p-value < 0.001).

Corrected CO₂ data thus is in accordance with GAW data quality objectives (±0.1 ppm) (WMO, 2016). However, very small differences between the two analyzers (0.01 ± 0.42 and 0.02 ± 0.06 ppm before and after the pump replacement, respectively) remain and may be due to biases from spectral broadening and interferences (Rella, 2012). 

Data Coverage and Selection

The temporal resolution of CRDS raw data should be in the order of seconds. In our case, however, the data periodically included measuring intervals far beyond the normal range (55 s and more). Since long measuring intervals are indicative of a slowed ring-down frequency and potential problems with the instrument (such as laser ageing or deficiencies in the optics, detector or data acquisition system), intervals > 55 s duration were removed from the dataset (personal communication by Dr. Renato Winkler, Picarro Inc.), which amounted to about 4.6% of data. The device returned to normal measuring intervals after restarts or laser readjustments. 1.8% of the data were missing due to the pump defect and replacement in May and June 2013. Another 3.8% of the data were missing due to the calibration procedures and/or power failures in the lab. At the end, this resulted in 89.8% of valid data over the entire measuring period.

The first 29 measurements of each calibration (ca. 3.5 minutes, less than 6% of each calibration) were discarded to ensure that no air sample from the previous measurement was left due to the transient response after valve switching (Vogel et al., 2013). Then, from the remaining data, the average of each calibration was calculated. A smoothing spline was then fitted to all the averages to account for the residual variation in the calibration data and to fill a gap in the calibration measurements (8 months, July 2012–March 2013). This reduced the residual deviation between calibrated standard measurement and the real standard value from 0.05 ± 0.75‰ and 0.03 ± 0.53‰ (linear interpolation) to 0.0 ± 0.44‰ and 0.0 ± 0.32‰ (smoothing spline) for Standards 1 and 2, respectively. Splines were calculated using the sm.spline function from the pspline package in R (Heckman and Ramsay, 1996). This function fitted a natural polynomial smoothing spline to the calibration values with 25 degree of freedom (out of 108 calibrations). Deviations between spline interpolated calibrations and the real standard values for each calibration are shown in Fig. 5. The smoothed coefficients were used for calibration. The two-point mixing ratio gain and offset calibration strategy of Bowling et al. (2003) was used for each measurement cycle as described above. Data were aggregated to 30-min averages by the statistical program R (R Core Team, 2016), which together with several packages was used for further analyses. At the end, uncertainty of calibration is estimated to be ±0.56 ppm for CO₂ and ±0.53‰ for δ¹³C and combined measurement uncertainty is estimated to be ±0.56 ppm for CO₂ and ±0.56‰ for δ¹³C.

Fig. 5. Histogram of the deviation between spline and true standard value for (a) Standard 1 and (b) Standard 2.

HYSPLIT Trajectory Model

One common method to establish source/sink-receptor relationships is to combine a calculated trajectory path of an air parcel with measured data at the time when the air parcel arrives at the site and consequently to determine locations of sources and sinks from these observations (Stohl, 1996).

In this work, the HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) Model was used to calculate backward trajectories to estimate air mass pathways to the UFS. Trajectories were calculated hourly for 96 h backward ending at UFS for the whole measurement period. The backward trajectory calculation was started at an altitude of 1500 m above ground level (a.g.l.) with respect to the model elevation at the coordinates of UFS, thus 3000 m a.s.l., roughly matching the real site altitude of 2650 m a.s.l. (UFS) and 2670 m a.s.l. (sample inlet).

The GFS (Global Forecast System) model “Grid 4” forecast weather data with 0.5° resolution was used as an underlying meteorological model. Forecasts from the 18:00 UTC cycle were obtained from the NOAA HAS data portal each day for the next 24 hours and were converted to ARL input model format. In case of the very few gaps (around 1.2%) in the GFS forecast data, forecasts from the previous day(s) were used as a replacement. For each trajectory time step, coordinates, altitude and mixing layer height (from the input model) were extracted. The backward trajectories were further processed by the openair package (Carslaw and Ropkins, 2012) within R. Here, clusters were calculated from all trajectories by applying an angle-based distance matrix method with the k-medoids algorithm (Sirois and Bottenheim, 1995). The number of clusters was set to five since for this number at least 10% of the trajectories were represented within each cluster. In order to filter the trajectories concerning their background air characteristics, they were split into two classes: (1) FT (free troposphere) trajectories when the trajectory height was higher than the mixing layer height for all time steps and (2) PBL influenced trajectories when the trajectory was below the calculated mixing layer height for at least one time step.

The potential source contribution function (PSCF) method using the residence time probability (Ashbaugh et al., 1985; Seinfeld and Pandis, 2016) describes the spatial distribution of probable geographical source locations derived by trajectories. PSCF is defined as the probability that an air parcel with atmospheric component concentrations higher than a specific threshold arrives at the receptor site after having been observed to reside in a certain grid cell. The PSCF value (P_ij) for a given grid cell is then calculated as P_ij = m_ij/n_ij, in which n_ij is the total number of trajectory segment endpoints terminating within the grid cell(i,j) over the entire time of measurement and m_ij is the number of trajectory segment endpoints terminating within the grid cell(i,j) corresponding to trajectories associated with concentration values at the receptor site greater than a specific threshold (75^th percentile of CO₂ and δ¹³C in this study).

Robust Extraction of Baseline Signal

In order to calculate the long-term background values, the statistical method REBS (Robust Extraction of Baseline Signal) was applied, which is based on robust local regression (Ruckstuhl et al., 2012). REBS is very flexible to derive the background levels of various trace gases at background measurement sites simultaneously (in our case both CO₂ and δ¹³C) due to its non-parametric basis. Sturm et al. (2013) used REBS for high alpine CO₂ and δ¹³C measurements at Jungfraujoch, with results very similar to those obtained from the commonly applied data filtering method by Thoning et al. (1989). The background concentration was extracted using local regression (60 day windows) implemented in the rfbaseline function of the IDPmisc package (Locher and Ruckstuhl, 2012) in R. Missing data points were interpolated by a simple linear interpolation using the values from the hours before and after. The mean peak-to-peak seasonal amplitude was calculated from the background curves for an entire 2-year period (2012–2014). 

Remote Sensing Data

In order to assess the vegetation properties, information on the vegetation type and activity in summer (June–August) were required. Therefore, we used remotely sensed phenological data based on MODIS Land Cover (LC) data and the Normalized Difference Vegetation Index (NDVI), respectively. The latest LC was issued in 2012 at 0.5° × 0.5° resolution (Friedl et al., 2010; Channan et al., 2014). More information can be found in the Global Land Cover Facility (http://glcf.umd.edu/data/lc/). The LC classes were aggregated to the ones shown and described in Fig. 11(b). NDVI is derived from multiple Advanced Very High Resolution Radiometer (AVHRR) measurements and is an indicator of the greenness of vegetation in each pixel of the satellite image. It ranges from –1.0 to +1.0. NDVI values below zero are excluded from this study because they indicate no vegetation, such as rock, sand or snow. Sparse vegetation such as grassland or cropland results in moderate NDVI values (0.2–0.5) and dense vegetation such as tropical forests results in high NDVI values (0.6–0.9). In this study, we averaged the respective NDVI values for the three months of summer (June to August) in the year 2012. The NDVI data were downloaded from the ECOCAST directory of the NASA Ames Ecological Forecasting Lab, version 3g.v0 (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v0/) for 2012.

RESULTS AND DISCUSSION

Atmospheric CO₂ and δ¹³C

Time series of CO₂ mixing ratios and δ¹³C recorded from May 1, 2012, to Nov. 2, 2014, are displayed in Fig. 6. This includes both hourly data and fitted background concentrations, using the REBS technique. Fig. 6 shows a clear seasonal variation of hourly mean CO₂ mixing ratios and δ¹³C values. CO₂ mixing ratios simultaneously increase as δ¹³C values decrease, which is associated with seasonal vegetation activity. The mean peak-to-peak amplitudes are 15.5 ± 0.15 ppm for CO₂ and 1.97 ± 0.53‰ for δ¹³C. The minimum CO₂ mixing ratios, as well as the maximum δ¹³C values, occurred in August due to (preferentially ¹²CO₂) CO₂ terrestrial uptake being dominated by the biosphere. The maximum in CO₂ mixing ratios occurred in March and the minimum in δ¹³C values in February when respiration is dominating. In contrast, the peak-to-peak amplitudes of CO₂ mixing ratios and δ¹³C at the High Altitude Research Station Jungfraujoch (JFJ) on the northern ridge of the Swiss Alps (46°32ʹ53ʺN, 7°59ʹ2ʺE, 3580 m a.s.l.) are 11.0 ppm for CO₂ and 0.60‰ for δ¹³C (Sturm et al., 2013). These differences in seasonal amplitudes between UFS and JFJ are most likely due to altitude since JFJ is located 930 m higher than UFS. Consequently, air masses at JFJ are much less impacted by PBL and mostly are from the lower free troposphere. Although amplitudes of CO₂ mixing ratios and δ¹³C from these two sites were different, the minimum and the maximum of CO₂ mixing ratios and δ¹³C were recorded nearly at the same time of the year.

Fig. 6. Time series of hourly mean values of CO₂ and δ¹³C for the period May 2012 to November 2014. Gray dots and solid lines represent the measurements and fitted background curves, respectively.

Diurnal Cycles

The mean diurnal cycles of CO₂ and δ¹³C for each season are shown in Fig. 7. The CO₂ level started to increase at about 06:00 and reached the maxima around 09:00 to 10:00 during spring, summer, and autumn and at around 13:00 during winter, possibly due to local tourist activity, regional respiration, and regional anthropogenic emission potentially due to traffic. The daily peak-to-peak amplitude is 1.4 ppm in winter and 1.6 ppm in spring and autumn. The maximal diurnal change in CO₂ is found during summer months with a peak-to-peak amplitude of about 2.9 ppm. The strong afternoon drop during summer months is due to upward transport of PBL air from the valley. In summer, the CO₂ of the PBL air in lower levels is depleted due to photosynthetic uptake, indicating the influence of vegetation activity in the lower local and regional area around (below) the UFS. During summertime, the PBL influence at UFS has been identified for other parameters as well. Namely, the diurnal variation and high standard deviation of formaldehyde (HCHO) mixing ratio (Leuchner et al., 2016) and the dependency of aerosol concentrations from the altitude of the mixing layer (Birmili et al., 2009) could be detected. The mean diurnal cycle of δ¹³C is very pronounced in summer with an amplitude of 0.4‰. The respective diurnal amplitudes in spring and autumn are 0.2‰ each, showing similar patterns as in summer. In contrast, winter months do not display any distinct diurnal cycles of δ¹³C.

Fig. 7. Diurnal variations of the mean hourly CO₂ mixing ratios and δ¹³C (solid line) relative to the respective daily means for the different seasons. Dashed lines show the 95% confidence intervals of the hourly mean calculated by bootstrap re-sampling.

For the Swiss site Jungfraujoch (JFJ), similar start of increase and maxima have been observed in spring, summer and autumn. However, amplitudes of diurnal variations in CO₂ mixing ratios and δ¹³C were different (Sturm et al., 2013). The diurnal peak-to-peak amplitude of CO₂ at JFJ in summer was 2 ppm (0.9 ppm less than at the UFS) and was 1 ppm in the other seasons (0.6 ppm less than at UFS in spring and autumn and 0.4 less in winter). In contrast to UFS, there were almost no diurnal variations of δ¹³C at JFJ in spring, autumn and winter. In summer, JFJ exhibited a diurnal variation of about 0.3‰ less than at UFS with a comparable minimum in the morning and maximum in the afternoon. The larger amplitudes at UFS are likely due to a lower elevation (i.e., the station being closer to sinks and sources in the valley) and more tourist activity in the vicinity of UFS. Compared to urban sites, which are influenced much more strongly by local biogenic and anthropogenic activities, UFS shows much lower diurnal variability, e.g., by factors of 10 and 34 for δ¹³C and CO₂, respectively, for a site in Kraków, Poland (Zimnoch et al., 2004).

Back Trajectory Clusters

Characterization of the Clusters

The back trajectory analysis described in the section “HYSPLIT trajectory model” is shown in Fig. 8. Cluster 3 (C3) comprises the fastest and the most frequent (41.7%) flow at UFS and corresponds to a westerly flow from the mid-Atlantic Ocean. Cluster 4 is the second most abundant (19.5%), comprising flow from the northwest. Cluster 2 (14.9%) represents a moderately fast southwesterly flow. Slow-moving air masses coming from the central and west Mediterranean Basin are grouped in Cluster 1, representing 12.1% of the data. Cluster 5 corresponds to northeasterly flows and accounts for the smallest share (11.9%) of trajectories.

Fig. 8. Clustering of trajectories arriving at UFS in the entire period of May 2012 to November 2014 with their percentages of trajectories.

Variation of CO₂ Mixing Ratios and δ¹³C among the Clusters

In order to study whether measured data of CO₂ and its isotopic composition systematically varied with air mass origin, we merged hourly CO₂ mixing ratios and δ¹³C data with the respective trajectory information of the same time step (hourly) and grouped them into the five clusters described above (Fig. 8). Afterwards, data were filtered in two steps in order to separate the short-term deviations from the background concentrations and to derive clear information on potential background air masses. The background concentrations derived from the local regression were first subtracted from CO₂ and δ¹³C values, providing ΔCO₂ and Δδ¹³C. Based on the 96-hour back trajectories of the HYSPLIT Model (see the section “HYSPLIT trajectory model”), the latter data were grouped into air masses in contact with PBL or of completely free troposphere origin. ΔCO₂ and Δδ¹³C by trajectory clusters and seasons are shown in Fig. 9. The highest frequency of air masses contacting the PBL occurred in summer, with 21.3%, while the lowest frequency occurred in winter, with 10.3%. Spring and autumn had contact with the PBL in 14.6% and 12.6% of the cases, respectively.

Fig. 9. Boxplots of ΔCO₂ mixing ratio and Δδ¹³C for all seasons classified for the five clusters of Fig. 8 for the entire period. ΔCO₂ and Δδ¹³C are filtered for free troposphere air and background concentrations are subtracted.

In summer, CO₂ mixing ratios of clusters from the south significantly differed from the respective median of clusters from the north and west (significant levels of the Student’s t-test are p < 0.001). In wintertime, the two clusters from the north (C4 and C5) showed the most pronounced differences in the measured ΔCO₂ and Δδ¹³C. Clusters 4 and 5 had the highest median ΔCO₂ values and, correspondingly, the lowest values in Δδ¹³C.

Using footprint clustering analysis for ΔCO₂ and Δδ¹³C in wintertime for the JFJ site, Sturm et al. (2013) revealed that a cluster representing air masses with surface contact mostly over northern European land masses had the highest median ΔCO₂ value and the second lowest value in Δδ¹³C. Another cluster with high CO₂ mixing ratio had its origin in Eastern Europe. These two clusters according to residence time maps had the same direction as our C4 and C5 clusters, respectively. C4 was always associated with the lowest Δδ¹³C. This can be due to anthropogenic emissions from this direction (as shown in Fig. 10) being prevalent in all seasons. Cluster 3 had the least differences across seasons, both for ΔCO₂ and Δδ¹³C. In summer, C5 exhibited the highest Δδ¹³C value and the lowest ΔCO₂ value, respectively.

Fig. 10. Potential source contribution function plot for ΔCO₂ in wintertime (December to February). (a) PSCF map of Cluster 4 (NW) and (b) PSCF map of Cluster 5 (NE). The position of the site is shown by black circles.

Case Study: Combination of Trajectory Clusters and PSCF

The potential source contribution function (see the section “HYSPLIT trajectory model”) was applied for the detection of geographical areas with an influence on the measured CO₂ concentration. This was done explicitly for air masses causing extremely high CO₂ concentrations (> 75^th percentile) in Clusters 4 and 5. PSCF maps (Fig. 10) depict that the air masses from C5 can be traced back to the anthropogenic emissions from coal mining areas in Lusatia and from coal districts in East Germany, while air masses from C4 originate from northwestern Europe, including the high emission regions in the Netherlands and the German Ruhr area.

A comparable PSCF map for Δδ¹³C (> 75^th percentile), equally inferred by trajectories, describes the spatial distribution of probable geographical source locations for δ¹³C (i.e., sinks of CO₂) in summer (Fig. 11). The combination of this PSCF map with MODIS land cover (LC) and NDVI map improve our understanding of the influences of air origin on seasonal variations in CO₂ andδ¹³C. Fig. 11 clearly illustrates that high values of Δδ¹³C are influenced by Western, Central, and Southwestern Poland, which according to LC are predominantly croplands. Fig. 11(c) displays that this area is characterized by NDVI values greater than 0.6, indicating dense vegetation in these area with CO₂ uptake.

Fig. 11. (a) Potential source contribution function plot for Δδ¹³C in the summertime (June to August) and (b) MODIS Land Cover for the year 2012. Land cover types on this map are grouped and labeled as follow: 1: water; 2: evergreen, deciduous as well as mixed forests; 3: closed and open shrublands, woody savannas, savannas, grassland; 4: croplands; 5: cropland/natural vegetation mosaic and 6: urban and build-up, (c) NDVI map. The position of the UFS site is shown by white circles.

CONCLUSIONS

We presented a 2.5-year measurement time series of CO₂ mixing ratios and δ¹³C at the high altitude GAW station UFS on the northern Alpine ridge. Since the Picarro G1101-i instrument was not upgraded to account for water vapor and methane interferences, comprehensive external data corrections and selections had to be implemented. These corrections would not have been feasible had there not been access to parallel measured carbon dioxide, methane and water vapor. Therefore, from the knowledge obtained with CO₂ measurement devices, adding a reliable drying system on the sample inlet line is strongly recommended.

patterns as at the High Altitude Research Station Jungfraujoch. However, the higher altitude of JFJ led to receiving more frequent free troposphere air masses, and therefore, amplitudes at UFS were generally larger than at JFJ in all seasons. The most pronounced diurnal variabilities were seen in summer at both sites, with an amplitude of 2.9 ppm and 2 ppm for the CO₂ mixing ratio at UFS and JFJ, respectively, and of 0.4‰ and 0.1‰ for the δ¹³C at UFS and JFJ, respectively. The smallest diurnal variabilities occurred in winter, when the peak-to-peak amplitude was 1.4 ppm and 1 ppm for the CO₂ at UFS and JFJ and there was no discernible cycle for δ¹³C.

HYSPLIT classification of the air mass origins at UFS indicated predominant air masses from the west (41.2%), followed by the northwest (19.7%), southwest (14.8%), southeast (12.5%) and northeast (11.8%). The potential source contribution function, using back trajectories as well as atmospheric measurements, provided helpful indications of the origins of air masses potentially influencing a measuring station. So far, this method has been restricted to the identification of sources in air pollution studies (Begum et al., 2005; Pekney et al., 2006; Kaiser et al., 2007; Pongkiatkul and Oanh, 2007; Zhu et al., 2011), as the probability function in PSCF filters out concentrations below a(n) (arbitrary) threshold. However, in this study, δ¹³C measurements and their inverse relationship to CO₂ allowed us to identify atmospheric CO₂ sinks. Among all clusters, high CO₂ mixing ratios in winter (anthropogenic sources) were associated with air masses originating in the Netherlands, the German Ruhr area and Lusatia, whereas high δ¹³C values and low CO₂ mixing ratios in summer, representing a terrestrial biosphere sink influence, mostly originated in Poland. For an improved understanding of the contribution of the sources and sinks of carbon dioxide and the exchange of CO₂ between the terrestrial ecosystem and atmosphere, greenhouse gas emission models (e.g., the WRF Greenhouse Gas Model) can be used to simulate high-resolution transport of carbon dioxide.

Sinks and sources’ isotopic signatures of continuously measured CO₂ (such as the data presented in this publication) or records of distinct CO₂ variations (e.g., a sudden enhancement in atmospheric CO₂) could be further identified and quantified using the Keeling plot method (Keeling et al., 1989; Vardag et al., 2016). Access to other related tracers, e.g., carbon monoxide and radon-222 measured at the same site (as is done at UFS), may be a great help in this venture (Hirsch, 2007; Tuzson et al., 2011).

ACKNOWLEDGMENTS

This study was funded by the Virtual Alpine Observatory project of the Bavarian State Ministry of the Environment and Consumer Protection. The authors would like to acknowledge Dr. Steven C. Wofsy and Maryann R. Sargent for fruitful discussions in data selection and correction, Rachel Chang, for her suggestions regarding source detections, and UFS staff for their kind support especially Dr. Rehm and Dr. Couret for essential help, including maintenances and technical supports, and UBA for the space in their lab and access to the facilities. The simulations of back trajectories were performed by Dr. Stephan Hachinger on the Compute Cloud of the Leibniz Supercomputing Centre (LRZ), Garching, Germany, and thanks for his important and valuable contribution to this research.