Statistical Analysis of the CO2 and CH4 Annual Cycle on the Northern Plateau of the Iberian Peninsula

Outliers are frequent in CO2 and CH4 observations at rural sites. The aim of this paper is to establish a procedure based on the lag-1 autocorrelation to form measurement groups, some of which include outliers, and the rest include regular measurements. Once observations are classified, a second objective is to determine the number of harmonics in order to suitably describe the annual evolution of both gases. Monthly CO2 and CH4 percentiles were calculated over a six-year period. Linear trends for most of the percentiles were around 2.24 and 0.0097 ppm year−1, and the interquartile ranges of residuals calculated from detrended concentrations were 6 and 0.02 ppm for CO2 and CH4, respectively. Five concentration groups were proposed for CO2 and six were proposed for CH4 from the lag-1 autocorrelation applied to detrended observations. Monthly medians were calculated in each group, and combinations of harmonics were applied in an effort to fit the annual cycle. Finally, adding annual and semi-annual harmonics successfully described the cycle where one step was observed in the concentration decrease in spring, not only for high CO2 percentiles but also for low CH4 percentiles.


Introduction
Most greenhouse gas measurements at remote sites and the related research focuses on determining gas concentration, annual and daily cycles, and trends, as well as comparing among observatories as well as between urban and rural sites [1][2][3][4]. This research is determined by these gases' links to climate change. However, both natural and anthropogenic local contributions appear as outliers that may have a marked impact [5,6]. Transport from pollution sources, such as cities, or low dispersion in the atmospheric boundary layer under stable stratification may determine these outliers, whose noticeable influence on the concentration trend is revealed when no robust statistics are used to calculate this trend.
In order to obtain the long-term concentration trend and to analyse the sources and sinks close to the observatory, which are measurements that need to be filtered [7]. Following Zhang et al. [8], statistical filters identify observations that differ from the smooth fit of measurements. The first objective of this paper is to establish the role played by outliers both in the trend and in the annual CO 2 and CH 4 cycle measured at a semi-urban site. Satar et al. [9] used a simple procedure based on the 2-σ filter followed by a 30-day moving average to obtain annual CO 2 and CH 4 evolutions. Other procedures involve two windows to filter fast observations: one for the annual increase rate, which is sensitive to long-term changes, and a second window that responds to seasonal amplitude [10,11], although additional windows may be included in order to detect other trends, such as the interannual component [12].
Following previous studies, this paper distinguishes between the trend of measurements and their annual cycle. However, the current research also focuses on the contrast between regular observations and outliers. Since the autocorrelation function provides information concerning the persistence of a time series [13], a simple form of this function is used to form groups of observations, whose trend is obtained for both regular measurements and outliers.
Once observations are detrended, their annual evolution is investigated. Although the different patterns observed are linked to the measurement site, all of them resemble a harmonic function [14], which has frequently been used to explain how these gases evolve [15]. However, to date, a range of choices has been considered that differ in the number of harmonics used, although the reasons for such choices are frequently not given. The second objective of this paper is to analyse the number of harmonics that should be used to describe the annual cycle in a simple way, while preserving seasonal evolution singularities, such as the number of maxima, the period between them, or the seasonal range.
The final aim of this paper is to provide fresh insights into measurement treatment in order to obtain a satisfactory description of CO 2 and CH 4 evolution and the features of their annual cycle by adopting procedures that simplify the methods traditionally used.

Experimental Description
CO 2 and CH 4 dry concentrations were measured with a Picarro G1301 analyser that uses the cavity ring down spectroscopy technique. In this method, a laser beam describes a ring in a cavity. When the beam is turned off, the quantity and time of the decreasing signal at the detector depends on the gas concentration at the cavity. This analyser is located at the CIBA station, Low Atmosphere Research Centre (41 • 48'50 N, 4 • 55'59 W, 852 m a.s.l., Figure 1a) in the centre of the northern plateau of the Iberian Peninsula. Measurements were provided each 2-3 s. However, hourly means were calculated, and these values were used as initial observations for further calculations, such as monthly percentiles. The measurement campaign was extended over six years from 15 October 2010. Calibrations were made each two weeks with three standard concentrations prepared by the NOAA (National Oceanic and Atmospheric Administration), and measured concentrations were corrected following the calibration equations. Although measurements were considered at three levels (1.8 m, 3.7 m, and 8.3 m, Figure 1b), only the lowest was used in the current analysis, since differences among the observations were very low [16]. Following the Köppen Climate Classification, the climate at the site is Cfb, and it is temperate with a dry season and a temperate summer [17]. Vegetation comprises Mediterranean shrubland that develops mainly in spring and dies or sheds its leaves in summer. However, evergreen shrubs, such as thyme, are also present. Moreover, some rainfed crops are found in the surrounding areas.

