Statistical regularization for trend detection: An integrated approach for detecting long-term trends from sparse tropospheric ozone proﬁles

. Detecting a tropospheric ozone trend from sparsely sampled ozonesonde proﬁles (typically once per week) is chal-lenging due to the noise in the time series resulting from ozone’s high temporal variability. To enhance trend detection we have developed a sophisticated statistical approach that utilizes a geoadditive model to assess ozone variability across a time series of vertical proﬁles. Treating the proﬁle time series as a set of individual time series on discrete pressure surfaces, a class of smoothing spline ANOVA (analysis of variance) models is used for the purpose of jointly modeling multiple correlated 5 time series (on separate pressure surfaces) by their associated seasonal and interannual variabilities. This integrated ﬁt method ﬁlters out the unstructured noise through a statistical regularization (i.e. a roughness penalty), by taking advantage of the additional correlated data points available on the pressure surfaces above and below the surface of interest. We have applied this technique to the trend analysis of the vertically correlated time series of tropospheric ozone observations from 1) IAGOS (In-service Aircraft for a Global Observing System) commercial aircraft proﬁles above Europe and China, and 2) NOAA GMD’s 10 (Global Monitoring Division) ozonesonde records at Hilo, Hawaii and Trinidad Head, California. We illustrate the ability of this technique to detect a consistent trend estimate, and its effectiveness for reducing the associated uncertainty in the noisy proﬁle data due to low sampling frequency. We also conducted a sensitivity analysis of frequent IAGOS proﬁles above Europe (approximately 120 proﬁles per month) to determine how many proﬁles in a month are required for reliable long-term trend detection. When ignoring the vertical correlation we found that a typical sampling strategy of 4 proﬁles-per-month results in 15 7% of sampled trends falling outside the 2-sigma uncertainty interval derived from the full data set, with associated 10% of mean absolute percentage error. We determined that an optimal sampling frequency is 14 proﬁles per month when using the integrated ﬁt method for calculating trends; when the integrated ﬁt method is not applied, the sampling frequency had to be increased to 18 proﬁles per month to achieve the same result. While our method improves trend detection from sparse data sets, the key to substantially reducing the uncertainty is to increase the sampling frequency. separated ﬁt results in trends that can vary considerably from one layer to the next.


