Homa Ghasemifard 1, Ye Yuan1, Marvin Luepke1, Christian Schunk1, Jia Chen2,3, Ludwig Ries4, Michael Leuchner1, Annette Menzel1,3

Professorship of Ecoclimatology, Technische Universität München, 85354 Freising, Germany
Professorship of Environmental Sensing and Modeling, Technische Universität München, 80333 Munich, Germany
Institute for Advanced Study, Technische Universität München, 85748 Garching, Germany
German Environment Agency (UBA), 82475 Zugspitze, Germany

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

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

  • Download: PDF

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


  • 2.5-year dataset of continuous measurements of atmospheric CO2 and δ13C.
  • Diurnal variations of CO2 and δ13C in different seasons were influenced by PBL.
  • The arriving air masses at the site were separated into 5 clusters by HYSPLIT model.
  • Influences of air masses origins on CO2 mixing ratio and δ13C were investigated.
  • PSCF method was applied to identify source and sink locations.


This study presents continuous atmospheric CO2 and δ13C measurements by wavelength-scanned cavity ring down spectroscopy (Picarro G1101-i) at the high-mountain station Schneefernerhaus, Germany. δ13C values were post-corrected for methane and water spectral interferences using accompanying measurements of CH4 and H2O, and CO2 in dried air, respectively. The best precision of ±0.2‰ for δ13C and of ±4 ppb for CO2 was obtained with an integration time of about 1 hour for δ13C and 2 hours for CO2. The seasonality of CO2 and δ13C was studied by fitting background curves for a complete 2-year period. Peak-to-peak amplitudes of the averaged seasonal cycle were 15.5 ± 0.15 ppm for CO2 and 1.97 ± 0.53‰ for δ13C, respectively. Based on the HYSPLIT Model, air masses were classified into five clusters, with westerly and northeasterly flows being the most and the least frequent, respectively. In the wintertime, northwest and northeast clusters had a higher median level for ΔCO2 and a lower median level for Δδ13C (the difference between observed and background concentrations), likely caused by anthropogenic emissions. In the summertime, air masses from the northwest had the lowest ΔCO2 and the highest Δδ13C. Potential source contribution functions (PSCFs) were used to identify the potential source and sink areas. In winter, source areas for high CO2 mixing ratios (> 75th percentile) were mainly located in northwestern Europe. In summer, areas with high δ13C ratios (> 75th percentile), indicating a carbon sink, were observed in the air from Eastern and Central Poland.

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


The sources and sinks of atmospheric carbon dioxide (CO2) 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 CO2 mixing ratio (IPCC, 2014; WMO, 2016; Xu et al., 2017). Nearly half of the emitted CO2 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 CO2 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 CO2 growth rates, it is crucial to monitor atmospheric CO2 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 CO2 13C) 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 12CO2 and thus discriminate against 13CO2 (~18‰) (Farquhar et al., 1989). In contrast, during oceanic uptake of CO2, there is almost no discrimination (~2‰) (Mook et al., 1974). Since the discrimination by oceans is small, the carbon isotope signature (δ13C) 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 CO2 and δ13C 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 CO2 concentration and its stable isotope ratios (δ13C) at the regional background station Lin’an (LAN) in China to identify the isotopic signature of CO2 sinks and sources. They concluded that during the winter season (Dec–Feb) coal combustion is the major CO2 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 CO2 isotopes (δ13C and δ18O) 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 CO2 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 CO2, δ13C, and δ18O values depending on the origin of air masses and surface residence time. Moore and Jacobson (2015) differentiated sources of atmospheric CO2 by records of CO2 concentration and carbon isotope composition (δ13C) at Evanston, an urban site north of Chicago, the third largest city in the United States.

Measurements of CO2 mixing ratios in combination with stable CO2 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 CO2 mixing ratio and δ13C 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 CO2 and δ13C (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 CO2 mixing ratios and δ13C and (3) identify potential CO2 source and sink locations with the help of backward trajectories and the potential source contribution function.


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.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.


A wavelength-scanned cavity ring down spectrometer (G1101-i, Picarro Inc. USA) was installed in the laboratory on the 4th 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). CO2 mixing ratios and δ13C 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 13C/12C are expressed in per mil (‰), defined by the following Eq. (1) as:


Isotopic values were given relative to the international standard VPDB-CO2 (Vienna Pee Dee Belemnite) (Brand et al., 2010). Measurements of CO2 and δ13C were calibrated using two working standard gases with high and low concentrations (DEUSTE Steininger GmbH) of the certified CO2 mixing ratios and the isotopic compositions (Standard 1: 350.1 ± 0.50 ppm CO2 mixing ratio and –3.28 ± 0.164‰ δ13C; Standard 2: 503.4 ± 0.50 ppm CO2 mixing ratio and –20.03 ± 1.002‰ δ13C) 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 δ13C values (Standard 1: –3.19 ± 0.009‰ δ13C; Standard 2: –20.26 ± 0.01‰ δ13C) 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 N2O (Standard 1 and 2: 0.1 ppm) and CO (Standard 1: 6 ppb; Standard 2: 5.47 ppb), no other species (CH4, H2O, SF6) 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 CO2 and CH4 data measured with a Picarro Envirosense 3000i instrument to correct our measured data (see the section “Correction of CO2”). 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 CO2 and WMO X2004a scale for CH4 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 CO2 measurements from UBA are reported to the WMO World Data Centre for Greenhouse Gases, the respective CO2 data in the section “Correction of CO2” will be referred to as CO2 data of Global Atmosphere Watch (CO2, 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 (CO2: 350.1 ± 0.50 ppm; δ13C: –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 CO2 (left) and δ13C (right) as a function of the integrating time τ, based on a long-term measurement of Standard 1 (CO2: 350.1 ± 0.50 ppm; δ13C: 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) = f0 and S(f)= f2, respectively. The Allan deviation follows a slope of –1/2 up to an integration time of 1 hour for CO2 and of 2 hours for δ13C and then turns over to a slope of 1/2 which defines a drift.Fig. 2. Allan deviation plots for CO2 (left) and δ13C (right) as a function of the integrating time τ, based on a long-term measurement of Standard 1 (CO2: 350.1 ± 0.50 ppm; δ13C: 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) = f0 and S(f)= f2, respectively. The Allan deviation follows a slope of –1/2 up to an integration time of 1 hour for CO2 and of 2 hours for δ13C 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 CO2 and around 2 hours for δ13C. The best achievable precision (1 sigma) was 4 ppb for CO2 and 0.2‰ for δ13C, 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 δ13C and 5–6 ppb for CO2. 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 CH4 that can bias δ13C 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 δ13C on the CO2 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., N2/O2 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 δ13C

The isotopic analyzer reports the peak height of the near-infrared absorption spectrum for the rovibronic transition of 12C16O2 and 13C16O2 with arbitrary labeled C12peak-BookAve and C13peak-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 C13peak-BookAve and C12peak-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).

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


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

Correction of CO2

To quantify and correct water vapor effects on the CO2 measurements, we compared the CO2 mixing ratio for the wet gas stream measured by Picarro G1101-i analyzer (CO2,meas) with CO2 data of the accompanying instrument (CO2,GAW). Although similar seasonal cycles were observed in Fig. 3, CO2, meas mixing ratios only showed a good agreement with CO2, 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 CO2 mixing ratio measured by Picarro G1101-i (in red) without drying system to CO2 measured by Picarro Envirosense 3000i (in black) dried by cold traps.Fig. 3. Comparison of CO2 mixing ratio measured by Picarro G1101-i (in red) without drying system to CO2 measured by Picarro Envirosense 3000i (in black) dried by cold traps.

We corrected the dilution effect using ordinary least squares regression of the CO2 mixing ratio difference between CO2,GAW and CO2, 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

[CO2]calculated = 6.5[H2O] – 2.7 + [CO2]meas                  (4)

[CO2]calculated = 5.6[H2O] – 0.85 + [CO2]meas                 (5)

where [H2O] is water vapor concentration in percent, [CO2]meas is CO2 mixing ratio measured by our analyzer and [CO2]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 CO2 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): R2 = 0.97, p-value < 0.001 and (b): R2 = 0.92, p-value < 0.001).Fig. 4. Scatterplot of water vapor concentrations and difference of GAW and Picarro CO2 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): R2 = 0.97, p-value < 0.001 and (b): R2 = 0.92, p-value < 0.001).

Corrected CO2 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 CO2 and ±0.53‰ for δ13C and combined measurement uncertainty is estimated to be ±0.56 ppm for CO2 and ±0.56‰ for δ13C.

Fig. 5. Histogram of the deviation between spline and true standard value for (a) Standard 1 and (b) Standard 2.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 (Pij) for a given grid cell is then calculated as Pij = mij/nij, in which nij is the total number of trajectory segment endpoints terminating within the grid cell(i,j) over the entire time of measurement and mij 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 (75th percentile of CO2 and δ13C 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 CO2 and δ13C) due to its non-parametric basis. Sturm et al. (2013) used REBS for high alpine CO2 and δ13C 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.


Atmospheric CO2 and δ13C