Statistics
Monthly percentiles were calculated and fitted to linear equations, which were used to detrend the percentiles. Different statistics were calculated to describe the residuals. These statistics include the location indicator; the median of absolute values, the interquartile range (IQR)-in other words, the difference between the third and the first quartiles-which is used as a spread estimator, and the Yule-Kendall index (YK), which is a symmetry statistic [13], where Q 1 , Q 2 , and Q 3 are the first, second, and third quartiles, respectively.

The lag-1 Autocorrelation as a Randomness Indicator
The concentration evolution is formed by two contributions, one of which is deterministic and the other random. By way of an example, let us suppose the function where t is the time in months and ε is a random variable uniformly distributed with an average null and variable range, R(ε), from 0 to 3, whose values are taken in 0.1 steps. Figure 2a presents Equation (2) with increasing values of R(ε). Small values of this range determine slight differences against the deterministic part of the function. However, large values of R(ε) hide this deterministic contribution. Time series analysis is a frequently used tool to investigate variables that display periodic behaviour. Methods such as Markov chains and autoregressive models have successfully been applied in this research line. The autocorrelation function has also been used in atmospheric research to provide information concerning time series persistence [18]. Recently, Anderson and Gough [19] considered the lag-1 autocorrelation when studying the influence of the previous value on the current one. Moreover, since this correlation coefficient, r, depends on the values of the random variable ε, Equation (2) was calculated 1000 times in the six-year period, and the lag-1 autocorrelation was obtained with each series. Figure 2b represents the quartiles of r for the lag-1 as a function of R(ε). Small values of this range are linked with the "memory" of the deterministic contribution of Equation (2). However, when R(ε) increases, r drops to low values, revealing a "memory loss" due to the increasing weight of the random part of Equation (2).

Statistics
Monthly percentiles were calculated and fitted to linear equations, which were used to detrend the percentiles. Different statistics were calculated to describe the residuals. These statistics include the location indicator; the median of absolute values, the interquartile range (IQR)-in other words, the difference between the third and the first quartiles-which is used as a spread estimator, and the Yule-Kendall index (YK), which is a symmetry statistic [13], where Q1, Q2, and Q3 are the first, second, and third quartiles, respectively.  (2). However, when R(ε) increases, r drops to low values, revealing a "memory loss" due to the increasing weight of the random part of Equation (2). Correlation coefficient for the lag-1 autocorrelation as a function of the range of the random variable ε. One thousand functions were proposed for each range. The median is presented in black, and the interquartile range is shown in grey.

Harmonic Composition
Following [20][21][22], let us consider the function formed by the addition of two harmonics, y = a0 + a1 cos(2π/T1 t − θ1) + a2 cos(2π/T2 t − θ2) where y is the concentration, T1 and T2 are periods established for both harmonics, and t is the time in months. Its coefficients may be fitted by multiple linear regression. In particular, the addition of two harmonics, one with an annual cycle, T1 = 12 months, an amplitude equal to one, a1 = 1 and a phase constant equal to zero, θ1 = 0, and the second with a semi-annual cycle, although with a different amplitude and phase constant, which provides the result presented in Figure 3. The third and fourth harmonics (ai cos(2π/Ti t − θi), where T3 = 4 months for the third harmonic and T4 = 3 months for the fourth harmonic) were also added to Equation (3), and its response against specific values of the parameters is also presented in Figure 3.  (2) with four selected ranges of the random variable ε. (b) Correlation coefficient for the lag-1 autocorrelation as a function of the range of the random variable ε. One thousand functions were proposed for each range. The median is presented in black, and the interquartile range is shown in grey.

