key: cord-0073792-5hkhza7f
authors: Deshcherevskii, A. V.; Sidorin, A. Ya.
title: Iterative Algorithm for Time Series Decomposition into Trend and Seasonality: Testing Using the Example of CO(2) Concentrations in the Atmosphere
date: 2022-01-20
journal: Izv
DOI: 10.1134/s0001433821080028
sha: 2df053c1422a809004bfbb7f4dd265699f4a96d7
doc_id: 73792
cord_uid: 5hkhza7f

An iterative algorithm for the decomposition of data series into trend and residual (including the seasonal effect) components is proposed. This algorithm is based on the approaches proposed by the authors in several previous studies and allows unbiased estimates for the trend and seasonal components for data with a strong trend containing different periodic (including seasonal) variations, as well as observational gaps and omissions. The main idea of the algorithm is that both the trend and the seasonal components should be estimated using a signal that is maximally cleaned of any other variations, which are considered a noise. In estimating the trend component, seasonal variation is a noise, and vice versa. The iterative approach allows a priori information to be more completely used in the optimization of models of both trend and seasonal components. The approximation procedure provides maximum flexibility and is fully controllable at all stages of the process. In addition, it allows one to naturally solve the problems in the case of missing observations and defective measurements without filling these dates with artificially simulated values. The algorithm was tested using data on changes in the concentration of CO(2) in the atmosphere at four stations belonging to different latitudinal zones. The choice of these data is explained by the features that complicate the use of other methods, namely, high interannual variability, high-amplitude seasonal variations, and gaps in the series of observed data. This algorithm made it possible to obtain trend estimates (which are of particular importance for studying the characteristics and searching for the causes of global warming) for any time interval, including those that are not multiples of an integer number of years. The rate of increase in the CO(2) content in the atmosphere has also been analyzed. It has been reliably established that in around 2016, the rate of CO(2) accumulation in the atmosphere became stabilized and even tended to decrease. The extent of this stabilization will become clear in the next 2–3 years as new data are available.

Continuous homogeneous observations over natural processes are extremely important, since they allow one to reveal various cyclicities and evolutionary changes (trends) in their development (Vernov, 1990) . Therefore, the methods of geophysical monitoring are more and more widely used. Observational data are processed and analyzed using statistical models that describe experimental data briefly, but with sufficient accuracy. As a result, two main characteristic time scales can often be distinguished: intra-annual and interannual variability.