Time series of CO2 mixing ratios and δ13C 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 CO2 mixing ratios and δ13C values. CO2 mixing ratios simultaneously increase as δ13C values decrease, which is associated with seasonal vegetation activity. The mean peak-to-peak amplitudes are 15.5 ± 0.15 ppm for CO2 and 1.97 ± 0.53‰ for δ13C. The minimum CO2 mixing ratios, as well as the maximum δ13C values, occurred in August due to (preferentially 12CO2) CO2 terrestrial uptake being dominated by the biosphere. The maximum in CO2 mixing ratios occurred in March and the minimum in δ13C values in February when respiration is dominating. In contrast, the peak-to-peak amplitudes of CO2 mixing ratios and δ13C 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 CO2 and 0.60‰ for δ13C (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 CO2 mixing ratios and δ13C from these two sites were different, the minimum and the maximum of CO2 mixing ratios and δ13C were recorded nearly at the same time of the year.

were recorded nearly at the same time of the year.Fig. 6. Time series of hourly mean values of CO2 and δ13C 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 CO2 and δ13C for each season are shown in Fig. 7. The CO2 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 CO2 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 CO2 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 δ13C 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 δ13C.

Fig. 7. Diurnal variations of the mean hourly CO2 mixing ratios and δ13C (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.Fig. 7. Diurnal variations of the mean hourly CO2 mixing ratios and δ13C (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 CO2 mixing ratios and δ13C were different (Sturm et al., 2013). The diurnal peak-to-peak amplitude of CO2 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 δ13C 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 δ13C and CO2, 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.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 CO2 Mixing Ratios and δ13C among the Clusters

In order to study whether measured data of CO2 and its isotopic composition systematically varied with air mass origin, we merged hourly CO2 mixing ratios and δ13C 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 CO2 and δ13C values, providing ΔCO2 and Δδ13C. 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. ΔCO2 and Δδ13C 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 ΔCO2 mixing ratio and Δδ13C for all seasons classified for the five clusters of Fig. 8 for the entire period. ΔCO2 and Δδ13C are filtered for free troposphere air and background concentrations are subtracted.Fig. 9. Boxplots of ΔCO2 mixing ratio and Δδ13C for all seasons classified for the five clusters of Fig. 8 for the entire period. ΔCO2 and Δδ13C are filtered for free troposphere air and background concentrations are subtracted.

In summer, CO2 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 ΔCO2 and Δδ13C. Clusters 4 and 5 had the highest median ΔCO2 values and, correspondingly, the lowest values in Δδ13C.

Using footprint clustering analysis for ΔCO2 and Δδ13C 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 ΔCO2 value and the second lowest value in Δδ13C. Another cluster with high CO2 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 Δδ13C. 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 ΔCO2 and Δδ13C. In summer, C5 exhibited the highest Δδ13C value and the lowest ΔCO2 value, respectively.

Fig. 10. Potential source contribution function plot for ΔCO2 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.Fig. 10. Potential source contribution function plot for ΔCO2 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 CO2 concentration. This was done explicitly for air masses causing extremely high CO2 concentrations (> 75th 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 Δδ13C (> 75th percentile), equally inferred by trajectories, describes the spatial distribution of probable geographical source locations for δ13C (i.e., sinks of CO2) 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 CO2 andδ13C. Fig. 11 clearly illustrates that high values of Δδ13C 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 CO2 uptake.

Fig. 11. (a) Potential source contribution function plot for Δδ13C 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.Fig. 11. (a) Potential source contribution function plot for Δδ13C 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.


We presented a 2.5-year measurement time series of CO2 mixing ratios and δ13C 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 CO2 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 CO2 mixing ratio at UFS and JFJ, respectively, and of 0.4‰ and 0.1‰ for the δ13C 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 CO2 at UFS and JFJ and there was no discernible cycle for δ13C.

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, δ13C measurements and their inverse relationship to CO2 allowed us to identify atmospheric CO2 sinks. Among all clusters, high CO2 mixing ratios in winter (anthropogenic sources) were associated with air masses originating in the Netherlands, the German Ruhr area and Lusatia, whereas high δ13C values and low CO2 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 CO2 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 CO2 (such as the data presented in this publication) or records of distinct CO2 variations (e.g., a sudden enhancement in atmospheric CO2) 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).


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. 

Don't forget to share this article 


Subscribe to our Newsletter 

Aerosol and Air Quality Research has published over 2,000 peer-reviewed articles. Enter your email address to receive latest updates and research articles to your inbox every second week.

Latest coronavirus research from Aerosol and Air Quality Research

2018 Impact Factor: 2.735

5-Year Impact Factor: 2.827

SCImago Journal & Country Rank

Aerosol and Air Quality Research (AAQR) is an independently-run non-profit journal, promotes submissions of high-quality research, and strives to be one of the leading aerosol and air quality open-access journals in the world.