Harmonic Composition
Following [20][21][22], let us consider the function formed by the addition of two harmonics, y = a 0 + a 1 cos(2π/T 1 t − θ 1 ) + a 2 cos(2π/T 2 t − θ 2 ) where y is the concentration, T 1 and T 2 are periods established for both harmonics, and t is the time in months. Its coefficients may be fitted by multiple linear regression. In particular, the addition of two harmonics, one with an annual cycle, T 1 = 12 months, an amplitude equal to one, a 1 = 1 and a phase constant equal to zero, θ 1 = 0, and the second with a semi-annual cycle, although with a different amplitude and phase constant, which provides the result presented in Figure 3. The third and fourth harmonics (a i cos(2π/T i t − θ i ), where T 3 = 4 months for the third harmonic and T 4 = 3 months for the fourth harmonic) were also added to Equation (3), and its response against specific values of the parameters is also presented in Figure 3.  Atmosphere 2020, 11, 769 6 of 16

Evolution of Monthly Percentiles
Percentiles were calculated with monthly observations in order to explore their evolution. Figure 4 presents the results for selected percentiles. The period considered is short enough to consider that it may be successfully described by a linear evolution. However, irregular shapes are obtained for extreme percentiles, especially for the highest, where the greatest outliers occasionally determine noticeable concentrations. Although outliers are also present in the low percentiles, their influence is less marked than in the high percentiles.

Evolution of Monthly Percentiles
Percentiles were calculated with monthly observations in order to explore their evolution. Figure  4 presents the results for selected percentiles. The period considered is short enough to consider that it may be successfully described by a linear evolution. However, irregular shapes are obtained for extreme percentiles, especially for the highest, where the greatest outliers occasionally determine noticeable concentrations. Although outliers are also present in the low percentiles, their influence is less marked than in the high percentiles. Linear regressions were calculated for percentiles. Intercepts and slopes are presented in Figure 5. Moreover, limits under which correlations are statistically significant at a 0.001 level are marked. Intercepts increase from 387.06 to 412.82 ppm for CO2 and from 1.8343 and 1.9251 ppm for CH4. However, slopes are steadier, around 2.24 ppm year −1 for CO2 and 0.0097 ppm year −1 for CH4 below the percentiles marked. Linear regressions were calculated for percentiles. Intercepts and slopes are presented in Figure 5. Moreover, limits under which correlations are statistically significant at a 0.001 level are marked. Intercepts increase from 387.06 to 412.82 ppm for CO 2 and from 1.8343 and 1.9251 ppm for CH 4 . However, slopes are steadier, around 2.24 ppm year −1 for CO 2 and 0.0097 ppm year −1 for CH 4 below the percentiles marked.
Monthly percentiles were detrended, and robust statistics of residuals (i.e., differences from the linear regression) introduced above were calculated and presented in Figure 6. Below the 80th percentile, the medians of the absolute value and the interquartile range of the residuals is very similar, around 3 and 6 ppm, respectively, for CO 2 and around 0.01 and 0.02 ppm, respectively, for CH 4 . The Yule-Kendall index of residuals was negative for most percentiles, particularly between the 40th and 70th percentiles for CO 2 . However, this index was negative below the 50th percentile for CH 4 residuals. Monthly percentiles were detrended, and robust statistics of residuals (i.e., differences from the linear regression) introduced above were calculated and presented in Figure 6. Below the 80 th percentile, the medians of the absolute value and the interquartile range of the residuals is very similar, around 3 and 6 ppm, respectively, for CO2 and around 0.01 and 0.02 ppm, respectively, for CH4. The Yule-Kendall index of residuals was negative for most percentiles, particularly between the 40 th and 70 th percentiles for CO2. However, this index was negative below the 50 th percentile for CH4 residuals.  Monthly percentiles were detrended, and robust statistics of residuals (i.e., differences from the linear regression) introduced above were calculated and presented in Figure 6. Below the 80 th percentile, the medians of the absolute value and the interquartile range of the residuals is very similar, around 3 and 6 ppm, respectively, for CO2 and around 0.01 and 0.02 ppm, respectively, for CH4. The Yule-Kendall index of residuals was negative for most percentiles, particularly between the 40 th and 70 th percentiles for CO2. However, this index was negative below the 50 th percentile for CH4 residuals.