By seasonal variation S, we mean changes that recur almost similarly from year to year but have some variability. In general, the amplitude, form, and phase of seasonal effects can change, and these changes are studied in different fields of knowledge (Deshcherevskii et al., 1993 (Deshcherevskii et al., , 1997a Gubanov, 2010 Gubanov, , 2017 Dagum and Bianconcini, 2016; Deshcherevskii and Sidorin, 1999; Polonsky and Shokurova, 2010; Timofeeva and Yulin, 2020 ). However, for example, in meteorology, a narrower interpretation of the concept (namely, the climatic norm) is traditionally used. In fact, it implies that the seasonal effects in different years are identical (Rukovodyashchie…, 2017) .

In studies of interannual changes, much attention is given to the identification of long-term unidirectional changes (trends). The concept of a trend likewise allows for different interpretations. In the simplest case, the trend is considered linear, but most often the law of change in the average signal is more complex and more adaptive trend models should be considered. Differences in the interpretation of the concepts (such as seasonal periodicity and trend) make it significantly difficult to obtain reliable results when building statistical models of time series data.

Additional difficulties arise due to the fact that experimental data series normally contain gaps and defects (Gubanov, 2010; Love, 2009; Mandea and Korte, 2011; Deshcherevskii et al., 2016b Deshcherevskii et al., , 2019 Gavrilov et al., 2016; Raffinetti et al., 2007; Rosenberg et al., 1994; Scargle et al., 2013) . Even if there were no gaps in the process of measurements, they will appear when defects are rejected (Deshcherevskii et al., 2016b; Plotnikov et al., 2014) . There are various approaches to processing these signals: for example, filling in the gaps (Kontsevaya, 2012; Shamrik, 2018; Vinogradov and Selyanichev, 2013; Vorob'ev and Vorob'eva, 2018; Plotnikov et al., 2014) . We believe that the approach based on the use of only actually performed measurements is much more correct (Deshcherevskii et al., 2016b (Deshcherevskii et al., , 2020 .

The time series found in different subject areas are very diverse. Therefore, one can hardly offer a unified universal method for building a statistical model of an arbitrary series. In our opinion, such a model should be built with the participation of a researcher who optimizes the model structure and refines its parameters with the help of formalized procedures. However, in doing so, the researcher always relies on some basic (typical) elementary procedures. These include a procedure for identifying and describing the trend. As noted above, the trend can be understood differently. Depending on the properties of the series, the optimal trend model can be built in different ways.

This study aims to describe an iterative algorithm developed by the authors for the decomposition of time series data with gaps into trend and residual (including the seasonal effect) components. The results of using this algorithm for analyzing real data from long-term geophysical observations are also presented.

The algorithm could be reasonably tested by using time series with a clearly expressed and nonlinear trend complicated by periodic variations of the observed parameter. In this study, we focused primarily on sharply nonstationary series with seasonal periodicity. This type includes not only a number of geophysical signals, but also various series of econometric, social, and other similar statistics. Ultimately, we took a series of observations of the concentration of carbon dioxide (CO 2 ) in the atmosphere to test the algorithm. These series have a significant trend with clearly expressed seasonal variations changing from year to year and a sufficient length of observations. In addition, the observations were conducted at a network of stations, and the results obtained at each station had their own specific features. At the same time, these signals have a fairly simple statistical structure, which allows the tested model to be devoid of unnecessary complication.

In this study, we used daily average values of hourly observed changes in the concentration of atmospheric carbon dioxide at four stations of the observation network: Barrow (BPW) and Alaska (1973 Alaska ( -2019 ; Mauna Loa (MLO) and Hawaii (1973 Hawaii ( -2019 ; Tutuila (SMO) and American Samoa ; and South Pole (SPO) and Antarctica (1975 Antarctica ( -2019 (Thoning et al., 2020) .

The calculations were performed using the ABD software package (Deshcherevskii et al., 2016a (Deshcherevskii et al., , 2017a (Deshcherevskii et al., , 2017b . This software package was developed specifically for processing series of long-term monitoring data and takes into account all the specific features of these signals; its architecture provides convenient tools for constructing various algorithms for data analysis and modeling from a variety of basic (elementary) signal processing procedures available in the package (Deshcherevskii et al., 2016b (Deshcherevskii et al., , 2017c .

Initially, the need to remove a trend in analyzing signals arises because it is much easier and more convenient to assume that the signal is stationary during the processing (for stationary signals, there are many highly effective statistical methods with a well-developed theory). However, in practice, many real data contain obvious nonstationary effects (Gazizov, 2016; Gubanov, 2003 Gubanov, , 2012 Osipov, 1992; Nadtoka and Sedov, 1994; Orlov and Osminin, 2011; Bitter et al., 2010; Martyushev et al., 2003; Nguyen et al., 2020) .

In the simplest case, the average level of the series is stabilized by subtracting the linear trend from it (Magritskii and Kenzhebaeva, 2017) . However, when processing real experimental signals, it quickly became clear that this approximation of the trend is often insufficient-the trend differs from its simplest linear model. In view of this, a large number of models have been developed to describe more complex trends. The most common are either analytical approximations (such as by polynomials, exponential, logarithmic, and power functions) or various signal smoothing methods (Gazizov, 2016) . Depending on the available a priori information about the signal and on the filtering goals, a certain model that seemed optimal in the particular case was used.

At present, it is known that the overwhelming majority of natural processes contain flicker noise as one of the main components (Borisov, 2015; Kolody, 2010; Kuznetsov, 2017; Makoviichuk, 2020; Mikhailov, 2010; Zakharov, 2010 Zakharov, , 2013 Boronoev and Trubacheev, 2008; Deshcherevskii and Sidorin, 2003; Deshcherevskii et al., 1997b; Lukk et al., 1996; Reshetnikov et al., 2000; Vstovskii et al., 2005) . Often, the signal is a superposition of quasi-periodic (for example, diurnal, seasonal and/or tidal) components and flicker noise (Deshcherevskii and Sidorin, 1999, 2001) . This means that the remainder from filtered deterministic and quasi-periodic signal components is flicker rather than white noise. Against the background of trend and seasonal fluctuations, the GPS measurement data involve a high-frequency jitter consisting of a combination of white noise and flicker noise (Zakharov, 2013; Prawirodirdj and Bock, 2004; Wang et al., 2012) .

We recall that flicker noise is characterized by a power-law dependence of the spectral power W on the period T: W ~ T β , where the power-law parameter β is between 0.5 and 2.0 (Feder, 1988; Schroeder, 1991; Schuster, 1987) . For flicker noise, nonstationarity is the basic content, the basis of the signal structure. The increase in the power of variations with increasing period means that, among variations with any available observation periods, the variations with the largest possible period always have the largest amplitude (in the statistical sense). If a time interval L is available for observation, the maximum contribution to the observed variations is made by components with a characteristic duration of at least L.

Clearly, in the first (roughest) approximation, these variations can be approximated by a linear function. However, it should be understood that the linear component always describes only a small part of the variations in the flicker-noise component. Removing the linear trend from the flicker noise, we remove only part of the nonstationarity, while the remainder still looks like random walks, but without the longest periods. This can be clearly seen from the graphs of modeled flicker noises with different power-law parameters given in (Deshcherevskii et al., 1997b) .

Thus, the nonstationary component (or trend) for flicker-noise series can be approximated neither by a linear model nor by any functional dependence in general (exponent, polynomial, etc) . This is mainly because the trend is fundamentally random. Even if we formally use some analytical dependence on some interval, it will be violated immediately after going beyond its limits. Therefore, for flicker-noise nonstationarities, the trend cannot be described analytically-it must be identified by nonparametric methods (for example, the moving average).

However, even the analysis of flicker noise processes is still sometimes based on the linear approximation of the trend, although it is clear that this approximation is always very rough and does not completely eliminate the nonstationarity. However, the linear trend model has an important advantage: the number of parameters in the approximating model is a minimum. That is why this approximation may well be useful in some cases.

For example, a linear model can be used with a small amplitude of the trend when the variance of the flicker-noise component is small compared to other signal components and there are no sufficient grounds for choosing a more complex model with a larger number of parameters, since they simply cannot be estimated with sufficient reliability due to the fact that the flicker-noise amplitude is small. A linear model can also be appropriate if the power-law parameter of the flicker noise is small (around 0.5) and the nonstationarity associated with low frequencies has a moderate amplitude compared to relatively higher frequency fluctuations (Deshcherevskii et al., 1997b) .

In addition, the linear model can be used as an initial approximation when the trend amplitude is very high. Identifying and eliminating (using a linear model) the main components of the trend, one can reduce the signal nonstationarity and use for its processing those procedures and methods (including trend construction methods) that could not be applied to the original signal due to its very strong nonstationarity.

In climatology, the trend is often approximated by a linear function. However, this is operational only for sufficiently short series and in cases when the trend amplitude is small. However, in many practically important situations, this approximation is obviously insufficient. For example, its application is completely unsatisfactory in climatology when describing almost any variable on time scales of at least 1000 years. One example is the temperature change in the Northern Hemisphere given in (Mann et al., 1999) , which has become widely known and called a "hockey stick graph" due to its shape.

There are often situations in practice when the trend can be approximated by a linear function over a limited time interval. In this case, the linear approximation can be regarded as the first approximation for describing the quasi-deterministic trend (although this does not exclude the need to additionally eliminate nonstationarities on smaller time scales). However, with an increase in the interval, one has to switch in this case to a piecewise linear approximation of the trend. The question arises about the choice of the number of linear approximation sections, about the determination of their boundaries, etc. Clearly if these boundaries cannot be given from simple a priori considerations (a sharp change in the known external parameters affecting the signal), the nonparametric model becomes preferable in this case as well, since it allows a more reasonable and justified trend model to be constructed.

Another consideration indicating the preference of nonparametric trend models over analytical ones is that analytical models require some a priori grounds allowing one class of functions to be preferred over another. However, a statistical model of the signal is most often built at the initial stage of the observational analysis, when the physical model has not been clear yet. At this stage, the law of change in the average signal level is usually still uncertain; therefore, it is impossible even to determine the class of functions that could be used to search for the most suitable trend model. In this case, it is most logical to estimate the trend as some previously unknown function Q(t), rather than introduce the trend model arbitrarily:

(1)

where D is the original signal (observed series); Q is the trend function, Z is the filtered residue (which may contain, for example, a seasonal component and some other variations), and t is time. Hereafter, the vector quantities denoting respective time series or their fragments within the window (Deshcherevskii et al., 2016b) are boldfaced in the formulas and body of the text.

In view of the above considerations, trend Q is meant to be sufficiently slow (smooth) changes characterizing relatively long-term changes (trends) of the series. Apparently, in general, it is this formulation that is most acceptable, can be a compromise, and most closely matches all the meanings that are embedded in the concept of "trend."

In that case, the main question in the trend estimation (simulation) is how to determine the optimal degree of trend variability? Where is the exact boundary-what variations are included in trend Q and what are included in the residue Z? The most natural approach is to specify the boundary in terms of the characteristic time of variation or its frequency limit.

If only variations with an infinitely large characteristic time (the variation periods exceed the length of the series) are included in the trend, the linear function Q = A L t, where A L is the estimated coefficient (the linear trend parameter) and t is time, can be used as a trend model. Another typical approach is when variations with a characteristic time exceeding the duration of those processes that are of interest currently are attributed to the trend. For example, the long-term slow drift of instrument zero can be considered a trend when studying terrestrial diurnal and semidiurnal tides. In analyzing the seasonal effects, the trend can involve variations with a characteristic time of at least several years.

If we assume that the critical attribute in the trend identification is its frequency composition, the simplest way of the trend identification (construction) will be the frequency filtering of the signal. In this case, the digital signal filtering looks like a convolution with a filter function. Filters of this kind are well known. In addition, there is a well-developed theory for the construction of filters of this kind, optimized for different tasks (Hemming, 1983; Lyons, 2009; Oppenheim and Schafer, 2009; Rabiner and Gold, 1975) .

Usually, the frequency filtering for filter construction predominantly involves the optimization of amplitude-frequency response parameters such as linearity and filter cutoff slope; sometimes the filter sequence

length is additionally minimized. In this case, the filtering sequences normally have alternating signs.

The use of alternating filtering sequences leads to big problems in processing signals with a significant number of gaps. Strictly speaking, if the time step between samplings is not equidistant (this happens if some observations are missed), one should recalculate the filter coefficients for each set of gaps again, i.e., actually at each filtration step (Deshcherevskii et al., 2016b) . This is not only very expensive computationally, but can also lead to various instabilities. If one tries simply to disregard the gaps rather than to recalculate the filter coefficients (calculate the convolution by nongaps with a corresponding renormalization), the filtered signal may contain unpredictable artifacts under unfavorable conditions.

We believe (which is confirmed by the practice of data processing) that the optimal compromise in processing experimental signals with gaps is to use only positive-definite filter sequences. For example, a Gaussian (nuclear sliding smoothing with a Gaussian kernel) has proven itself very well as a filtering sequence (Lyubushin, 1993 (Lyubushin, , 1996 (Lyubushin, , 1998 (Lyubushin, , 2003 (Lyubushin, , 2007 Deshcherevskii et al., 2016b) . The filter cutoff slope decreases in this case, but this is usually noncritical for further analysis of the signal decomposition into a trend and a residual.

The choice of positive-definite filtering sequences makes it possible to use the following compromise approach when gaps are present (Deshcherevskii et al., 2016b) : the filtering sequence is specified once and is not recalculated at each step; further, a criterion of the permissible percentage of gaps is introduced to prevent unacceptable distortions of the filtered signal. As long as the number of gaps does not exceed the acceptable level, the convolution is calculated in the standard way, except that the missing observations are excluded from the summation and the filter weight coefficients used in the convolution are renormalized so that their sum be unity. If the number of gaps (more precisely, the total weight of filter coefficients for missed observations) exceeds the permissible level, the result of the operation is declared a gap. Obviously, this condition can be checked by calculating the filter renormalization coefficient and comparing it with the given critical value.

A detailed description of this approach can be found in (Deshcherevskii et al., 2016b) . The history of processing of various experimental signals shows that this approach, on the one hand, is very simple and can be efficiently implemented, and, on the other, it provides maximum flexibility in processing signals with gaps and can be used to reliably solve problems associated with the filtering of such signals. It is this approach that we use in this study.

In constructing trend function Q (regardless of whether Q is specified by a nonparametric model or analytically), all other signal components should be considered interference. First and foremost, this certainly applies to significant components with high amplitudes, simply because these noises are stronger. Although this fact is quite obvious, we illustrate it using a couple of examples.

Let us consider a sinusoidal signal that simulates the seasonal fluctuation S (Fig. 1a) . Despite the fact that the model of the series is represented by a pure sinusoid (without any trend!), the linear approximation shows a trend for this signal with an amplitude A L = -0.47 relative units per year, which corresponds to a quarter of the amplitude of the seasonal function! Clearly, this trend is "detected" due to the fact that the values of the seasonal function S (sinusoid) in the summer period are lower than in the winter, and the beginning of the interval coincides with the beginning of the "winter" half-period. Although Fig. 1a shows a linear trend, it is clear that the effect will be similar for any other trend model; i.e., the estimates will reveal a nonzero trend, although the model signal is given by a pure sinusoid.

In the following example, the window for estimating the trend parameters is chosen so that the seasonal variation within the window is fully symmetric (Fig. 1b) . Obviously, this choice makes it possible to guarantee that there is no bias in coefficient A L (see Fig. 1a ). However, in this example, we have excluded a number of observations from the signal, which simulates data gaps. As a result, the linear approximation shows a trend for this signal with an amplitude A L = +0.22 relative units per year, which corresponds to 10% of the seasonal function amplitude (see Fig. 1b ).

Clearly, the effect shown in Fig. 1a can always be eliminated by choosing a time interval for trend estimation with starting and ending points when the phase is π/2 or 3π/2. However, this requires some observations at the beginning and at the end of the series to be ignored. If the seasonal function is not fully symmetric (sinusoidal), a special analysis of the seasonal function is required to accurately calculate the special initial phase at which the trend is estimated without a bias. This means that in any case we will not be able to correctly estimate the trend without first calculating the seasonal function. If, in addition, the series contains gaps in observations, the correct estimation of even a linear trend becomes problematic. It should be noted that the situation shown in Fig. 1b is not yet the worst: the interval of gaps is sufficiently far from the right edge of the series. Nevertheless, as mentioned above, the magnitude of the coefficient bias reaches 10% of the amplitude of the seasonal function per year.

Of course, as the length of the series grows, the influence of both effects will decrease. However, with a moderate length of observation series, the influence still remains significant. Indeed, with a series length of 20 years, when the beginning of the signal coincides with the zero phase of the sinusoid (similar to Fig. 1a) , the calculated trend amplitude is A L = 0.0027 relative units per year, or 0.05 relative units per 20 years. For a typical series of air temperatures, the amplitude of the seasonal effect is around 30°, which corresponds to a bias in the estimate of the trend amplitude by 0.8° per 20 years. Since the rate of global warming is estimated at 0.2-0.4° per 20 years, the error in estimating the trend amplitude due to the use of a signal with unfiltered seasonal variation can be twice as large as the measured value.

To exclude the influence of the effects described above, the trend estimation in climatology is customarily based on deviations from the climatic norm or data for only 1 month in different years. This allows one to obtain correct estimates. However, if the series is nonstationary in terms of both the mean and seasonal variation (which is indeed the case, for example, for СО 2 series), this approach seems to be nonoptimal. First, the concept of the climatic norm in this case can change from era to era. Second, the trends estimated for different months may not coincide. As a result, several different trend variants are identified for the same signal. However, the signal decomposition into a single trend function Q(t) and a time-varying seasonal variation looks more natural in this situation (for example, see the M3-M5 models described in (Deshcherevskii et al., 2017c) ).

Another possible approach is to estimate the Q(t) trend from average annual values of the series. However, this estimation has several disadvantages. First, the time sampling in this case is rather coarse. Second, if the beginning and end of observations do not coincide with the boundaries of the calendar year, the "extra" data disappear (they are not used in the trend estimation). Third, the presence of gaps (see Fig. 1b ) makes it impossible to estimate the average annual value; with a formal averaging of all observations made during the calendar year, that value will be shifted and filling in the gaps has the drawbacks described in (Deshcherevskii et al., 2016b) .

In general, it is important that no algorithm for trend construction can ignore the presence of strong seasonality. For example, a moving average with a window multiple of an integer number of years "skips" a part of the seasonal variations into a trend function if the series contains observation gaps. A trend distortion due to seasonality also occurs in the case of the interannual variability (inequality) of seasonal effects.

The most direct way to completely exclude the influence of all the effects mentioned above is to filter in advance the seasonal component S(t) (built by the previously formulated seasonality model) from the series.

In both examples presented above (see Figs. 1a, 1b), the signal is represented by a constant. Regardless of the calendar interval chosen for the estimation and/or the presence of data omissions and observational gaps, any reasonable estimation of Q(t) in this case yields a correct result, i.e., shows the complete absence of trend components.

Summing up, we come to the conclusion that, during the estimation of Q(t), all other components of the signal should be considered a noise. In geophysical signals, the seasonal component S(t) usually has a very high amplitude. It is this component that appears as the most significant noise during trend estimation. Very strong distortions arise if the series simultaneously contains both a seasonal effect and data omissions (observational gaps).

In this case, one cannot correctly estimate the trend function without knowledge of the properties of seasonal variation. If the function S(t) is known, the most natural and direct way to estimate the trend is to first filter out function S(t) from the signal and then estimate Q(t) by the filtered residual.

The trend estimation using the sliding smoothing method usually assumes that the smoothing window should be completely located within the series. In this case, the total run of the window is equal to the series length minus the window width. This leads to a decrease in the length of the filtered signal by the window width, which is undesirable. Deshcherevskii et al. (2016b) proposed an approach that allows one to avoid the signal length reduction with this filtering. If gaps are allowed in the data model, the series is simply supplemented with gaps on the right and left sides, and then algorithms are applied adapted to work with a signal containing observational gaps. In this case, the sliding window completely runs through the entire series; however, at the initial time, the left border of the window is located left of the observation start time, and, at the final time, it is located right of the observation end time.

This approach works well in most situations, but it leads to a distortion of the identified nonparametric trend if there is very strong nonstationarity of the mean (Fig. 1c) . The reason is clear: in the absence of information about values outside the series (when they are specified by gaps), the averaging is performed only over data that are within the window, which breaks the balance between the right and left parts of the window. In fact, this algorithm is equivalent to the implicit addition of a series with a constant equal to the boundary value of the series. For other weighting functions of the window (for example, with frequency filtering), the distortion type will change, but the essence of the effect-the trend bias at the series edges-will remain the same. Obviously, the very same edge effects will be observed within the signal at the boundaries of very long series of gaps.

To exclude the bias described here, one should ensure the quasi-stationarity of the series near its boundaries, as well as in the segments with long breaks in observations. This can be done, for example, by subtracting the linear component of the trend from the signal, estimating the nonparametric trend function by the method of moving averaging or frequency filtering, and then constructing the full trend by a superposition of the linear and nonparametric components of the trend.

Thus, when estimating the nonparametric trend using the window collapse algorithm (the series is implicitly supplemented with gaps to ensure that the sliding window goes beyond the series boundaries (Deshcherevskii et al., 2016b) ), the very strong nonstationarity of the signal should also be considered a noise. In this case, one can identify and subtract the linear component of the trend before estimating the nonparametric trend (moving average, frequency filtering, etc.). Then, a cumulative trend is formed as the sum of the linear component and the nonparametric trend.

For estimating the seasonal function S(t), all the remaining signal components should be considered noise. Deshcherevskii et al. (2019) give an example of catastrophic distortion of S(t) when it is estimated by the method of epoch superposition from a signal with level shifts and observational gaps.

The situation is completely analogous when estimating S(t)) from a signal with a strong trend Q(t) if the length of the series is not a multiple of an integer number of years. Figure 2 shows an example of a model (synthetic) series derived from a sinusoid and a linear trend supplemented to it (curve 1). Graph 2 shows the seasonal mean function S(t) estimated from this series by the method of epoch superposition. It can be seen that it differs from the sinusoid by the presence of the two jumps that occurred on June 1 and November 1. The times of these jumps coincide with the times of the beginning and end of the series. Clearly, if there are gaps, similar jumps will be observed on times that coincide with the times of the beginning and end of the series of gaps.

It should be emphasized that this effect is associated precisely with the presence of the trend. At zero trend amplitude Q(t) ≡ 0, the introduction of any number of gaps is noncritical: S(t) is identical with the model sinusoid. It should also be noted that if the signal contains a significant trend component, any estimate of the seasonal function S(t) obtained by any models under the assumption that the nonseasonal The above description means that the signal should be cleared from the trend before the estimation of the seasonal mean function (SMF). Otherwise, the seasonal variation will be estimated with large distortions. On the calculated SMF, they appear as jumps (level shifts) on calendar dates that coincide with the beginning and end of the series, as well as with the beginning and end of observational gaps.

Thus, the trend should be estimated using the signal with a filtered seasonal component, the estimation of which requires the trend be first estimated and filtered. In that case, how can a season and trend model be built? After all, both are unknown to us.

In this case, the common approach is to simultaneously estimate the parameters of Q(t) and S(t). To this end, the multiparametric model W(t) = Q(t) and S(t) is built, which includes both components at once. The & symbol indicates that the combined model W(t) is not necessarily built as the sum of Q(t) and S(t), but may include a more complex combination of these functions. For example, in the case of a multiplicative model (Deshcherevskii and Sidorin, 1998, 2000) , this can be a product. Then, the parameters of the model W(t) are somehow estimated, chosen so as to minimize the mismatch between W(t) and the observed signal D(t).

However, this approach has several obvious disadvantages. First, the estimation of the parameters of the combined model W(t) requires the class of functions that simulate both the trend and residue to be determined in advance, since otherwise the minimization problem cannot be formalized. However, the choice of a class of functions that contain Q(t) and S(t) is not a formal procedure: it always requires the participation of an expert. The best way to make a reasonable choice is to analyze the required component of the signal, clearing it from other components as much as possible. The researcher not only analyzes the behavior of the cleared signal, but also takes into account all available a priori information that cannot be formalized. For example, when choosing the type of S(t), an expert must decide whether the seasonality model should allow variability of the amplitude and phase of seasonal effects, whether the seasonal variation form can change from year to year, whether it should be approximated by some simple function (for example, a sinusoid), or if a different estimate should be used.

Similarly, when building a model of the trend Q(t), it is the expert who decides on the model type and the degree of its adaptability. It is hardly possible to formalize this choice. For example, the trend can be specified by a linear function, but sometimes (when additional a priori information is available and the signal type is taken into account) a piecewise linear approximation can be used, when several fragments of the signal are selected and the trend is approximated separately within each fragment. In other cases, the nonparametric trend estimation may be favorable. In this case, when describing the trend Q(t), an expert must not only choose a specific trend model (the type of nonparametric estimation), but also specify the degree of its adaptability. For example, in frequency filtering, the adaptability is regulated by the choice of the filter cutoff frequency, and when the trend is estimated by moving smoothing, the adaptability is regulated by the width of the smoothing window and the type of the smoothing kernel.

When using the complex model W(t), the trend Q(t) and the seasonal effect S(t) are estimated simultaneously. In this case, it is much more difficult for a researcher than in component-wise analysis to determine the degree of optimality of Q(t) and S(t). In this respect, the complex (combined) model W(t) looks like a "black box": the expert specifies the initial settings and immediately obtains the result at the output, instead of completely controlling all the details. Certainly, depending on the goals of the analysis, this can be not only a disadvantage, but also an advantage. However, for data-researching activities, full control and expert involvement in all stages of the process is clearly preferable.

The second disadvantage of the approach, which is characterized by a simultaneous estimation of the parameters Q(t) and S(t), is the lower stability of the multiparametric estimation. This problem is exacerbated if there are gaps in the data. In this case, it is almost impossible to analytically analyze the properties of the minimization algorithm and prove its stability and convergence. The numerical simulation and expert analysis of the suitability of the algorithm turn out to be very laborious because the situation is multifactorial.

However, if Q(t) and S(t) are estimated separately, the problem is highly simplified. The number of variants that need to be considered in numerical modeling (testing) is reduced, and the assessment of filtering results (the quality of the approximation of the corresponding component) during expert analysis is also facilitated. This is very important since it is the expert who chooses the optimal variant of the model in these studies. It is almost impossible and hardly necessary to formalize this choice, since the optimization criteria are not clear in advance. In addition, these criteria may differ for various models, since the filtering goals are often specified when the optimal model is chosen-this is a quite normal for an exploratory (research) analysis of observations. Thus, the approach based on the simultaneous estimation of parameters of the trend Q(t) and the sea-sonal function S(t) has a number of disadvantages. However, let us find an alternative regardless of the fact that the signal simultaneously contains both a significant seasonal variation and a high-amplitude trend. In addition, in this case, they cannot be described analytically. On the contrary, optimal models of Q(t) and S(t) are searched during the analysis of observations; i.e., the problem of analysis is to choose the model type for each signal component, rather than to estimate the parameters of a preliminary formulated model of W(t). It has been shown above that, in the presence of a high-amplitude trend in the signal, the seasonal function S(t) is estimated with errors, and, conversely, in the presence of significant seasonal variations, the estimated trend Q(t) will be biased.

However, the problem has a solution.

First, it should be noted that the distortions (biases) of estimated parameters of the functions Q(t) and S(t) decrease under the influence of various noises appropriately with the decrease in the amplitude of these noises. Thus, the estimation errors for Q(t) and S(t) can be reduced without completely clearing the signal from all other variations. It is quite enough to decrease their amplitude to an acceptable level.

Second, even if significant noises are available, the bias of estimates in many cases can be reduced by increasing the signal length (estimation window). Let us consider this in more detail.

It can be seen from the examples given in Fig. 1 that the estimate of Q(t) will be biased in estimating the linear trend by a series containing a strong seasonal effect. Recall that the initial signal in Fig. 1 does not contain a trend at all. However, it is quite obvious that this bias decreases rather quickly as the number of seasonality periods increases. In the first approximation, the calculation error for Q(t) is comparable with o(А S /N), where А S is the amplitude of seasonal variation and N is the number of seasonality periods.

When the nonparametric trend Q(t) is estimated using running window algorithms, the influence of seasonality on the estimates will also be smaller the more seasonality periods fall into the window. Therefore, if the amplitude of seasonal effects is high, the trend distortions can be reduced by using the widest possible trend estimation window. However, if the main portion of the seasonal variation is filtered out, a narrower window can be used. Even rather coarse models of S(t) normally allow the amplitude of seasonal effects in the filtered series to be repeatedly suppressed. After filtering, the estimation errors of Q(t) can be reduced to quite acceptable values even if a window with a width of several years is used. Since the trend Q(t) usually includes only longer period variations, this is quite enough in most cases.

When the seasonal function S(t) is estimated from a signal with a strong trend Q(t), the function S(t) can be very strongly distorted (see Fig. 2 ). The amplitude I of these distortions is greater the shorter the series is (the smaller is the number of seasonality periods) and the higher the trend amplitude: I ∼ o(А Q /N), where А Q is the trend amplitude and N is the number of seasonality periods.

If the main component of the trend is filtered out (for example, using a linear approximation of Q(t)), the distortions of S(t) usually become small at a width of the estimation window of S(t) of around 10 years. In the case of a more accurate trend filtering, smaller windows can also be used.

To construct an iterative algorithm with guaranteed convergence, it is sufficient that the bias of the estimates of Q(t) and S(t) are smaller at each step than these functions themselves. To do this, the iterative process must begin with the estimation and filtering of the component that has the maximum variance.

Thus, the main idea of the iterative algorithm for estimating the functions Q(t) and S(t) is to estimate each of these functions using a signal that is maximally cleared of any other variations that are considered a noise. In this case, the optimal model for Q(t) and S(t) is chosen by an expert assessment, since this procedure cannot and should not be formalized. Clearly, the quality of the expert assessment strongly depends on the competence of the researcher, but it is equally important that all the necessary information about the structure of the signal component is presented in the most convenient form, with a minimum influence of noise. This is exactly what is ensured when the iterative algorithm is used.

ALGORITHM FOR ESTIMATING THE СО 2 CONCENTRATION TREND 5.1. Original Data Using the mechanism described in the previous section, we processed data on variations in the atmospheric CO 2 concentration at all stations chosen for this study: Barrow (BRW), Mauna Loa (MLO), Tutuila (SMO), and South Pole (SPO). The initial series of observations are shown in Fig. 3 . It can be seen that all the series have a high-amplitude trend and a seasonal component, while the graphs are sufficiently close in terms of their general course, although the amplitude of the seasonal variation is highly diverse at different latitudes. The latter is due to the fact that the atmospheric mixing along the meridian is rather slow and the main sources of CO 2 flow into the atmosphere are distributed very nonuniformly by the latitude. The effect due to the vegetation variation is almost entirely tropic-related, and anthropogenic emissions are also concentrated in a relatively small latitudinal zone.

The data series are around 45 years in length. All series have a significant number of gaps in observations (Table 1) 

The СО 2 concentration trend at all stations is close to linear (Table 2) . However, an analysis of the graphs of series with a filtered linear trend, or (1) LT (Fig. 4) , shows that the trend characteristics clearly change over time. The left index (1) indicates that the estimate refers to the first iteration step.

This becomes even more obvious if we filter out both the linear trend and the seasonal variation estimated by the method of epoch superposition (estimate (1) S). Figure 5 shows the CO 2 concentrations filtered from the linear trend and seasonal effects at the stations listed above. It can be seen that the growth rate of the CO 2 concentration clearly varies with time. The strongest break (increase in the rate) occurred in 1995-1998, and the second jump in the growth rate (with a smaller amplitude) can be fixed nearly to 2010. Clearly, the linear approximation of the trend is too rough, since it does not reflect all these features of the curves, which are rather clearly recorded in the expert analysis.

The discussion in the previous section clearly shows that the linear approximation of the trend is unsatisfactory. In addition, this approximation fails to trace the trend variations-this requires a more adaptive model. Here, the specification of such a model by any analytical function seems doubtful both due to the absence of an obvious correspondence between the observed variations of any known function and due to the absence of a priori physical considerations for its choice. Thus, it makes sense to estimate the trend component of CO 2 variations nonparametrically.

It has been noted above (see Fig. 1 ) that the signal nonstationarity can lead to a bias in the estimates of the nonparametric trend (NPT). These distortions occur at the ends of the series and at the boundaries of large intervals of gaps. In the case of CO 2 series, the NPT estimation window should be at least several years. Therefore, observational gaps with a length of less than 1 year (which are present in the data of all four stations) can be considered small. Thus, it is necessary to prevent only the distortions that may occur at the ends of the series.

To this end, we built a refined piecewise linear approximation of the trend (PLT), flattening the nonstationarity of the mean value at the beginning and end of each series. It seems that this goal can be achieved by merely splitting the series into two nearly equal fragments (22-23 years in length). In this case, the nonstationarity of the residual signal D-PLT at the boundaries of the series will be smaller than the amplitude of its high-frequency fluctuations and, therefore, the PLT bias will also be small (see Fig. 5 ).

To construct the PLT, the following calculations were performed for each series:

(1) the estimate of the seasonal component (1) S was subtracted from the original signal (0) D:

(2)

(2) the characteristics of the linear trend were estimated at intervals from the beginning of the series to December 31, 1997, and from January 1, 2000, to the end of the series;

(3) the resulting linear functions were extrapolated to 1998-1999 up to the point of their intersection.

The piecewise linear function obtained by the superposition of two fragments specifies a refined (in comparison with the linear function) trend model, or (2) PLT. Analysis shows that the difference between (2) PLT and measurements in most cases does not exceed ±5 ppm for the series of BRW station, ±3 ppm for MLO and SMO stations, and ±2 ppm for SPO station.

At the same time, the analysis of difference series 

(2) PLT(t) (Fig. 6 ) made it possible to detect several outliers that strongly break the normal course of the curves. First and foremost, for the series of BRW station, this is the value of 324.98 on May 13, 1976, which is less than the neighboring values by 10 ppm (the difference exceeds 10σ, where σ is the normal amplitude of concentration fluctuations). In view of

the very large declared measurement error for this date, this value was rejected, and the further calculations used an observation gap instead.

In addition, the BRW station was characterized by three more outliers: 381.54, 380.37, and 390.98 observed on July 19, 2013, July 20, 2013, and July 12, 2016, respectively. They all also differ from the normal local curve by almost 10 ppm, which is more than 10σ. Although the declared measurement errors in this case are much lower (corresponding to the common values), all the indicated values are surrounded by series of gaps before and after the given day. For this reason, they were also rejected before further analysis.

At MLO and SMO stations, no abnormal values deviating from the normal course of the curve by more than 10σ were found. However, the check for SPO station with the same criterion revealed a whole range of abnormal behavior of the CO 2 concentration from November 25, 1976 , to February 27, 1977 . During these months, the CO 2 concentration change at SPO station was in antiphase with the common seasonal rhythm; at the beginning and end of the fragment, the series underwent sharp breaks (jumps) with a strong level shift. For this reason, the specified fragment was also rejected. Now, we consider the behavior of the filtered residue (2) R(t) in 1998-1999. Since the CO 2 concentration growth rate changes gradually rather than abruptly and the piecewise-linear trend has a jump in the derivative at that time, this method of trend construction leads to a rather significant discrepancy between observations Years ΔCO 2 , ppm and (2) PLT, which can be clearly seen from the series (2) R(t) (see Fig. 6 ). In this regard, it should be emphasized that the construction of (2) PLT is a purely auxiliary technique, a preliminary step before evaluating the nonparametric trend. All discrepancies between (2) PLT and the experimental signal are automatically included in the NPT. The real trend is built as the sum PLT + NPT. Therefore, it makes no sense to analyze them separately.

In addition, it can be seen from the CO 2 concentration variations shown in Fig. 6 that the deviation of (2) PLT from the observed curve in 1998-1999 does not exceed the similar deviations in amplitude in other time periods. Therefore, it can be expected that the discrepancy will not lead to an additional bias in estimates for the total trend and there is no need for (2) PLT to be constructed from a larger number of fragments (three, rather than two).

The idea of the iterative algorithm implies that the nonparametric trend should be estimated accordingly with respect to the series cleared of seasonal effects, the main factors of nonstationarity, and all signal defects identified at previous steps. Thus, the estimation of (4) NPT included the following preliminary steps:

(1) a piecewise linear trend built at the previous iteration step was subtracted from the original signal. The resulting difference series (3) R(t) = (0) D -(2) PLT are shown in Fig. 7 . It should be emphasized that the series (3) R(t) contain a seasonal effect, which will be subtracted at the next step using a refined model of seasonal variation;

(2) difference signals (3) R(t) are used to estimate the seasonal component (3) S. However, analyzing the series in Fig. 7 shows that the amplitude of seasonal effects has some tendency to increase over time and the differences between the beginning and end of the observation period are not so large to be taken into account in the construction of the nonparametric trend with a smoothing window of a width of several years;

(3) the constructed seasonal functions (3) S were subtracted from the difference signal. The residuals, i.e., the difference series (4) R(t) = (0) D -(2) PLT -(3) S, are shown in Fig. 8 .

Operations (1)-(3) results in CO 2 variations cleared from the piecewise-linear trend and seasonal effects. It is these series that are used to estimate the nonparametric trend. However, before turning directly to the estimation of (4)NPT, we consider the resulting series to formulate the optimal trend model. To this end, we analyze the graphs of the series (4) R(t) (see Fig. 8 ) and their spectra (Fig. 9) .

Analyses of the series (4) R(t) shown in Fig. 8 makes it possible to identify fluctuations with a characteristic time of 3-5 years. The special role of fluctuations with this time period is confirmed by the analysis of the spectra (see Fig. 9 ). No other features can be seen on the spectra for periods exceeding 1 year. (Clearly, the decay of the spectra over periods exceeding 10 years is explained by the fact that the contribution of corresponding harmonics was removed from signal D when filtering the piecewise-linear trend.) Interestingly, for periods of less than 10 years, all spectra in a double logarithmic (log log) scale seem to be almost perfectly linearized. It should be noted that this nature of the spectra is typical for geophysical time series with a filtered seasonal effect. Table 3 gives the values of the power-law parameter β estimated over the range of periods of 2-4000 days (Deshcherevskii and Zhuravlev, 1996) . It can be seen from the data presented in this table that β is in the range between 0.7 and 1.1, which is typical for geophysical time series (Lukk et al., 1996) . Also, some spectra have singularities in the nearannual period. Obviously, their presence is associated with incomplete filtering of seasonal effects when the model (3) S is used. Generally speaking, a more accurate estimation of seasonal variations requires more adaptive models of seasonality (Deshcherevskii et al., 2017c) . However, it will be shown below that there is no need to refine the seasonality model, since the errors arising due to incomplete filtering of seasonal effects at this step are sufficiently small. A more correct estimate of these errors is given in Section 5.5.

Thus, variations with periods of 3-4 years are of primary interest in refining the trend model. In the remaining range of periods, the series are self-similar (have no identified scales); therefore, the cutoff frequency for attributing variations to a trend or to the residual can be chosen from a priori considerations.

The attribution of variations with periods of 3-4 years to the trend or residual component of the decomposed signal depends on the decomposition goals. In this study, we attributed these variations to the trend component of the signal, bearing in mind that, if necessary, this component can be studied in more detail, including by the additional decomposition of the trend into components.

In view of this, we estimated the nonparametric trend (4) NPT by the method of nuclear sliding smoothing with a window of a width of 4 years and a Gaussian weighting function of the window. An example of the resulting trend is shown in Fig. 10a . It should be noted, that despite the nonparametric estimation, the series (4) NPT are free from the edge distortions shown in Fig. 1c . This result is achieved due to the preliminary removal of the piecewise linear trend from the signal. Figure 10b compares the series (4) NPT for different stations. It can be seen that they are all sufficiently similar (Table 4 ); some differences can be found only at the BRW station. For all other stations, the trends correlate between 0.94 and 0.96. This is a sufficiently high indicator, since the main nonstationarity has already been removed from the signals and variations with a characteristic time of 3 to 10 years are compared. The formal 95% significance level for the reduced correlation coefficients is 0.55-0.58 (the estimate was made from the number of degrees of freedom of the signal).

Trend Accuracy Now, we estimate the error arising from the construction of the nonparametric trends (4) NPT.

The largest contribution to the error is obviously made by high-frequency fluctuations of the signal (4) R(t). The results of numerical simulations show that the nuclear smoothing with a window of 1461 days and a Gaussian kernel allows the high-frequency noise to be suppressed by nearly 30 times. This means that the standard deviation of the resulting trend will be close to 0.03 when the white noise is smoothed with a single standard deviation. The variance respectively decreases by 3 orders of magnitude; however, for practical purposes, it is more convenient to operate with standard deviations, since they are expressed in the same units as the signal under consideration.

The standard deviation of the series (4) R(t) is given in Table 5 (column σ(res)). Since the amplitude (standard deviation) of high-frequency noises for the signals under consideration is 1.4-0.2 ppm, the error in estimating the series (4) NPT is from 0.05 to 0.01 ppm (column σ(npt)) for different stations. Now, we return to the analysis of Fig. 7 . It can be seen that the amplitude of seasonal effects has a definite tendency to increase over time. Although the differences between the beginning and end of the time interval are not so great, the seasonality model (3) S ignores these differences, which can lead to the "infiltration" of seasonal effects into the signal used to estimate the nonparametric trend. The singularities in the near-annual period in the spectra of these series (see Fig. 9 ) confirm that this infiltration does indeed occur. Let us estimate the magnitude of the error that can arise because of this when the nonparametric trend (4) NPT is estimated.

First of all, we estimate the amplitude coefficient of the filter used for the annual period (coefficient P in Table 5 ). This coefficient shows the attenuation of variations with the annual period when sliding smoothing is used. Normally, this coefficient is calcu- Table 3 . Estimates of the power-law flicker-noise parameter β for series (4) R(t) over periods of 2-4000 days. The regression of the spectral power logarithm to the period logarithm was estimated lated for a harmonic signal (the standard amplitudefrequency response of the filter is constructed in this way). However, we estimated it using the numerical simulation method, supplying the seasonality function

S constructed for the given series to the filter input. In this case, the coefficient of periodicity suppression is even higher than in the case of a sinusoid (see Table 5 ). Obviously, this is associated with the complex form of seasonal variation, which, when expanded in a Fourier series, contains not only the annual periodicity, but also multiples of higher frequency periodicities, which are suppressed by the smoothing filter more strongly. Now, we note that the trend estimation is hindered by only the part of the seasonal variation that was not involved in the model component (3) S, rather than by the seasonal variation itself. We estimated its amplitude by analyzing the spectra of series (3) S and (4) R(t). The ratio of the amplitudes of respective spectra for the annual period was calculated (K S in Table 5 ). It can be seen from Table 5 that this ratio varies from around 3 to almost 30. It is very obvious that the stronger the seasonal variation in the signal is, the larger the effect of its filtering and the higher the K S ratio.

Knowing the amplitude of seasonal function (3) S, coefficient P of suppression of seasonal effects in signal (4) R(t), and the coefficient of suppression of the annual periodicity during nuclear smoothing, one can estimate the amplitude of trend fluctuations arising due to the fact that the seasonal effects in (4) R(t) are incompletely suppressed. The values of this amplitude are given in the last column of Table 5 (column σ(Sz_err)). Comparing these values with the errors arising from high-frequency signal fluctuations (column σ(npt)), we finally obtain that the noise introduced into the nonparametric trend due to incomplete filtering of seasonal effects is almost 100 times lower than the "natural" trend estimation error associated with high-frequency fluctuations of (4) R(t).

Thus, with the chosen parameters of seasonal and trend filters, the annual variations practically do not penetrate into the trend component of the signal. Therefore, it can be concluded that the additional step of the iterative process with a more detailed description of seasonal effects will not increase the accuracy of the nonparametric trend estimation, since the influence of these noises is already negligible.

At the same time, it should be noted that, if the nonparametric trend is estimated from series with an unfiltered seasonal component, the errors induced by this component become comparable by magnitude with other errors in trend estimation. Since the errors have different magnitudes at different times and peak at the boundaries of the series, as well as at the boundaries of the gaps, it can be concluded that the subtrac-tion of (3) S is a necessary procedure and, if it is not used (when using a noniterative algorithm), the trend estimation error can increase significantly.

The cumulative trend can be obtained by adding the nonparametric trend calculated above to the piecewise-linear trend:

(3)

The series (4) Q(t) are shown in Fig. 11 . It should be noted that the constructed model trends of (4) Q(t) do not include the seasonal effect, but they contain variations with periods of more than 1 year. These curves can also be used to estimate the average growth rate of the CO 2 concentration for any period. These estimates are correct for any time period above 2.5-3 years (cutoff frequency of the smoothing filter), even if the length of the interval is not a multiple of an integer number of years, and the beginning and end of the interval fall on different seasons of the year. The estimation of the CO 2 concentration growth rate is also correct if there are gaps in observations, since all the models used allow unbiased estimates to be obtained in this case as well. For this reason, the estimation of even the linear component of the CO 2 concentration growth rate most preferably should be based on the model trends Q(t) constructed by us, rather than on initial data (the series D(t)) and much less on the average annual values of D(t).

In many problems, the CO 2 concentration growth rate V D , rather than the concentration values, are of interest. Estimates for the growth rate can be obtained in different ways, including through direct processing of the original series D(t). However, in view of the

. t Q PLT NPT Table 5 . Infiltration of seasonal variation into the series of the CO 2 concentration trend, constructed by the nuclear sliding smoothing method with a window of 4 years in width and a Gaussian weighting function of the window (edge effects are excluded) σ(res) is the standard deviation of the residual component of the series obtained from the original signal D by the successive subtraction of the piecewise linear and nonparametric trend, as well as the seasonal function (3) S; σ(npt) is the error (standard deviation) of nonparametric trend estimation; σ( (3) S) is the standard deviation of seasonal variation (3) S; σ(trend) is the standard deviation of the trend built by the indicated method from the series (3) S; P is the noise infiltration into the trend series, calculated as the ratio σ(trend)/σ( (3) S); K S is the coefficient of suppression of the amplitude of variations at periods close to a year, when the component (3) S is filtered from the signal; σ,rel is the ratio of the standard deviation of the noise infiltrated into the trend series to the standard deviation (3) S; and σ(Sz_err) is the standard deviation of the noise that may be present in the trend due to incomplete removal of seasonal effects before trend estimation. abovementioned reasons, much more reliable and stable estimates can be obtained from the analysis of Q(t) trends: the CO 2 concentration growth rate V D can be estimated through the derivative of interpolated series (4) Q(t).

To suppress random short-term fluctuations, we additionally smoothed the increments of the trend (4) Q(t) with a sliding 10-year Gaussian kernel. The resulting smoothed series of V Q (CO 2 concentration growth rate) are shown in Fig. 12 . Each point of the graph in Fig. 12 shows the weighted average growth rate at a given station within a symmetric 10-year window with a Gaussian kernel weighting function. Since this weighting gives more weight to the values in the 5-7-year vicinity of the point, changes in the CO 2 concentration growth rate can be seen on scales of at least 5 years. Let us consider some features of these curves.

First of all, attention should be drawn to the fact that the CO 2 accumulation rates almost perfectly coincide for all stations. This confirms the reliability of the measurement results. There is a significant discrepancy between the curves only for 1973-1976. Most likely, this is associated with the difference in the dates of the beginning of observations at the stations. The incomplete suppression of edge effects leads to discrepancies between trends in 1976-1977, when observation data are available for all four stations. At all other times, the discrepancy between the V Q curves for all stations usually does not exceed ±0.1 ppm/year.

The amplitude of differences between stations was numerically estimated from calculated difference series. Each series shows data deviation for the current station from the average over the four stations. The statistics of difference series are given in Table 6 . It can be seen that the largest discrepancies between the estimate for a particular station and the average over the four stations in 1976-2019 do not exceed 0.15 ppm/year. This roughly corresponds to the error expected in view of the (4) Q(t) estimation error, which is 0.005-0.030 ppm/year for different stations. The standard deviation of the difference series does not exceed 0.05 ppm/year. Here, the greatest differences are demonstrated by the BRW station (the most northerly among the stations considered), which is consistent with the results shown in Table 5 (column σ(npt)).

The closest to the average over the four stations is the value given by the MLO observatory; i.e., these are the most representative data.

For a comparison with global variables, it is most interesting to estimate the value of V Q averaged over the four stations (see Fig. 12b ). This curve was obtained with sufficient accuracy: its estimation error (standard deviation) after 1979 does not exceed 0.10 ppm/year, and the typical errors are close to 0.05 ppm/year, which is 2-3% of the estimated value.

It can be seen from Fig. 12b that the rate of CO 2 accumulation in the atmosphere in 1980-1990 was 1.5-1.6 ppm/year. In the early 1990s, there was an indubitable decrease in this rate to 1.2-1.3 ppm/year.

Since the second half of the 1990s, there has been a rapid and steady increase in the CO 2 concentration growth rate. From 2000 to 2015, the smoothed rate of CO 2 accumulation increased from 1.8 to 2.5 ppm/year. However, around the middle of 2016, the rate of CO 2 Years CO 2 , ppm accumulation in the atmosphere stabilized and even began to slightly decline.

Since the constructed V Q curve is rather strongly smoothed, the inflection of this curve can be fixed with an accuracy of no larger than ±3 years. However, the nature of the V Q curve over the 4.5-year interval from mid-2016 to late 2019 leaves no doubt that the trend will change.

It should be emphasized that this study did not use data for 2020, when the СО 2 emissions began to decline due to the global lockdown associated with the coronavirus pandemic. Thus, the rate of СО 2 accumulation in the atmosphere decreased for reasons unrelated to the pandemic.

In general, it can be concluded that the stage of rapid (and ever accelerating) accumulation of СО 2 in the atmosphere in the second half of the 2010s ended and was replaced by a stage with a fixed or even decreasing accumulation rate.

The calculations allowed us obtain trends in the СО 2 concentration change at four stations. The series shown in Fig. 12 can be used to estimate the changes in atmospheric CO 2 accumulation trends.

6. DISCUSSION OF RESULTS Thus, the idea of an algorithm of iterative trend estimation is as follows. At the first step iterations, the simplest and small-parametric model is chosen to estimate the seasonal effects and trend over the entire series using the simplest nonadaptive models. Next, the main part of the seasonal effects is removed from the signal and the refined trend is estimated; then, the refined seasonality model is built, and so on. These iterations can be continued until the estimation errors of Q(t) and S(t) are reduced to acceptable levels. Since the errors decrease with each step, this algorithm converges quickly.

It is important that only a limited number of parameters are evaluated. This reduces the dimension of the problem and increases the stability and reliability of the estimates. More importantly, this algorithm is completely controllable. At each iteration step, the expert receives the estimation results most vividly and can make a decision to correct the model or stop the iterative process, as well check the validity of this decision on the basis of appropriate formal criteria.

However, it should be remembered that the error estimates obtained from the criteria for the internal convergence of the method are always overly optimistic (Deshcherevskii et al., 2016a (Deshcherevskii et al., , 2017a (Deshcherevskii et al., , 2017b . Specifically, these errors are in fact conditional: they are based on the assumption that the chosen method (the model of trend and/or seasonality) adequately describes the process under consideration. In reality, this cannot be known for sure. That is why a careful researcher always analyzes different mathematical models to describe the signal properties using statistical approaches, when such a model cannot be deduced from physical considerations. In this case, a much more objective criterion for the description quality is the difference between models, rather than the internal error estimates within each model considered. Actually, these estimates involve uncertainty associated with the absence of knowledge about the true model of the process.

It should be noted that the nonparametric trend constructed by the sliding nuclear smoothing method (Deshcherevskii et al., 2016b) has exactly the same length as the original signal (Figs. 10, 11 ). This would be impossible with standard approaches, when the trend is calculated only for those window positions that allow the window to be fully inside the series. Clearly, the trend length in this case will be shorter than the length of the original signal by the sliding window size. The disadvantages of this approach are very obvious: each filtering at the signal boundaries "cuts off" increasingly more observational fragments; as a result, the amount of filtered data available for analysis decreases very significantly after a few computational procedures.

The range of window positions can be "expanded" using the two technological techniques proposed in (Deshcherevskii et al., 2016a (Deshcherevskii et al., , 2016b (Deshcherevskii et al., , 2017b . First, the trend is built using the centered smoothing method (when the smoothed value corresponds to the window center), which makes it possible to estimate its values at the ends of the series when the window goes beyond the signal limits by half. Second, the nonparametric trend is built from a signal with a piecewise linear trend that has already been removed. This makes it possible to estimate the trend values at the edges of the series without a significant bias.

This study proposes an iterative algorithm for the decomposition of monitoring data series into trend and residual (including the seasonal effect) components. This algorithm is based on the approaches proposed in (Deshcherevskii et al., 2016b) and makes it possible to obtain unbiased estimates of the trend and seasonal components for rapidly growing series containing observational gaps and omissions. The main idea of the algorithm is that both the trend and seasonal components should be estimated using a signal Table 6 . Statistical characteristics of the difference between the CO 2 concentration growth rate averaged over the four stations and the concentration growth rate at the corresponding station, ppm/year. All curves are smoothed using a sliding window of 10 years and a Gaussian kernel weighting function. that is maximally cleared of any other variations considered noise. For example, in estimating the trend component, seasonal variation is a noise, and vice versa.

At the first step, the linear trend (1) Q L (t) is filtered out from the original signal (0) D(t). The left index (1) indicates that the estimate refers to the first iteration step. This trend model is the roughest, and its parameters are estimated for the entire series as a whole. Therefore, it can be expected that the bias in the estimate of A L (see Fig. 1a ) due to the influence of unfiltered seasonality effects will not be overly large.

Next, the resulting difference signal (0) R(t) = (0) D(t) -(1) Q L (t) is used to obtain the simplest estimate of seasonal effects; for example, an estimate of the smoothed average seasonal function (1) S(t). Despite the roughness, this estimate allows one to take into account largely the variance of the seasonal component of the series and reduce the residual effects of seasonality by an order of magnitude.

Then, the function (1) S is subtracted from the signal. We emphasize that the seasonal effect is removed precisely from the original series that is not cleared of the trend:

(4) This ensures that the results of subsequent calculations are independent of possible estimation errors for (1) Q L (t). It is this technique, when the original series At the second step, the difference signal (1) R(t) with the filtered seasonal variation (1) S(t) is analyzed. This allows us to obtain a more accurate trend model (2) Q(t). In practice, several variants of the function (2) Q(t) are often considered to the best one among them. In this study, we constructed (2) Q(t) as a superposition of piecewise-linear (denoted as PLT) and nonparametric (NPT) trends.

At the third step, the refined trend model (2) Q(t) is filtered out from the signal:

For the series of atmospheric CO 2 concentrations, the trend model (2) Q(t) is already sufficiently accurate to turn to the construction of the final model of seasonal effects (3) S(t). We emphasize that the expert at this stage works with the difference function (3) R(t) cleared of the trend, rather than with the initial signal. Clearly, it is always useful to clear the signal from noise regardless of the estimation method. However, this is doubly important in expert analysis, when one needs to correctly choose the class of models.

At the fourth step, the signal is almost completely cleared of seasonal effects with the help of model (3) S(t). We also removed the piecewise-linear trend (2) Q(t) from the series:

(6) Then, the difference component (4) R(t) was use to estimate the nonparametric trend (4) Q A (t) and obtain the final trend estimate:

The analysis showed that at this stage, the error in estimating the trend (4) Q(t) reaches a possible minimum due to the contribution of the fluctuation (highfrequency) component of the signal. This means that four iterations are enough to completely eliminate the influence of seasonal effects and data gaps on trend parameters. For series with a less pronounced trend, three or even two iterations are usually sufficient.

The iterative approach allows the researcher to be naturally involved in the optimization of models of both trend and seasonal components. The approximation technique provides maximum flexibility and is fully controllable and controllable at all stages of the process. In addition, this technique allows the problems that arise in the presence of observational gaps and defective measurements to be solve in a natural way, without filling them with artificially modeled values.

The algorithm proposed here was tested using the example of the series of atmospheric CO 2 concentration data obtained at four stations of different latitudinal zones. The trends of changes in the СО 2 concentration were obtained and the rates of its accumulation in the atmosphere were calculated. Unlike standard approaches, the method proposed by us allows the estimates to be made for any time intervals, including those that are not multiples of an integer number of years.

An analysis of the rate of CO 2 accumulation in the atmosphere is complicated by the fact that this indicator has a sufficiently high interannual variability. The optimal smoothing level for these data is a debatable issue. Most often, 5-or 10-year averages calculated over calendar intervals are considered. When the CO 2 concentration trend is calculated using the proposed technique, the analysis of changes becomes free of a given discrete time scale. Instead, one can trace the changes occurring with any temporal resolution limited only by the chosen range of frequencies (periods).

The graph of the rate of CO 2 accumulation in the atmosphere (see Fig. 12b ) was constructed by smoothing the station graphs with a sliding 10-year window with a Gaussian kernel, which limits the frequency range under consideration to fluctuations with a characteristic time of at least 5 years. However, the nature of this graph leaves no doubt that the rate of CO 2 accumulation in the atmosphere stabilized around 2016 and even tends to decrease. Before that, a similar stabilization period was observed in 2004-2008. It should be emphasized that this result was obtained from data limited to December 2019, i.e., before the time when anthropogenic CO 2 emissions began to decline due to the pandemic-related effects. The extent of persistence of the revealed tendency and the effect of economic recovery on it will become clear in the next 2-3 years as new data are available.

Radiation variations in the northern part of the Azov Sea

A technique for the preprocessing of nonstationary time series

Models of the spectral density of the flickernoise power

Evaluation of pulse waves as a physical process

Seasonal Adjustment Methods and Real Time Trend-Cycle Estimation

Additivnaya i mul'tiplikativnaya modeli sezonnykh variatsii geofizicheskikh polei (Additive and Multiplicative Models of Seasonal Variations in Geophysical Fields)

Nekotorye voprosy metodiki otsenki srednesezonnykh funktsii dlya geofizicheskikh dannykh (Some Problems of the Method for Estimating the Average Seasonal Functions for Geophysical Data)

Two models of seasonal variations in geophysical fields

A two component model of geophysical processes: Seasonal variations and flicker noise

A flicker-noise problem in the study of cause-and-effect relationships between natural processes

Testirovanie metodiki otsenki parametrov flikker-shuma (Testing the Technique for Flicker-Noise Parameter Estimation)

Algorithm of seasonal variation filtering for geophysical time series

Spectral-temporal features of seasonal variations in apparent resistivity

Flicker noise structure in the time realizations of geophysical fields

WinABD-A universal software package for analyzing monitoring data

Problems of analysis time series with gaps and methods for their solution by the WinABD software

Technologies for analyzing geophysical time Series: Part 1. Software requirements, Seism. Instrum

Technology for analyzing geophysical time series: Part 2. WinABD-A software package for maintaining and analyzing geophysical monitoring data

Analysis of rhythms in experimental signals

Assessment of the seasonal variation in the apparent resistivity of rocks in a well on a Chirkey HPP dam with level shifts and interruptions in observation series

Comprehensive methodology for describing and filtering exogenic effects in monitoring data, taking into account the type of observations and experimental data defects

Technologies of preliminary data processing for multidisciplinary geophysical monitoring and a case study of their application in the Kamchatka geoacoustic observation system

Review of methods for statistical analysis of time series and problems arising in the analysis of nonstationary time series

Identification of the nonstationary cyclic component from time series

Comparison of methods for the seasonal refinement of time series

Extrapolation of nonstationary time series with cyclic components

The spectrum and phases of conjuncture cycles of economic indicators

Digital Filters

Flicker-noise of electronic equipment: Sources, ways of reduction and application

Analysis of methods for filling data gaps in time series of financial market indices, Vestn. Voronezh

Climate fluctuations as flicker-noise

Understanding Digital Signal Processing

Missing data and the accuracy of magnetic-observatory hour means

Variatsii geofizicheskikh polei kak proyavlenie determinirovannogo khaosa vo fraktal'noi srede (Geophysical Field Variations as Manifestation of Deterministic

Multidimensional analysis of time series of geophysical monitoring systems

Analysis of the interaction of geophysical processes in monitoring problems

An aggregated signal of low-frequency geophysical monitoring systems

Geophysical monitoring: Noises, signals, and precursors

Analiz dannykh sistem geofizicheskogo i ekologicheskogo monitoringa (Analysis of Geophysical and Environmental Monitoring Data)

Regularities, characteristics, and origins of the variability of annual and seasonal river runoff in the Ural River basin

Flicker-noise processes in structurallydisordered silicon systems

Geomagnetic Observations and Models

Northern hemisphere temperatures during the past millennium: Inferences, uncertainties, and limitations

Separating a weak periodic component from a nonstationary time series

Flicker-noise and its simulation

Adaptive models for the prediction of nonstationary power consumption time series

Non-stationary time series prediction using one-dimensional convolutional neural network models

Nestatsionarnye vremennye ryady: Metody prognozirovaniya s primerami analiza finansovykh i syr'evykh rynkov (Nonstationary Time Series: Prediction Methods with Cases of Analyzed Financial Commodity Markets

Adaptive analysis of nonstationary time series in studies of seismic oscillations in the range of periods from 0.5 to 5 hours

Local polynomial approximation with a variable-size moving window for the recovery of remote sensing time series data

Variations of the seasonal behavior of geostrophic circulation in the Black Sea

Instantaneous global plate motion model from 12 years of continuous GPS observations

Flicker noises in astronomy and elsewhere

Theory and Application of Digital Signal Processing

Procedures for numerical analysis of circadian rhythms

1/f noise in oscillatory modes of combustion

Ekologicheskoe prognozirovanie: Funktsional'nye prediktory vremennykh ryadov (Ecological Prediction: Functional

Studies in astronomical time series analysis. VI. Bayesian block representations

Base methods for the recovery of gaps in data arrays

Minutes from an Infinite Paradise

Seismoelectromagnetism and spatiotemporal structures

Deterministic Chaos

1/f noise as a nonstationary process: Experimental evidence and some analytical conditions

Atmospheric carbon dioxide dry air mole fractions from continuous measurements at Mauna Loa

Variability of the seasonal course of the Laptev Sea ice cover in the 1940s, in Sovremennye problemy gidrometeorologii i monitoringa okruzhayushchei sredy na prostranstve SNG: Tez. Mezhdunar. nauch.-prakt. konf., posvyashchennoi 90-letiyu Rossiiskogo gosudarstvennogo gidrometeorologicheskogo universiteta (Current Problems of Hydrometeorology and Environmental Monitoring in the Union of Independent States

Continuous uniform observations and scientific progress

Algorithm for the recovery of data gaps in time series in the power consumption prediction system for agglomerating production

Inductive method for the recovery of geomagnetic data time series

Search for electric earthquake precursors by the method of flicker noise spectroscopy

Noise analysis of continuous GPS coordinate time series for CMONOC

WMO Guidelines on the Calculation of Climate Normals

Dynamic and fractal characteristics of time series of seismic energy release

Dynamic characteristics of GPS time series and their relation to seismotectonic regional features

ACKNOWLEDGMENTS Some capabilities of the algorithm described here were illustrated using the example of long-term observation data for atmospheric carbon dioxide concentrations prepared and put on the Internet by K.W. Thoning, A.M. Crotwell, and J.W. Mund. We thank them and all their colleagues who participated in these activities. 

The authors declare that there is no conflicts of interest.