Introduction
The vertical profile is a type of time series data that reports the composition (e.g. ozone, water vapor, carbon monoxide) or thermodynamic properties (e.g. temperature, relative humidity, wind speed) of the atmosphere from the surface to an altitude that can range from a few tens of meters (e.g tethered weather balloons) to more than 30 km (e.g. ozonesondes). Conventionally, the analysis of trends based on vertical profile data is conducted on particular altitude bins or pressure levels which are treated 5 as independent time series (e.g., Miller et al. (2006); Harris et al. (2015); Lossow et al. (2019); Petropavlovskikh et al. (2019)), but due to the low sampling frequency of vertical profiles the vertically distributed trend estimates may be noisy or inconsistent from one layer to the next. Trace gas or thermodynamic properties on adjacent pressure or altitude levels are often correlated and the treatment of the levels as independent time series ignores the potential vertical interconnection. In order to reduce the estimation uncertainty of a trend on a particular level, methods have been developed that produce an area average from multiple 10 time series in close spatial proximity (Wigley et al., 1984;Chang et al., 2017), or that employ a dimension reduction technique (e.g. principal component analyses) to project multiple layers of data into a reduced set (Ghil et al., 2002;Meiring, 2007;Park et al., 2013).
In terms of producing an area average based on multiple stations dispersed across a given region, a typical procedure is to create a surface gridded product, which usually aggregates all available monitoring stations within the grid cell without con- 15 sidering spatial sampling and irregularities (Simmons et al., 2010); another example is trend analysis of ozonesondes profiles, based on units of partial pressure with observations integrated into a limited number of layers for the simplification of analysis and calculation of column sum (Tiao et al., 1986;Miller et al., 2006). The dimension reduction technique, e.g. principal component analysis, also known as empirical orthogonal function, is commonly used to reduce correlated multivariate data into uncorrelated vectors that maximize the explained variance with as few vectors as possible. This approach is used to calculate 20 certain climate indices, such as Atlantic and Antarctic Oscillations (https://www.esrl.noaa.gov/psd/data/climateindices/list/). However, this technique is a purely numerical algorithm that is not based on physical principles, thus a meaningful interpretation of the outcome can be subjective or prohibitive.
The aim of such aggregations, either through averaging or dimension reduction, is to enhance a certain common underlying signal (or signals), or to produce regional means (Wigley et al., 1984;Ghil et al., 2002). While the issue of spatial irregularity 25 can now be tackled by a wide range of sophisticated spatial interpolation techniques (Stein, 1999), the averaging approach has not been thoroughly investigated in an objective manner with respect to enhancing the signal to noise ratio. In terms of calculating long-term trends from sparse time series comprised of vertical ozone profiles, this study seeks a statistical framework to accommodate a broad set of correlations, rather than trying to maximize the signal-to-noise ratio (e.g. trend value and its 2-sigma uncertainty) by deliberately selecting a subset of correlated time series. Specifically, this paper aims 30 to extend the traditional multivariate linear regression model by applying a smoothing spline ANOVA (analysis of variance) framework (Wang, 1998;Gu, 2004Gu, , 2013, to account for the vertical profile variability without a dimension reduction, therefore a physical meaning in each pressure layer could be maintained. A smoothing spline is a non-linear curve or surface fitting technique commonly used to model longitudinal data (for example a large cohort study in biostatistics) (Wahba et al., 1995;https://doi.org/10.5194/acp-2019-959 Preprint. Discussion started: 19 March 2020 c Author(s) 2020. CC BY 4.0 License. Verbyla et al., 1999;Shaddick and Zidek, 2015), and it has been widely applied to a variety of spatial pattern and/or temporal trend studies (Dominici et al., 2002;Schlink et al., 2006;Augustin et al., 2009;Park et al., 2013;Wood et al., 2015;Chang et al., 2017). We thus propose a geostatistical approach to fit the vertical profiles jointly as a smooth field for understanding nonlinear interactions and to determine the long-term vertical variability.
Our method uses a two-dimensional geoadditive model for carrying out the multiple correlated time series analysis on 5 every level of the vertical profiles. This approach borrows the strength of spatial correlation models, but instead of working on dimensions of latitude and longitude, we replace the geographical coordinates with a temporal index and altitude bins.
This approach also enables us to identify the common structures among correlated time series (e.g. ozone observations on adjacent layers of the atmosphere) through their seasonality and interannual variability, and to filter out unstructured noise from the underlying signal by a regularization technique. We demonstrate that this method reduces the noise in an irregular 10 and/or sparse time series by taking advantage of the additional information in adjacent layers of the vertical profiles, and yields a trend with reduced uncertainty. While our method was developed for tropospheric ozone profiles it can also be applied to the stratosphere, profiles of other trace gases or atmospheric thermodynamic properties, or to profiles of oceanic or geologic properties.
Section 2 begins with a brief introduction of the time series decomposition and smoothing spline ANOVA framework, 15 we then propose a statistical methodology to assess the profile variability. Section 3 applies this methodology to two tropospheric ozone profile datasets measured by 1) IAGOS (In-service Aircraft for a Global Observing System) commercial aircraft above Europe and China (including a probability test of how the trend precision and accuracy change according to sampling frequency); and 2) ozonesondes launched from Hilo, Hawaii and Trinidad Head, California, by NOAA's Global Monitoring Division (GMD). In Section 4, we discuss the potential extension and benefit of this trend technique, and provide a summary 20 of our ozone profile trend analysis.
2 Model and data

Smoothing spline decomposition
Conventional long-term trend detection of a time series is influenced by many factors such as the number of years of data, the magnitude of the uncertainty and autocorrelation of the noise (Weatherhead et al., 1998). We consider a long-term time series 25 where observations are averaged into monthly means over a number of years in discrete altitude bins or pressure levels. The basic statistical model for an individual time series can be written as: where y(t) is the monthly mean ozone value at time t. A regular seasonal cycle can be estimated by Fourier harmonic series, cyclic spline functions or by simply averaging the monthly values from the long-term climatology. The specification of a trend 30 component is, however, dependent on the goals of the analysis, and it must be flexible enough to represent the underlying structure in the data (e.g. Niu (1996); Augustin et al. (2009);Park et al. (2013); Wood et al. (2015); Chang et al. (2017)); the residual noise series is generally assumed to be an autoregressive process in this scenario.
Ozone variations are in general strongly dependent on altitude or pressure level. For a given altitude or pressure level, our method borrows the information from the neighboring altitudes, at approximately the same time, for a more stable trend estimate (e.g. observations at nearby altitudes are more highly correlated with the layer of interest than those further away). In 5 order to explicitly account for common structures among correlated time series, a class of generalized additive models (GAMs, Hastie and Tibshirani, 1990) for the functional decomposition can form the general framework for joint trend estimates. Additive models extend traditional multivariate linear models by replacing linear covariates with smooth functions, each represented by various types of spline functions (depending on the nature of the covariate) (Wahba, 1990;Wood, 2006;Sangalli et al., 2013;Wood et al., 2015;Chang and Guillas, 2019). Nonlinear variability and complex interactions can thus be modeled by using the 10 combination of multiple spline basis representations within the GAM framework.
With an extension of a single time series to a vertical profile series (let y(t, h) be the monthly mean ozone value at altitude or pressure level h and time t), the statistical model can be extended to a two dimensional signal decomposition: where 15 ω = 1, . . . , 12 monthly index; τ = 1, . . . , n y interannual index (period of coverage); t = ω × τ = 1, . . . , 12 × n y total temporal index; h = 1, . . . , n h altitude bin or pressure level, where the temporal index t is decomposed into monthly index ω and interannual index τ , and f seasonal (ω, h) and f interannual (τ, h) 20 are two dimensional smooth functions representing the seasonal and interannual ozone variations across different altitudes.
Specifically, f seasonal (ω, h) represents a smooth matrix including a single mean for every month and every altitude (each cross section is an annual cycle for the corresponding altitude), and f interannual (τ, h) can be seen as a deseasonalized annual mean anomaly series at each altitude, representing the year-to-year variability. The theory of smoothing spline ANOVA is to decompose a multivariate function into a sum of multiple independent functions with main components and all possible two way (or 25 higher) interactions (up to as many variables as required) (Wang, 1998;Gu, 2004Gu, , 2013, in order to understand the contribution from each single component. In terms of our trend study, we only consider the vertical distribution of two basic elements in the time series: seasonal and interannual variability. A higher order of interactions can be included, but it might not be helpful or informative, and could require a substantial amount of computational complexity. The two dimensional penalized regression splines were chosen for representing both smooth components (Wood, 2003), which can be expressed through a basis representation: This model can be fit to the data and the coefficients can be estimated by minimizing the penalized least square criterion 10 (Wood, 2000) y where · is the Euclidean norm. The first term in Equation (3) minimizes the residual sum of squares (quantifying the goodness of fit). The second term is a nonnegative measure of the complexity of F, which penalizes the roughness of the curve or surface (avoiding over-fitting to the noise). This "regularization" or "roughness penalty" term is defined in terms of 15 an integral of second-order partial derivatives of F.
On the geometry of a function, the second derivative corresponds to the curvature or concavity of the field, thus this penalty term is designed to prevent a "standalone" sharp variability from its neighborhood system. For example, except for the top and bottom layer, the variability in each layer should be bounded by its upper and lower layers, and any variability exceeding these bounds is potential noise. Specifically, the penalty term in Equation (3) can be defined as (Wood, 2003): where {λ 1 , λ 2 } are the tuning parameters for seasonal and interannual components, respectively. A major feature of the roughness penalty term in the regression spline approach is the isotropic property, which means that the roughness penalty is invariant in all directions (i.e. across altitude and time). This isotropic property is natural and straightforward when considering smooth 25 spatial structures on the same geographic coordinates. However, since time and space are measured in different units, we cannot interpret or distinguish the relative importance between smoothness in time and smoothness in space. A pragmatic but efficient approach for avoiding such scaling issues is to transform all variables into the unit hypercube (Sobol, 2001). We considered a potential extension for our particular problem by specifying a different smoothing parameter for each derivative term in Equation (4) (Wood, 2006), however this step proved ineffective as the specification of many smoothing parameters resulted in an over-smooth field. The primary concern in this study is to borrow the strength of common variabilities among neighboring time series for a better trend estimate, with one penalty for the monthly scale and another penalty for the annual scale which are reasonable for our application, as long as the relative vertical distance for the different time series is appropriately documented.
The solution for the penalized least square in Equation (3) is provided by the joint (K 1 + K 2 ) coefficient vector c = (c k1 , c k2 ) , as follows: where X is a (12 × n y × n h ) × (K 1 + K 2 ) covariate matrix that stacks up all smooth basis functions {ψ k1 } and {ψ k2 }, and S 1 and S 2 are known penalized matrices consisting of the smoothing spline basis functions from Equation (4). If λ 1 = λ 2 = 0, the solution is equivalent to the ordinary least square estimate. Extremely large values of the tuning parameter depreciate the curve features and could lead to a plane estimate for F. The estimation of tuning parameter λ = {λ 1 , λ 2 } can be solved by 10 minimizing the generalized cross validation (GCV) score, which is defined as (Golub et al., 1979): where is the trace (sum of the elements on the main diagonal) of A(λ), and N = n h × n y × 12 is the total number of data points. Furthermore, the number of basis functions has to be selected, which is not a part of the optimization problem. The number of basis functions for each 15 component is adjusted cautiously so that the impact on the result is negligible from a further increase of basis functions (here K 1 = 120 and K 2 = 300, a full spline representation would be up to K 1 = n h × 12 and K 2 = n h × n y , but it is considered to be computationally wasteful). Further technical details on spline functions in the penalized least square can be found in Wood (2006).
The specification of our model enables us to decompose the vertical profile into two geoadditive fields: seasonal and interan-20 nual components, according to their associated variabilities. These surfaces are defined by approximating the overall monthly means and annual mean anomalies in each layer or altitude through a statistical optimization. Rather than directly specifying a linear component into a multivariate linear regression model, our functional approach aims to partition monthly mean ozone profiles into components attributable to different sources of variation before any attempt at trend detection.
We choose the thin-plate smoothing spline as the spatial basis function because it is computationally efficient, and because 25 it avoids the problem of choosing "knot locations" (Wood, 2003). Note that the focus of this study is jointly modeling multiple correlated time series at different heights above the same region, which differs from a trend study of a series of measurements with uneven temporal sampling at a single altitude (Tiao et al., 1990;Rehfeld et al., 2011;Weatherhead et al., 2017), or irregularly distributed observations from a monitoring network (Stein, 1999;Chang et al., 2015Chang et al., , 2017. Our problem can be seen as decomposing a two dimensional signal according to the associated variability on a regular grid cell (i.e., based on the 30 monthly and annual index at a fixed altitude bin). Specifically, we implemented the penalized regression splines decomposition using the thin-plate splines available from the R package mgcv (Wood, 2006; R Core Team, 2013).

Ozone profile observations
Our new trend calculation methodology is applied to the long-term tropospheric ozone vertical profile datasets measured by 1) IAGOS commercial aircraft above Western Europe and China, at 50 hPa intervals from 950 hPa to 250 hPa for the period 1994-2017; and 2) ozonesondes launched from Hilo, Hawaii (1982) and Trinidad Head, California (1997, with a vertical resolution of 100 meters from the surface to 15 km. For consistency the vertical coordinate of the ozonesondes is 5 converted from meters to pressure levels. The IAGOS program has measured ozone worldwide since 1994, using instruments onboard commercial aircraft of internationally operating airlines (Marenco et al., 1998;Petzold et al., 2015, http The IAGOS dataset in western Europe is a regional aggregation of profiles from several airports. The profiles are measured by several aircraft and on average there are approximately four profiles per day within the study region. With such a high sampling frequency, the data are expected to be far less noisy than a sparsely sampled ozonesonde time series (one to three profiles per week), and highly representative of the true ozone variability above Western Europe. The IAGOS dataset in China is also a 5 regional aggregation, but the time series has a period of low sampling frequency (2007)(2008)(2009)(2010), with many months containing no data. In contrast, the ozonesonde records at Hilo and Trinidad Head are expected to have greater uncertainty, due to the lower sampling rate of one profile per week. Our first task was an evaluation of our method against the standard linear model using the IAGOS data above Europe, to determine if our method yields trend values with reduced uncertainty. The next step was to examine the impact of sampling frequency on the trend estimate, based on the IAGOS profiles above Western Europe. Finally, 10 we applied our method to the sparsely sampled Hilo and Trinidad Head ozonesonde records to demonstrate the improved trend quantification.  Using the interannual component in Figure 2(b) we calculate the trend value for each individual pressure layer. This integrated fit is then compared to the separated linear fit in each layer based on Eq (1) and the original monthly mean time series 15 (with seasonal cycle and autocorrelation accounted for). Note that the trend result from the integrated fit (i.e. via the smoothing spline ANOVA model) takes advantage of the correlated observations from neighboring layers, thereby reducing the noise in the data set, while the separated fit is obtained by applying a simple linear fit to the original (and rather noisy) time series on each independent pressure level.

IAGOS trends above Western Europe
As shown in Figure 4, these two methods agree very well, because, as shown in Section 3.2, the monthly means are generated 20 from very large sample sizes (approximately 120 profiles each month), which minimizes the uncertainty on the simple linear trend in any given layer. A slight discrepancy can be seen above 400 hPa; to investigate the cause we plot the observed time series and model fitted values for 4 pressure levels in Figure 5.   Figure 4 shows that the ozone trends increase with altitude, ranging from 3-8 ppbv decade −1 above 400 hPa, as determined by the integrated fit method. These strong trends are driven by the ozone in the stratospheric air masses that are frequently sampled at these altitudes. To focus on the ozone trends in the upper tropospheric air masses we removed all of the stratospheric air mass samples from the time series, as determined from the potential vorticity values associated with each ozone profile, and available from the IAGOS data portal (https://doi.org/10.25326/20) (Cohen et al., 2018). Removal of the stratospheric air 35 masses reduces the monthly mean ozone values above 400 hPa by approximately 10 ppbv, and reveals a much more subtle year-to-year variability, especially in the upper levels (supplementary Figure S1). Overall the trends in the upper levels are reduced by more than a factor of two when the stratospheric air masses are removed, but they are still greater than the trends in the mid-troposphere (supplementary Figure S2).

5
Due to ozone's high temporal variability, accurate trend detection may not be possible if the ozone profile sampling rate is low. Three previous studies have estimated the optimal sampling frequency for accurate quantification of tropospheric ozone variations from profiles. Logan (1999) concluded that at least 20 profiles per month are required in tropical and extratropical regions to derive reliable monthly mean ozone values. Saunois et al. (2012) analyzed free-tropospheric IAGOS profiles above Frankfurt (1995Frankfurt ( -2008 and determined that sampling frequencies of 4 and 12 profiles per month approximately result in uncer-10 tainties of 9-29% and 5-15% on the seasonal mean, respectively. Leonard et al. (2017) determined that at least 3-6 ozonesonde profiles per month are required for the estimation of climatological ozone values, in terms of matching the 95% confidence level of the fully sampled monthly ozone means derived by the GFDL-AM3 model.
Here we determine the optimal ozone profile sampling frequency for trends calculated using either the separated fit or the integrated fit methods. In order to investigate the impact of sampling frequency on profile trends, we randomly selected a  Figure S4. 20 If we assume that the trend derived from the full data set (dozens of profiles per month) is the "true" value, the sampled trend is randomly located around the true value as a Gaussian distribution for both approaches. For both the integrated and separated fit methods, the increased sampling frequency results in a diminished distribution of trends. The integrated fit is designed to smooth out the noise in the vertical profiles, and as a result the integrated fit produces a narrower range of trend distributions than the simple separated fit. This case study clearly illustrates that a noisy trend structure results from infrequent sampling. 25 The summary statistics for the sensitivity analysis in terms of both precision and accuracy are reported in Table 1 For the trend accuracy metric, we calculate the mean absolute percentage error (MAPE) for each sampling profile: where p h and a h are the sampled and the assumed true trend values at height h, with the averaged metric reported in Table 1.
The separated fit method improves the accuracy but to a relatively smaller degree than its improvement on the precision. While our sophisticated statistical method clearly improves the accuracy and precision of the trends, greater improvements can be 5 achieved by simply increasing the sampling rate by one or two profiles per month. Figure 7 visualizes Table 1 and shows the sampling rates that allow the precision and accuracy to approach that of the complete data set. At least 4 or 5 profiles per month are required to arrive at a 5% outlier rate with respect to the 2-sigma interval for the integrated or separated fit, respectively. Improved precision becomes negligible at sample sizes greater than 15 profiles per month. In terms of accuracy, a monthly sampling frequency of 14 profiles is required to reduce the error to 5% 10 using the integrated fit, while the separated fit requires 18 profiles per month to meet the same goal.  Figures S5 and S6, respectively, to illustrate the statistical fitting quality.

Application of the integrated fit to IAGOS trends above NE China
The resulting trend distribution above China is displayed in Figure 9. The separated fit finds that the 2-sigma interval of trend 25 values is inseparable from zero at 600 and 300 hPa, while the integrated fit detects positive trends with lower uncertainties throughout the depth of the troposphere, with a larger trend in the boundary layer and a smaller trend at 500 hPa. Compared to the separated fit, which is a typical example of overfitting to the noise at a small scale, our statistical adjustment that borrows information from adjacent layers is effective in terms of a better representation between fidelity and complexity.
3.4 Application of the integrated fit to the Hilo and Trinidad Head ozonesonde records 30 Ozonesondes have been launched from Hilo, Hawaii, continuously since 1982, at an average rate of 3 profiles per month for the first 10 years and 4 profiles per month since 1993; the sampling rate at Trinidad Head has been weekly since August 1997. tropospheric enhancement above mid-latitude western North America in 1998 has been reported previously (Langford, 1999;Cooper et al., 2010), and this feature has been attributed to enhanced stratosphere-troposphere exchange following the strong 1997-1998 El Niño event (Langford, 1999;Lin et al., 2015). to the Hilo plot (Cooper et al., 2019). The separated fit results in trends that can vary considerably from one layer to the next.
For example, from 650 hPa to 350 hPa at Hilo the trend uncertainty range includes zero for some levels, but not for others. In contrast, the integrated fit produces narrower uncertainty ranges because it avoids the influence of extreme events and borrows information from adjacent layers, resulting in all uncertainty ranges exceeding zero above 650 hPa. Supplementary Table S1 25 reports the trend value and associated 2-sigma uncertainty above Hilo and Trinidad Head for the tropospheric column, along with the partial column in the lower, middle and upper troposphere, for the full records and since a common reference year of 2000.
For the trend distribution at Trinidad Head, the overall trend profile is shifted toward negative values due to the positive ozone anomaly in 1998/1999 (as described above), but above 900 hPa the 2-sigma uncertainty range still does not exclude 30 zero. Supplementary Figure S7 shows the trend profile when the time series begins in 2000, to avoid the positive anomaly in 1998/1999. Since 2000 trends are very weak, except for the boundary layer where trends are generally negative.
This paper provides a sophisticated statistical approach for carrying out a joint trend analysis of vertical profile data, as applied to tropospheric ozone observations from IAGOS commercial aircraft profiles above Europe and China, and NOAA GMD's ozonesonde records at Hilo, Hawaii and Trinidad Head, California. Our approach is designed to deliver a consistent trend and uncertainty estimate by functionally decomposing the overall vertical profile into seasonal and interannual components. Instead 5 of fitting the trend component as a single linear term, our technique adopts a two-step data driven approach that, 1) expresses the continuous change as a series of estimated annual mean anomalies in each altitude or pressure level, after accounting for deseasonalization and vertical autocorrelation; 2) and then derives the trend value from the annual mean anomalies according to its underlying structure (a linear approximation of the interannual variability is considered to be appropriate). In this study we introduce roughness penalties on both monthly and annual scales, which are the fundamental elements of time series 10 decomposition.This data-driven approach reveals the underlying trend structure without determining the linear, nonlinear or any other type of polynomial in advance. The success of this approach lies in the smoothing spline decomposition of ozone profile variability. Importantly, the regularization aspect diminishes the influence from extreme observations. Our method is easy to implement under low computational costs, and it offers a spatial visualization for investigating the trend and interannual variability typical of ozone profile data. 15 While this study focused on temporal and vertical variations at a specific location or region, an extension to a spatialvarying vertical process (i.e. incorporation of multiple sites with dimensions of longitude and latitude) is theoretically possible.
However, four-dimensional correlations could be difficult to disentangle, while still maintaining computational tractability. The scaling problem beyond three dimensions could be even more challenging, as it involves multiple measurement units in space and time.

20
The main results of the ozone trend analysis are summarized as follows: 1. The abundance of IAGOS ozone profiles above Europe is sufficient to provide a reliable trend result without the need for applying an advanced statistical technique. Nevertheless, a sensitivity analysis demonstrates that at low sampling frequencies the integrated fit outperforms the separated fit in terms of accuracy and precision. 4. The reference year of the trend estimate is particularly crucial at Trinidad Head, due to high anomalies at 1998-1999.
The trends are weak in the mid-and lower troposphere, and tend to be more negative in the upper troposphere. Finally, although this method was developed for tropospheric ozone profiles it can also be applied to the stratosphere, profiles of other trace gases or atmospheric thermodynamic properties, or to profiles of oceanic or geologic properties. Supplementary Figure S8 illustrates the application of the method to the quantification of trends in the lower stratosphere above Hilo, Hawaii.