Lag-1 Autocorrelation
The lag-1 autocorrelation was calculated for the time series of monthly percentiles and presented in Figure 7. Moreover, ∆ lag-1 for a specific percentile was calculated as the difference between lag-1 Atmosphere 2020, 11, 769 8 of 16 autocorrelation for the following percentile and this value for the preceding percentile. This quantity was used as an indicator of the change in the lag-1 autocorrelation with the percentile. Varied groups of percentiles were suggested following both the lag-1 autocorrelation and ∆ lag-1. Five groups were proposed for CO 2 , where the first corresponds to an increase in the lag-1 autocorrelation. This value is steady in the second group. However, it decreases in the third group and at a greater rate in the fourth, where oscillations of ∆ lag-1 were also noticeable. Finally, the last group was associated with major changes in ∆ lag-1. Six groups of percentiles were suggested for CH 4 , where the first group features fast oscillations of ∆ lag-1. These oscillations were wider in the second group. The lag-1 autocorrelation was steady in the third group and decreased slowly in the fourth. This decrease is more pronounced in the fifth group and is very sharp in the sixth.

Lag-1 Autocorrelation
The lag-1 autocorrelation was calculated for the time series of monthly percentiles and presented in Figure 7. Moreover, Δ lag-1 for a specific percentile was calculated as the difference between lag-1 autocorrelation for the following percentile and this value for the preceding percentile. This quantity was used as an indicator of the change in the lag-1 autocorrelation with the percentile. Varied groups of percentiles were suggested following both the lag-1 autocorrelation and Δ lag-1. Five groups were proposed for CO2, where the first corresponds to an increase in the lag-1 autocorrelation. This value is steady in the second group. However, it decreases in the third group and at a greater rate in the fourth, where oscillations of Δ lag-1 were also noticeable. Finally, the last group was associated with major changes in Δ lag-1. Six groups of percentiles were suggested for CH4, where the first group features fast oscillations of Δ lag-1. These oscillations were wider in the second group. The lag-1 autocorrelation was steady in the third group and decreased slowly in the fourth. This decrease is more pronounced in the fifth group and is very sharp in the sixth.

Annual Cycle
The initial observations were divided into the groups presented in Figure 7 and were detrended. Table 1 presents the parameters of linear fits. Correlation coefficients are high for central groups, but especially low for group 5 of CO2 and group 6 of CH4. Residuals were used to establish the annual cycle. Figure 8 presents the monthly medians and interquartile ranges for each group. The annual cycle is similar for both gases, since a high concentration is observed in winter and low values are reached in summer. However, groups 1 and 2 present a second minimum in May (which is more visible in Figure 9) that is not observed in the remaining groups for CO2. Moreover, the decrease in spring is slow, and the increase in autumn is fast. The behaviour of groups is similar for CH4, since interquartile ranges and skewnesses are nearly the same. However, both statistics depend on the group and month for CO2. The interquartile range increases from group 1 to 5 and, for group 5, skewness is mainly negative in January and February, but it is noticeably positive in May.

