Application of a Common Methodology to Select in Situ CO 2 Observations Representative of the Atmospheric Background to an Italian Collaborative Network

: We describe and implement a data selection algorithm aimed at identifying background atmospheric CO 2 observations from in situ continuous measurements. Several selection criteria for detecting the background data have been developed and are currently used: the main objective of this work was to deﬁne a common methodology to extract the atmospheric background signal minimizing heterogeneities due to the use of different selection algorithms. The algorithm used in this study, (BaDS, Background Data Selection) was tested and optimized using data (from 2014 to 2018) from four Italian stations characterized by markedly different environmental conditions (i.e


Introduction
The Mediterranean basin is considered a global hot-spot region for climate change and air-quality. CO 2 is the single most-important anthropogenic greenhouse gas (GHG) in the atmosphere, accounting approximatively for 62% of the anthropogenic radiative forcing by long-lived GHGs. Fossil fuel and cement production emissions are the main drivers for the increasing atmospheric CO 2 dry mole fraction in the global atmosphere [1].
In response to this, regulation and emission trading schemes have been adopted at international, national, and city levels to reduce GHG emissions, while responding to the needs for economic and societal development. The successful implementation of GHG emission reductions must be guided by reliable scientific evidence. While satellite observations and ground-based remote sensing measurements have gained importance in recent decades [2][3][4], global and regional networks of in situ GHG observations still largely represent the reference information for GHG monitoring, emission and variability studies [1,5,6]. These networks are characterized by significant efforts dedicated to the harmonization and increased comparability of GHG observations [7][8][9].
The analysis and interpretation of atmospheric mole fraction time series play a fundamental role in assessing long-term trends and interannual growth rates of CO 2 . To this aim, it is pivotal to extract the information related to the atmospheric background from the data series, i.e., data "representative of the global atmosphere not affected by local conditions" [10]. The identification of the background signal in the GHG time series is also necessary to isolate data, which yields information related to processes occurring at regional scales, i.e., regional fluxes [11]. Usually, remote islands or high-mountain sites are indicated as suitable locations to trace the background atmospheric levels of atmospheric species. However, as shown in previous studies [12][13][14][15][16][17][18][19], even at island or high-mountain sites, the influence of direct anthropogenic and natural emissions cannot be neglected.
For these reasons, selection criteria have been developed and used to identify measurements directly affected by emission or removal processes, and to retain data that are well representative of the atmospheric background [10,20,21]. For CO 2 , some of these approaches are based on the following common strategy to identify non-background data: a polynomial equation is used to remove the long-term CO 2 trend, then a harmonic function is applied to reproduce the seasonal variation, and (finally) residual variations are determined. To separate the seasonal component from long-term trend and irregular variations (i.e., time series decomposition), different approaches can be used: "HPspline" [22][23][24], CCGCRV (Carbon Cycle Group CuRVe) [25][26][27], Seasonal Trend decomposition (STL) [28] or least square fitting methods [29]. However, to infer information related to latitudinal gradients, regional fluxes, and large-scale temporal variability, it is necessary to apply these decomposition techniques to a subset of data that are representative of the atmospheric background. Different methodologies exist to extract the background data from the time series of atmospheric constituents: the robust extraction of a baseline signal (REBS and 2-D-REBS) [30,31], the adaptive selection of diurnal minimum variation (ADVS) [21], the coefficient of variation (COV) [32] and the standard deviation of the background (SD) [33]. A comparison study between COV, SD and REBS approaches applied to four stations with different environmental conditions contributing to the Integrated Carbon Observation System, ICOS (http://www.icos-cp.eu, accessed on 9 December 2020), European Union research infrastructure, was carried out by El Yazidi et al. [34]. According to their results, the SD and REBS methodologies are the most flexible for spike detection across different sites, but REBS tends to underestimate spikes when these are very high. Instead, the SD methodology appears to better detect data not belonging to the baseline. Giostra et al. [20] proposed the use of a statistical filter that identifies the data baseline, i.e., the subset of data in a time series which are not significantly affected by non-well-mixed contributions. By assuming that a well-mixed species follows a Gaussian distribution, the observed CO 2 molar fraction probability density function (PDF) should be fitted with the sum of a Gaussian PDF and a Gamma PDF (the latter describing the non-well-mixed data). All the data for which the corresponding Gaussian is higher than the Gamma PDF are assigned to the baseline. Among the different approaches, the method of Thoning et al. [25], initially developed for the measurements of Mauna Loa Station, has been widely used. It consists of two steps which analyze the variability of the hourly CO 2 data (i.e., their standard deviation) and the difference between consecutive hourly means. Then, an iterative algorithm removes any values which differ from the weighted spline curve by more than a given threshold. Due to its simple usability and the lack of any strong pre-requisite regarding the measurement site environmental conditions, it has been applied to other remote stations, for example, Schauninsland, a site in the southwest of Germany [35], Junfgraujoch (Switzerland), Puy de Dôme (France) [36], and four WMO/GAW (Global Atmosphere Watch programme of the World Meteorological Organization) stations in China [37].
In this study, we tested the application of the background data selection algorithm presented by Apadula et al. [19], called BaDS (Background Data Selection) to five different datasets of in situ continuous CO 2 observations in Italy, belonging to a collaborative network. This algorithm is based on the method of Thoning et al. [25], and has some common points with the methodology already used by Cundari et al. [38]. With respect to [19], we paid particular attention to applying a robust identification of the threshold value used to discriminate the CO 2 variability (i.e., impacted by large-scale vs. regional/local phenomena). The aim of this study is to evaluate the ability of the BaDS methodology in extracting the information related to the background variability from CO 2 time series recorded at fixed monitoring stations with substantially different characteristics: mountain, coastal and marine sites. This will permit verifying if, by using a common approach, we can obtain a consistent view of the background CO 2 variability at the different sites of the observational network. Moreover, the identification of measurement periods representative of the atmospheric background also has the co-advantage of permitting identification of the periods which are above or below the atmospheric background, thus allowing for the specific characterization of the measurement sites. This latter point is particularly important when assessing the capability of a measurement site in catching signals useful for the quantification and spatial allocation of anthropogenic or natural emissions/sinks.
Even if most of these observatories have longer datasets [19,29,39], only a subset of 5 years has been used to test the methodology. Indeed, a long-term analysis is beyond the scope of this paper. Section 2 reports a description of the measurement sites, of the experimental setup, and of the background data selection algorithm (BaDS). The daily and the seasonal variability of the different CO 2 datasets, with an analysis of the impact of the data selection algorithm are discussed in Section 3. In Section 3, we also provide an example of the BaDS application for the investigation of the interannual variability of the seasonal cycles. A discussion on the method, with implications and future research directions, is provided in Section 4. Conclusions are drawn in Section 5.

Measurement Sites
This study considered near-surface CO 2 measurements from four laboratories (CAMM, ENEA, National Research Council of Italy-Institute for Atmospheric Sciences and Climate (CNR-ISAC) and RSE S.p.A.) at the four atmospheric observatories, the location of which is shown in Figure 1. The geographical locations of these observatories (spanning from north to south across the Italian peninsula, separated by more than 1200 km and more than 10 • latitude, and from sea level up to 3480 m a.s.l.) provides the opportunity to investigate the spatial and temporal variability of greenhouse gases over Italy and the central Mediterranean basin. In the classification of Henne et al. [15], Lampedusa (LMP) was classified as a generally remote site and Mt. Cimone (CMN) as weakly influenced by regional emissions. Plateau Rosa (PRS) can be considered similar to the near site of Jungfraujoch (JFJ, at 3571 m a.s.l. on the Switzerland Alps), classified by Henne et al. [15] as a mostly remote site.
The main aim of this work was to test and implement a methodology for background data selection to verify its performance and usefulness for obtaining a coherent baseline dataset. We considered only a limited period of observations (see Table 1    The PRS station is part of the Testa Grigia Observatory (45 • 56 N 7 • 42 E, 3480 m a.s.l., http://oasi.rse-web.it/, accessed on 9 December 2020), a high mountainous site in the Italian western Alps far from potential pollutant sources. It is situated over an Alpine glacier, often above the planetary boundary layer (PBL), devoid of soil and vegetation. These characteristics make the station suitable for background measurements of greenhouse gases. Incursion events of polluted air masses from the European plains rarely occur at the PRS site [17]. More information on the measurement site can be found in [19]. PRS contributes to WMO/GAW as a regional station. This site is under the labelling processes to join the ICOS Research Infrastructure. The dataset used in this work is available at the World Data Center for Greenhouse Gases (WDCGG) by the WMO/GAW (see https://gaw.kishou.go.jp/, accessed on 9 December 2020). CMN (44 • 12 N, 10 • 42 E, 2165 m a.s.l.) is the highest peak of the northern Apennines. This mountain range separates the European continent from the Mediterranean regions. The site is characterized by a free horizon for 360 • and is far from big cities and important industrial areas (Bologna and Florence are located more than 50 km away and are 2 km below the station). CMN is considered to be well representative of the free tropospheric condition, even if direct transport of air masses from the regional PBL can be systematically observed, especially during the warm months [13,40]. In particular, as observed by Colombo et al. [13], during May-September, the CO 2 variability is significantly affected by the diurnal vertical transport of air masses, which produces a decrease in the atmospheric molar ratio due to photosynthetic uptake by vegetation. The atmospheric observatory, composed by the Italian Air Force meteorological station (located at the mountain top) and the "O. Vittori" CNR observatory (located about 50 m below the mountain top, see http://cimone.isac.cnr.it/, accessed on 9 December 2020), is a WMO/GAW global station [41]. The "O. Vittori" observatory was approved as an ICOS atmospheric class 2 site in November 2018. The dataset considered for the CAMM laboratory was extracted by WDCGG. The dataset for the CNR laboratory was extracted from the ICOS dataset for the period May-December 2018 (Cristofanelli et al., 2020 "ICOS CO 2 release", see https://data.icos-cp.eu/, accessed on 9 December 2020) [42], while for the period January-April 2018 an "interim" dataset (for which the ICOS protocols were not fully implemented) was considered.
The CGR station (37 • 34 N, 12 • 40 E, 15 m a.s.l.; http://www.i-amica.it/i-amica/ ?page_id=632 accessed on 9 December 2020) is managed by CNR and is located on the south-western coast of Sicily, overlooking the Sicily channel, in the central Mediterranean. The observatory has been operational since 2014 and contributes to WMO/GAW as a regional station. Thanks to its semi-background location, various studies about the atmospheric processes influencing the variability of traces gases and aerosol [43,44] and the characterization of the Mediterranean Sea background conditions have been carried out [45]. This measurement site is strongly affected by the land/sea breeze regime, which leads to the advection of background air masses from offshore during the daytime, and polluted air masses from inner Sicily during the nighttime [43]. This observatory is located within the CNR scientific campus at Torretta Granitola (a small village 12 km from Mazara del Vallo, 12,000 inhabitants). The shelter where the observatory is hosted is situated on a cliff (15 m high) overlooking the coastline. The dataset used for this analysis is available at the WCDGG by WMO/GAW.
The ENEA Atmospheric Observatory LMP (35 • 31 N, 12 • 38 E, 45 m a.s.l.; http://www. lampedusa.enea.it/ accessed on 9 December 2020) is located on the island of Lampedusa. Lampedusa is a small island (about 22 km 2 surface area, 6000 inhabitants) more than 100 km from the African continent, and more than 200 km South of Sicily. The island is flat, with sparse vegetation, and limited local emission sources. In the study by Henne et al. (2010) [15], the site was characterized as the most remote (together with Mace Head, in Ireland), out of 34 different European considered sites. The atmospheric observatory is located on a 45 m plateau, on the north-eastern coast of the island. Several parallel measurement programs aimed at the long-term monitoring of climate and investigating climate-related processes are carried out at the observatory [46]. Previous studies [47] have shown that the Lampedusa CO 2 record is characterized by a very limited diurnal cycle; this is due to the very limited impact of local sources and sinks, and the reduced daily variability of the planetary boundary layer in the open sea. Since 1992, this observatory has been contributing as a WMO/GAW regional station. LMP also contributes to the US National Oceanic and Atmospheric Administration Global Greenhouse Gas Reference network (https://www.esrl.noaa.gov/gmd/ccgg/ggrn.php accessed on 9 December 2020). On May 2020, the LMP station was labelled as a "class-2" atmospheric station within ICOS. The dataset considered for LMP was produced by ENEA and is available at the WCDGG by WMO/GAW.

Background Data Selection Algorithm (BaDS)
We tested and adopted a methodology (BaDS-Background Data Selection) originally applied at PRS [19] and based on the classical data selection used at the Mauna Loa observatory [48]. BaDS is based on the assumption that background observations are characterized by a low temporal variability because they should not be affected by nearby emissions or removal processes. In this work, BaDS was applied to time series of mean hourly CO 2 values. Before BaDs application, short data gaps due to calibrations or outages were filled using a simple linear interpolation. The procedure, summarized below, consisted of 3 steps.
Firstly, hourly mean values were selected as a function of their respective standard deviation, calculated from the 1 min mean values: only data with hourly standard deviation <1 ppm were retained. The selected data were then used to: (a) calculate a running median on 504 h (21 days), and (b) calculate the average difference between consecutive hourly mean values (defined as the "S" parameter).
Secondly, the differences (∆C) between the retained hourly values and the corresponding 504-h running median values were calculated. A confidence interval equal to n * S, where n was a variable parameter defined by the user, was calculated. All data points with |∆C|> n * S were discarded. Basically, through this step all the outliers were rejected from the dataset.
Thirdly, 48 h running means were calculated using the data retained by the previous step. The differences (δC) between hourly mean values (with σ < 1 ppm) and the corresponding 48 h running means were then calculated. All data points with |δC|> n * S were flagged as "non background".
It should be noted, that during Step 1 (i.e., calculation of the reference running median and "S" parameter), only nighttime observations (i.e., under mountain breeze regime, considered less affected by local/regional contributions) were considered from May to September for the CMN-ISAC and CMN-CAMM datasets. Similarly, only observations during daytime (i.e., under sea breeze regime) were considered for CGR. In this way, the impact of vegetation uptake under the thermal upward wind regime (at CMN) and the impact of anthropogenic emissions transported under land breeze circulation (at CGR) have been minimized. The complete datasets were considered in the analysis for PRS and LMP.
A critical aspect for this method was the determination of the value of n. In [19], the value of n is determined based on expert judgement and good knowledge of the site data. Since data from different sites are considered in this study, we have developed an improved (and less arbitrary) method for the determination of n.
A sensitivity study was carried out to define, for each dataset, the optimal values of the parameter n. The procedure was applied for n ranging from 1 to 15. Two separate tests have been performed on each dataset. The first test aimed to assess the influence of the parameter n on the number of non-background data, in order to verify the method selection capability and to avoid creating under representative datasets. The second test was performed to evaluate the effect of the parameter n to the diurnal CO 2 cycle amplitude under background conditions. As denoted by the analysis of the original dataset, each measurement site was characterized by a typical diurnal variation cycle whose amplitude is a function of the proximity of source/sink regions ( Table 2). The lowest median amplitude characterized the diurnal cycle at PRS and LMP (0.78 and 0.58 ppm, respectively). Despite its elevation, CMN is affected by a mountain breeze regime which favors the advection of air masses from the PBL during the daytime of warm months (May-September), resulting in a larger median amplitude of the CO 2 diurnal cycle (3.80 ppm for CMN-CAMM and 4.47 ppm for CMN-ISAC). The largest median amplitude of the diurnal cycle was found at CGR (8.9 ppm), due to the influence of the marine breeze regime which advects (relatively clean) marine air masses during the daytime and (polluted) air masses from inland during the nighttime. Since the CO 2 diurnal variability is expected to be minimized at the measurement sites for background conditions, the second test was performed to identify the optimal n values which reduce this amplitude. We believe this provided a useful criterion for the determination of the optimal n value to be used to detect the background data. To this aim, for each station, the median amplitude of the diurnal cycle (calculated over the whole dataset for each site/laboratory) was calculated as a function of the background dataset obtained by applying different n values. The results of both sensitivity tests are shown in Figure 2. Table 2. Median amplitude of the diurnal variation (expressed in ppm) and percent of data removed (in brackets) for each measurement site. The values of S and of the retrieved optimal n parameters are reported for each site. In general, for CMN-ISAC and CMN-CAMM, the minima in the median diel CO 2 cycle amplitude were found at n = 5 and n = 6, respectively. For LMP and PRS, only a limited dependency was found on the n value, while the strongest dependency was found for CGR with the lowest median CO 2 cycle amplitude at n = 1. Thus, setting n = 6 and n = 10 for PRS and LMP represented a good compromise to minimize the (already low) diel CO 2 cycle amplitudes and to retain a usable fraction of data as representative of background.

Dataset
Results for CGR showed a different behavior with respect to the other datasets, with a somewhat larger daily cycle amplitude, and a fraction of rejected data which denoted only a limited dependency on the n value. Interestingly, differences were evident between CMN-CAMM and CMN-ISAC as concerning the dependence of the diel CO 2 cycle amplitude and the fraction of data retained as "background" as a function of the adopted n values. This is related to the different temporal data coverages of these dataset. Indeed, the optimal value of n = 5 was obtained for both CMN-CAMM and CMN-ISAC when only data for one year (2018) were considered.
As is apparent from Table 2, it is evident that the application of BaDS was effective in reducing the diurnal cycle, by decreasing the impact of local/regional emissions or removal processes at all the measurement sites. In particular, the amplitude of the diurnal cycle was strongly reduced at CMN (for both ISAC and CAMM dataset a 4-fold reduction was obtained) and CGR (8-fold reduction). At these sites a significant fraction of the data (considering the whole datasets) was rejected by BaDS as non-background: 38 and 20% at CMN-ISAC and CMN-CAMM, 65% at CGR. The impact of BaDS on the median diurnal cycles at PRS and LMP was less evident (even if still discernible at PRS) due to the lower impact of local/regional emission/removal processes at these sites, as also testified by the low fraction of removed record (11% for PRS and 8% for LMP). Based on these results, LMP can be considered as being well representative of the baseline conditions of the central Mediterranean basin, and its comparison with the dataset from other sites can provide hints to discuss the CO 2 variability over Italy.
The optimal overall threshold used in the selection, equal to n * S, is very similar for PRS, LMP, and CMN (between 2.16 and 2.6 ppm). A much smaller threshold was needed at CGR to isolate the background conditions from the totality of observations.
A further sensitivity test was performed with the aim of assessing the impact to the BaDS results of applying, for each measurement site, different n values calculated on a seasonal basis. Table 3 reports, for each measurement site, the fraction of data discarded and not representative of the atmospheric background when the seasonally-derived n values were used in BaDS. In respect to the annually-derived n values (Table 2), the use of the seasonally calculated n values led to differences in the fraction of discarded data (Tables 3 and 4) that ranged from −9.7% at CGR to +3.8% at PRS. The agreement between the two approaches was 90% (LMP), 83% (PRS), 74% (CMN-CAMM), 55% (CMN-ISAC), and 30% (CGR), defined as the percentage of data that was retained as background simultaneously by using the annual and the seasonal n value. For all the measurement sites (but not for PRS), the median amplitude of the CO 2 diel cycle was minimized when adopting the annual n value. Only for PRS was there no difference between the two approaches: this was probably due to its remote location, more consistently representative of the atmospheric background and less affected by local source/sink through the seasons of the year. Thus, based on these results, we decided to use the n value obtained on the annual basis in the BaDS. Table 3. Median amplitude of the diurnal variation (expressed in ppm) and percentage of data removed (in brackets) for each measurement site using seasonal values of n. The values of S and of the retrieved optimal n parameters are reported for each site.

Analysis of CO 2 Diurnal Variation and Impact of BaDS Application
With the purpose of assessing the impact of the BaDS selection on the typical CO 2 diurnal cycle at each measurement site, we calculated the seasonal (winter: DJF, spring: MAM, summer: JJA, autumn: SON) average diurnal variations of the detrended CO 2 hourly values for the original dataset and for the background selection (after the BaDS application). For each dataset, the detrended values have been obtained by subtracting from the hourly values the corresponding value of the fit obtained applying the Equation (1) from [29]: where y is the monthly CO 2 molar fraction, c 1 and c 2 are constant values, m is the incremental number of months, T 1 and T 2 represent semi-annual and annual periods (6 months and 12 months, respectively), ϕ 1 and ϕ 2 are the semi-annual and annual phases, respectively. The parameters obtained from the fit are shown in Table 5. Results for the calculated seasonal average diurnal cycles are shown in Figure 3. For the original datasets ( Figure 3, top), a small diurnal cycle amplitude (less than 1 ppm) characterized PRS and LMP during all seasons. A significant diurnal cycle (with lower values during daytime) was evident for CMN in summer (average amplitude: 5.5 ppm) for both the ISAC and CAMM datasets, in good agreement with the summer variability observed by Ref. [13]. The diurnal cycle at CMN, as previously discussed, was attributed to the systematic thermal transport of PBL air masses up to the mountain peak during daytime; the valley breeze brings air masses that are depleted in CO 2 to the measurement site due to the vegetation sink, favored by the presence of woods on the slopes of the mountain [13]. The summer PRS diurnal cycle was about 2.5 ppm, in agreement with Schmidt et al. [35], who suggested that the amplitude of the summer diurnal cycle in mountain sites is inversely proportional to the altitude. Small differences appeared between CMN-ISAC and CMN-CAMM, which were primarily related to the different data coverage (only one full year for CMN-ISAC), even if an impact from the different sampling heights [50] cannot be completely neglected. The largest diurnal cycle amplitudes were observed for CGR, with seasonal amplitudes ranging from 8 ppm in winter to 12 ppm in summer. As discussed above, this was linked to the influence of the marine breeze regime and the influence of regional sources and sinks. Fires from agricultural residues and uncontrolled landfills that occurred within the territory might have provided additional anthropogenic CO 2 sources that contributed to produce a marked daily cycle. As also reported by Cristofanelli et al. [43] and, for other coastal sites in Sicily, by Perrino et al. [51], during all the seasons CGR showed strongly higher pollutants concentrations than LMP during nighttime (due to the advection of air masses from inner Sicily, affected by anthropogenic sources and by the effect of soil and vegetation as a sources of respiratory CO 2 ) and more comparable values during daytime (due to the advection of marine air masses representative of the Mediterranean basin PBL). The diurnal cycles at CGR were significantly reduced in amplitude for the background data selected by BaDS (Figure 3, bottom). This was true also for the summer data at PRS and CMN. For CGR, the shape and average values of the background diurnal cycle approached those observed at LMP. Even if the difference from LMP was strongly reduced, at CGR, signals of regional anthropogenic emissions were probably still reflected by the higher values (especially during summer nighttime). Especially during spring, the CGR background average diurnal cycle approached that observed at LMP, stressing the ability of BaDS in selecting measurement periods representative of the atmospheric background even at a measurement site strongly exposed to anthropogenic emissions.
For the background datasets, lower CO 2 values were typically observed at PRS (−1.5 ppm, on average) with respect to LMP. During winter the tropospheric CO 2 is characterized by a vertical gradient with the highest CO 2 levels near the surface due to biogenic and anthropogenic emissions [52], and it is reasonable to suggest that PRS is more affected by free tropospheric air masses compared to LMP. During the vegetative season, it is likely that the lower CO 2 values observed at PRS and CMN were tracing a background latitudinal gradient between continental and oceanic regions in the Mediterranean basin. With respect to LMP, slightly higher values (+0.5 ppm) were instead observed for CMN-CAMM during the central part of the day in winter, possibly indicating a residual role of upward thermal transport of air masses from the PBL. The opposite situation was observed in summer daytime, when the mountain sites (especially CMN) were more exposed to the transport of air masses from the PBL (Figure 3 bottom). Figure 4 shows the CO 2 monthly mean time series at the different sites from 2014 to 2018 by considering the original datasets, and the monthly time series of the background CO 2 dataset, as obtained by the BaDS application. Apart from the multi-year increasing tendency, the seasonal variations of CO 2 through the years is evident, with maxima in the first months of each year, a minima during summer, and an amplitude of the annual cycle of about 10-15 ppm, depending on the year and the measurement site. The seasonal CO 2 cycle was a combination of different contributions: biosphere emission and removal processes, anthropogenic emissions and atmospheric transport (on different temporal and spatial scales). The different time series selected with BaDS were more consistent with each other than the original datasets. Only the CMN-ISAC dataset was still characterized by the appearance of the high CO 2 values during February 2018, significantly higher than those observed for CMN-CAMM. Additionally, the application of BaDS based on the seasonalderived n parameter (Section 2.2) does not allow for the detection and removal of high CO 2 values in February 2018, which were related to the occurrence of specific multi-day peak events observed at CMN-ISAC (please note that these events were not present in the CMN-CAMM dataset due to a stricter flagging strategy adopted at this laboratory:

Analysis of CO 2 Time Series and Impact of BaDS Application
i.e., data points that consistently deviated from the rest of the dataset were rejected by a data screening process performed by station operators on the raw data). This would suggest that BaDS (at least with the parameter configuration indicated in Section 2.2) was not able to correctly identify the multi-day pollution episodes occurring in this month as non-background data. The methodology performed well during the winter season and also during pollution events; however, elevated CO 2 cases lasting for numerous days remain problematic with this algorithm.

Analysis of CO 2 Seasonal Cycle and Impact of BaDS Application
To better characterize the impact of the BaDS application on the determination of the CO 2 seasonal cycle at the different sites, we calculated the average monthly mean values of the detrended CO 2 (Figure 5, left) over the whole periods of data availability by using the same methodology used for the analysis of the diurnal cycle (see Section 3.1). This has been done for the original dataset and for background dataset. In the original dataset, we observed similar seasonal cycles for PRS and LMP: both sites showed an annual minimum in August and a maximum in March. CMN (both ISAC and CAMM dataset) was denoted by a different shape of the average annual cycle with a more marked (and fast) decrease in May-July with respect to PRS and LMP, reflecting a more marked signature of the continental photosynthetic activity. This led to the occurrence of a broader summer minimum at CMN than at PRS and LMP. Marked differences existed between the CO 2 average values during the cold months (February-March and November-December) between ISAC and CAMM dataset at CMN. As already explained in Section 3.2 this was related to the different flagging strategies adopted at the two laboratories. The CO 2 seasonal cycle was significantly larger at CGR than at the other sites. As discussed above, especially during nighttime (under the land-breeze regime), CGR was also more directly affected by anthropogenic CO 2 emissions due to fires from agricultural residues and uncontrolled burning of landfills that occurred within the regional territory. Although characterized by more vegetation than LMP, the monthly difference at CGR with respect to LMP did not decrease in summer. This may be due to the hot summer that characterizes these territories, favoring a more dry and arid soil, affecting the vegetation and thereby reducing CO 2 uptake. In addition to this lack of vegetative uptake, an increase in CO 2 emissions due to touristic activity cannot be neglected, as suggested by the seasonal increase observed for other anthropogenic co-emitted compounds (e.g., NOx, CO, SO 2 ) [43]. The differences among the different seasonal cycles were minimized after BaDS application ( Figure 5, right). If we exclude the CMN-ISAC dataset, for which high values were still observed in January and February (even if notably reduced in respect to the original dataset), the differences between the different measurement sites are reduced to 2.5 ppm (from October to March) and to less than 1 ppm (July). In particular, the decreasing phase of the annual cycle from April to August was similar for all sites, with the CMN dataset approaching PRS and LMP. The outstanding agreement of the CGR background to the other sites demonstrates the ability of the BaDS algorithm to extract the background signal even if the measurement site was strongly impacted by local/regional emissions.

Use Case of the BaDS Application: Investigation of Interannual Variability of the Seasonal CO 2 Cycle
With the purpose of providing an example of a possible application of the BaDS algorithm, we analyzed the interannual variability of the seasonal cycles of background CO 2 data at the four observatories belonging to the Italian collaborative network. In order to quantify annual and semi-annual components, Equation (1) from [29] was fit to the CO 2 monthly mean values for the background data selection with BaDS ( Table 5).
The largest annual cycle amplitude (calculated as |2xB|) was found for CMN-CAMM (12.0 ppm) and the lowest for LMP (8.8 ppm). Intermediate values were found for PRS (10.2 ppm) and CGR (9.8 ppm). This clearly confirmed the lower regional/continental influence at LMP (less affected by vegetation uptake in summer and anthropogenic emission in winter) as well as a larger exposure of CMN to regional/continental emission/removal processes. The values of the seasonal cycle amplitude showed a positive gradient with latitude, indicating a decreasing influence of continental sources/sinks moving from northern to southern sites. Figure 6 shows the annual cycles for the different datasets and the different years. An interannual modulation of the summer minimum is evident. Consistently, all the dataset shows a broader summer minimum in 2015 and 2017, while a narrower minimum appearing in August in 2014, 2016, and 2018. This phenomenon was related to a lower decrease in CO 2 values in August during years 2015 and 2017, which might suggest less efficient uptake by vegetation during these years and/or increased CO 2 emissions, possibly linked to forest fires. Murayama et al. [53], by examining the role of atmospheric circulation variability in affecting the interannual variability of the observed atmospheric CO 2 growth rate in the Northern Hemisphere, suggested that the variability of the summer atmospheric CO 2 decrease can be the result of the synergistic action of variability in biospheric sources and air-mass transport, as deduced by analysis of the North Atlantic Oscillation (NAO) and the Pacific-North America (PNA) modes.

Discussion
In order to discriminate CO 2 observations as representative and non-representative of atmospheric background conditions at the four measurement sites spanning from the Alps to the central Mediterranean, we implemented the methodology presented in Ref. [19] (BaDS). BaDS is a pre-conditioning algorithm for extracting time periods representative of the atmospheric background conditions from time series of CO 2 derived by in situ continuous observations. Due to differences in the environmental characteristics of the measurement sites, a specific strategy was used to identify the most appropriate tuning of the parameters involved in the data filtering.
We showed that the implemented methodology has the capability to reduce, from the original time series, the impact of local/regional CO 2 fluxes, thus extracting a signal more related to the atmospheric CO 2 background. In particular, it also performed well at CGR, a measurement site characterized by a complex interaction of anthropogenic emissions, vegetation uptake and a sea breeze wind regime. After the BaDS application, the consistency of the monthly mean CO 2 values strongly increased among the considered measurement sites. A notable exception was the CMN-ISAC dataset which showed a divergent average value in February 2018. For this specific month, BaDS did not appear to identify a long-lasting pollution event (corroborated by the simultaneous increase in CO and NO 2 ) which affected the measurement site. This episode was characterized by an overall increase in CO 2 for several days: it is likely that the lack of high-frequency variability prevented the attribution of this episode to the non-background conditions.
The comparison among the CMN-ISAC dataset (for which one full year of data are available) and the CMN-CAMM dataset (for which three full years of data are considered), suggested that the length of the time series can affect the BaDS output, thus creating a certain degree of uncertainty in the selection process. As shown by Table 2, the different lengths of the time series lead to different optimal "n" values and different amplitudes of the diurnal CO 2 cycle for the background dataset. When only 1 year of data were considered for CMN-CAMM, more consistent results were found with CMN-ISAC. Based on preliminary tests carried out by considering an extended CMN-ISAC time series (2018-2019, here not shown), we argue that the background data selection increases in robustness when longer time series are used.
With the aim of evaluating the impact of tuning the "n" parameter in respect to the original BaDS setting used in Ref. [19], we compared BaDS results for the PRS station using n = 7 (according with Ref. [19]) and n = 6 (according with our optimal tuning). As regards our optimal tuning, for the background selection, the original BaDS application lead to a larger amplitude of the diurnal CO 2 cycle (0.42 ppm) with a lower number of discarded data (8%). This would suggest that, in respect to the "expert" subjective setting proposed by Ref. [19], the optimal "tuning" is more effective in leaving out the residual signals that are unrelated to background conditions. However, the parameter "n" selected in this study, was obtained by applying the methodology to a shorter time series in respect to Ref. [19] and this might have led to a different optimal value for "n".
The BaDS capability to identify the CO 2 measurement period representative of the atmospheric background in different environmental conditions (from the high mountain to a coastal site influenced by anthropogenic emissions, and to a remote marine site) may allow for further studies concerning the long-term CO 2 variability by using a larger number of measurement sites. However, as demonstrated in this work, a fine-tuning method is recommended to set the proper value of a key parameter (n) for optimizing the selection procedure: this optimal value is dependent on the measurement site and on the length of the time series. This implies that the adoption of a full common set-up for all the measurement sites still remains a challenging goal.

Conclusions
In this work, we assessed the capability of a data selection methodology (BaDS) for identifying the observations representative of the atmospheric background from the time series of near-surface CO 2 (hourly mean values) at four permanent observatories in Italy from 2014 to 2018: Plateau Rosa, Mt. Cimone, Capo Granitola, and Lampedusa. These observatories are located in different environmental conditions and latitudinal/altitudinal ranges, and are managed by different institutions. In addition to providing hints about the ability of BaDS in detecting background data, the systematic comparison of the CO 2 dataset recorded at these observatories is also relevant to investigate the dependence of the observed CO 2 variability as a function of altitude and latitude in Italy. In respect to the original BaDS application performed by Apadula et al. [19], we implemented a more objective definition of the key parameter n which tunes the confidence interval to retain data points as representative of the atmospheric background. After this specific tuning of the key parameter n, 94 and 89% of the data were selected as representative of the background conditions at LMP and PRS, respectively. A larger fraction of data was determined as non-background at CMN, with a maximum occurrence of non-background data at CGR (65%). A more coherent diurnal and seasonal CO 2 evolution among the various datasets was obtained after BaDS application.
Considering the original dataset (i.e., without BaDS application), a low average diurnal amplitude (less than 1 ppm) characterized PRS and LMP during all seasons, while a significant diurnal cycle was evident at CMN during vegetative growing months (mostly due to the vertical transport of air masses depleted in CO 2 due to vegetation photosynthesis) and CGR (due to the transport of air masses affected by anthropogenic emissions during nighttime). After BaDS application, the diurnal CO 2 cycle was strongly reduced at CMN and CGR. At PRS and LMP the diurnal cycle amplitude remained lower than 1 ppm. Due to the high altitude at PRS, and the very remote location at LMP, these sites were less affected by emission or removal processes at regional scales. This was also shown by the larger seasonal cycle observed at CMN (12 ppm) with respect to the other sites, and the lowest seasonal amplitude at LMP (8.8 ppm). Based on these results, the LMP and PRS dataset can be considered well representative of the baseline conditions of the central Mediterranean basin, while CMN (especially in summer) and CGR are suitable locations to systematically assess the influence of emissions and sinks occurring within the regional PBL.
Based on our results, after a correct tuning of the setting parameter "n", the BaDS methodology appeared to be able to identify data affected by sources/sink of CO 2 and to extract the information about the background variability at the selected sites which display very different environmental characteristics. Here, we tested BaDs over a set of CO 2 datasets observed in Italy. Its usage at other measurement sites may be the object of further investigation.
Finally, as a further application of the method, we showed how we can retrieve robust information on the regional variability of the annual cycle. Despite large differences in the original datasets, a coherent evolution, with broader summer minima at all sites during 2015 and 2017, was obtained from the analysis of the background data. In this paper, the performance of BaDS methodology have been evaluated on a set of data produced by the Italian network, but it will be interesting to verify the application at different monitoring sites, such as the European infrastructure network ICOS, and evaluate the performance in respect to other methodologies present in the literature (e.g., Refs. [20,[30][31][32][33]) in a future work.  Data Availability Statement: Data available in a publicly accessible repository. The data presented in this study for CMN-ISAC are openly available in ICOS RI, licensed under CC4BY at PID 11676/a5Jn7fKEo4dz8f4pKmqrQhPM. Data for CMN-CAMM, CGR, LMP and PRS are available in a publicly accessible repository that does not issue DOIs. Publicly available datasets were analyzed in this study. These data can be found here: https://gaw.kishou.go.jp/ accessed on 9 December 2020.