Annual Cycle
The initial observations were divided into the groups presented in Figure 7 and were detrended. Table 1 presents the parameters of linear fits. Correlation coefficients are high for central groups, but especially low for group 5 of CO 2 and group 6 of CH 4 . Residuals were used to establish the annual cycle. Figure 8 presents the monthly medians and interquartile ranges for each group. The annual cycle is similar for both gases, since a high concentration is observed in winter and low values are reached in summer. However, groups 1 and 2 present a second minimum in May (which is more visible in Figure 9) that is not observed in the remaining groups for CO 2 . Moreover, the decrease in spring is slow, and the increase in autumn is fast. The behaviour of groups is similar for CH 4 , since interquartile ranges and skewnesses are nearly the same. However, both statistics depend on the group and month for CO 2 . The interquartile range increases from group 1 to 5 and, for group 5, skewness is mainly negative in January and February, but it is noticeably positive in May.   Some combinations of four harmonics have been suggested to fit the monthly medians of the groups formed. A look at the correlation coefficients in Table 2 reveals that the contribution of the fourth harmonic is extremely weak, with the third harmonic contribution being slightly better in general. As a result, the composition of annual and semi-annual cycles provides satisfactory agreement, since the inclusion of higher harmonics complicates the equation unnecessarily. The parameters of Equation (3) are presented in Table 3. The semi-annual amplitude is considerably lower than the annual amplitude for both gases. For CO2, the annual amplitude decreases when the percentile increases, although the semi-annual amplitude increases. Moreover, annual maxima are obtained in January-February, although semi-annual maxima are present around April-May and October-November. For CH4, the semi-annual amplitude decreases when the percentile increases. Annual maxima are obtained around January, although semi-annual maxima are calculated in March-September. The addition of the two harmonics produces the irregular evolution shown in Some combinations of four harmonics have been suggested to fit the monthly medians of the groups formed. A look at the correlation coefficients in Table 2 reveals that the contribution of the fourth harmonic is extremely weak, with the third harmonic contribution being slightly better in general. As a result, the composition of annual and semi-annual cycles provides satisfactory agreement, since the inclusion of higher harmonics complicates the equation unnecessarily. The parameters of Equation (3) are presented in Table 3. The semi-annual amplitude is considerably lower than the annual amplitude for both gases. For CO 2 , the annual amplitude decreases when the percentile increases, although the semi-annual amplitude increases. Moreover, annual maxima are obtained in January-February, although semi-annual maxima are present around April-May and October-November. For CH 4 , the semi-annual amplitude decreases when the percentile increases. Annual maxima are obtained around January, although semi-annual maxima are calculated in March-September. The addition of the two harmonics produces the irregular evolution shown in Figure 9, where the derivative with respect to time is also presented for a better description of the evolution. The concentration of the two gases increases from August to December. However, two types of decrease were observed from January to August, particularly for CO 2 . The decrease is gradual for groups 1 and 2, since the derivative of concentration with respect to time is steady. However, one maximum close to zero is noticeable in March-April, revealing a step in concentration, which is slight for group 3, but more marked in groups 4 and 5. The opposite behaviour emerges for CH 4 , since this step of the concentration is marked for group 1, but gradually decreases from group 2 to group 6.  Table 3. Derivatives of these functions with respect to time are presented in light grey.

Observations and Outliers
Pérez et al. [23] highlighted the different behaviour of selected percentiles and studied decile trends. The analysis presented in the current study expands that research by using all the percentiles. Although most of the measurements follow the evolution given by the trend and the yearly cycle [24], outliers are frequently observed due to local sources or sinks and transport from cities [25]. The shape of the annual cycle and the distribution of outliers are specific to the site. Figure 4 presents an annual cycle with low marked minima in summer for both gases when hemispheric CO2 uptake draws down global CO2 concentrations and CH4 is oxidised by OH throughout the hemisphere, although with noticeable outliers in spring for CO2 attributed to respiration when soils start to warm up but are still moist and vegetation develops, with low dispersion at night when atmospheric stable stratification prevails [26]. Moreover, outliers were frequent in winter for CH4 when the occasional impact of the Valladolid urban plume has also been reported [27]. This plume may be noticeably affected by natural gas heating during winter months and might contribute to CH4 outliers during stable stratification periods. However, hourly concentrations at Sammaltunturi, Finland, described by Lohila et al. [28] showed one very marked minimum in summer for CO2 and frequent high outliers for CH4 throughout the year. Higuchi et al. [29] studied CO2 at Fraserdale, Canada, where a noticeable range in the daily cycle was observed in summer. This pattern was also presented by Zhu and Yoshikawa-Inoue [30] at Rishiri Island in northern Japan. Frequent night-time outliers were attributed to the formation of stable boundary layers in this period together with plant and soil respiration. However, outliers of low CO2 daily means were measured in summer at Takayama, in central Japan [31].  Table 3. Derivatives of these functions with respect to time are presented in light grey.

Observations and Outliers
Pérez et al. [23] highlighted the different behaviour of selected percentiles and studied decile trends. The analysis presented in the current study expands that research by using all the percentiles. Although most of the measurements follow the evolution given by the trend and the yearly cycle [24], outliers are frequently observed due to local sources or sinks and transport from cities [25]. The shape of the annual cycle and the distribution of outliers are specific to the site. Figure 4 presents an annual cycle with low marked minima in summer for both gases when hemispheric CO 2 uptake draws down global CO 2 concentrations and CH 4 is oxidised by OH throughout the hemisphere, although with noticeable outliers in spring for CO 2 attributed to respiration when soils start to warm up but are still moist and vegetation develops, with low dispersion at night when atmospheric stable stratification prevails [26]. Moreover, outliers were frequent in winter for CH 4 when the occasional impact of the Valladolid urban plume has also been reported [27]. This plume may be noticeably affected by natural gas heating during winter months and might contribute to CH 4 outliers during stable stratification periods. However, hourly concentrations at Sammaltunturi, Finland, described by Lohila et al. [28] showed one very marked minimum in summer for CO 2 and frequent high outliers for CH 4 throughout the year. Higuchi et al. [29] studied CO 2 at Fraserdale, Canada, where a noticeable range in the daily cycle was observed in summer. This pattern was also presented by Zhu and Yoshikawa-Inoue [30] at Rishiri Island in northern Japan. Frequent night-time outliers were attributed to the formation of stable boundary layers in this period together with plant and soil respiration. However, outliers of low CO 2 daily means were measured in summer at Takayama, in central Japan [31]. Different behaviour was observed at Minamitorishima, a remote coral island in the western North Pacific, where few noticeable CO 2 outliers were recorded [32].

Annual Trend
Liu et al. [33] considered an equation with a linear evolution for the trend accompanied by two sinusoidal terms and calculated its coefficients for CO 2 in the period 1997-2006 following the vegetation classes from the International Geosphere-Biosphere Programme. They obtained a global value of 2.04 ppm year −1 , which was slightly higher than that recorded at the Mauna Loa Observatory, with the highest trend being 4.92 ppm year −1 for open shrublands against no trend for grassland or mixed forest. The trend value closest to this study site was 2.28 ppm year −1 for the woody savannah, which responds to the steppe feature of the vegetation. The same equation was used by Wu et al. [34] for CO 2 concentration at a tall forest in Northeast China in the period 2003-2010. They calculated 1.7 ppm year −1 , which was close to the value presented by Liu et al. [33] for an evergreen needleleaf forest.
Vermeulen et al. [35] studied daily trimmed means of CO 2 and CH 4 among other gases at Cabaw, in the Netherlands, using an equation with a linear trend and four harmonics, and obtained an increase of 2.00 ppm year −1 for CO 2 in the period 2005-2009 and 0.0074 ppm year −1 for CH 4 in 2005-2010, which is slightly below those recorded in this research. Timokhina et al. [36] used the same equation for CO 2 over central Siberia in the period 2006-2013 and obtained an increase of 2.02 ppm year −1 for CO 2 , which is similar to that of Vermeulen et al. [35]. Bergamaschi et al. [37] recorded a growth rate of around 0.004 ppm year −1 for CH 4 in early 2010 in the Northern Hemisphere. Similarly, Dalsøren et al. [38] reported an increase in global mean surface CH 4 of around 0.006 ppm year −1 in the period 1984-2012.
The low percentile trend is related with large-scale hemispheric and global background concentration evolution, whereas the high percentile trend is related with emission evolution. The results of this research agree with trends recorded at different observatories for low CO 2 percentiles, with the CH 4 trend for low percentiles being higher. However, the CO 2 trend of high percentiles is greater than the low percentile trend, and it reveals the ever-increasing role played by CO 2 emissions. The CH 4 trend is similar for low and high percentiles, indicating that the emission contribution remains steady.

Harmonic Analysis
Equations with only the annual frequency, such as Liu et al. [33] or Wu et al. [34], who used an addition of two sinusoidal functions, respond to the simplest approach for explaining a near-symmetrical annual cycle. However, Figures 8 and 9 reveal an asymmetrical evolution. Several harmonics have often been used to fit observations with enough detail. The most complex analyses consider four harmonics [35,36,39]. However, expressions with three harmonics have also been used [40][41][42]. As a result of this research, two harmonics may be considered as a suitable description of the annual cycle, which is in agreement with previous studies [20][21][22]. In fact, addition of the first and second harmonics has proved to be a suitable procedure to describe cycles with one maximum and one minimum that are unevenly distributed, such as the daily temperature cycle [43]. One result to emerge from the current paper is the extension of this procedure to explain the annual evolution of CO 2 and CH 4 .

Annual Cycle
The contrast among latitudes and the time evolution of the two gases has been described by the NOAA [44]. Pu et al. [45] presented the CO 2 annual cycle at Mauna Loa, Hawaii, which shows a very smooth evolution with one maximum in May and one minimum in September-October. The annual range was slightly higher at Waliguan, China, another remote station, although the maximum was observed in April and the minimum was observed in July-August. This latter evolution was similar to that observed at Jungfraujoch, Switzerland [46].
A similar cycle to that observed in this study was reported by Haszpra et al. [47], whose CO 2 annual cycle at Hegyhátsál, Hungary, reflected a rapid increase from August to December, whereas the decrease in the rest of the year was slower, with a slight step in concentration in spring. However, the range was considerably higher, around 30 ppm. The low annual range at the CIBA station compared to other observatories was also described by Curcoll et al. [48] and could be explained by the weak influence of the sparse vegetation at this site. Hatakka et al. [49] presented the annual CO 2 cycle in Pallas, northern Finland, where the concentration was highest in late winter and could be attributed to ecosystem respiration and the accumulation of this gas in the low atmosphere during this season. Similar cycles were presented by Higuchi et al. [29] for Fraserdale, Canada, with one maximum in April and one minimum in August; Alert, near the North Pole, where the annual cycle range was smaller and the maximum and minimum were delayed by approximately one month, and Cold Bay, Alaska, where the cycle was softer. Xia et al. [50] showed the CO 2 annual cycle at the Lin'an regional background station, some 150 km from Shanghai, China, where a step in spring similar to that presented in Figure 9 for the highest percentiles may be observed. Fang et al. [51] presented the annual CO 2 cycles at the same station, with different filtering procedures. The result was similar to that presented in Figure 9, although with a higher range and where the spring step depends on the filtering procedure followed.
With regard to CH 4 , Zhang et al. [52] compared the seasonal cycle among four observatories. Although two presented similar evolutions to those shown in this research, the seasonal cycle proved to be irregular for Zugspitze/Schneefernerhaus, southern Germany, with maxima in April and August, while the cycle at Waliguan, China, featured one maximum in summer.
Graven et al. [53] studied CO 2 seasonal amplitude change following the latitude in the Northern Hemisphere. They observed that CO 2 seasonal amplitude increases with latitude and with time above 35 • N. However, the shape of the annual evolution may be affected by local features, which is a topic that is currently subject to inquiry [54].

Conclusions
A six-year CO 2 and CH 4 measurement campaign was carried out at the CIBA station, which is a site near a steppe with a slight urban influence located in the centre of the upper Spanish plateau.
The structure of the low atmosphere has a marked impact on CO 2 and CH 4 outliers at this rural site, since they are mainly observed at night when the low atmosphere is stably stratified. High CO 2 values were observed in spring due to ecosystem respiration, whereas CH 4 outliers were observed in winter, when the boundary layer height was low.
Since concentration evolution is slow in the six-year period analysed, linear equations are suitable to describe the trend. Most of the monthly percentiles provided similar annual rates, and residuals are low.
The lag-1 autocorrelation applied to monthly percentiles evidenced the contrast not only between outliers and regular observations, but also between the lowest and the highest concentrations. Five groups of observations were proposed for CO 2 where the highest trend was observed for the highest percentiles, whereas six groups were proposed for CH 4 with the highest trend below the 42nd percentile, i.e., for intermediate and small observations. The annual cycle of residuals was similar for both trace gases, with a marked minimum in summer, and with group behaviour being similar for CH 4 . However, one concentration step in monthly medians was observed not only for high CO 2 percentiles, but also for low CH 4 percentiles. Harmonic analysis applied to these monthly medians revealed that annual and semi-annual harmonics satisfactorily described the seasonal evolution, with additional harmonics proving unnecessary.
Finally, the current research presents and applies a procedure to mark outliers and to examine how they impact average CO 2 and CH 4 concentrations and trends, which may be extended to other ecosystems and environments. Moreover, although the number of harmonics was set at two, further research around the cycle amplitude in the evolution equation is still open to detailed analysis at this site. Funding: This research was funded by The Ministry of Economy and Competitiveness and ERDF funds, project numbers CGL-2009-11979 and CGL2014-53948-P, and by the Regional Government of Castile and Leon, project number VA027G19.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.