Vertical Distribution of Arctic Methane in 2009–2018 Using Ground-Based Remote Sensing

: We analyzed the vertical distribution of atmospheric methane (CH 4 ) retrieved from measurements by ground-based Fourier Transform Spectrometer (FTS) instrument in Sodankylä, Northern Finland. The retrieved dataset covers 2009–2018. We used a dimension reduction retrieval method to extract the proﬁle information, since each measurement contains around three pieces of information about the proﬁle shape between 0 and 40 km. We compared the retrieved proﬁles against Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS) satellite measurements and AirCore balloon-borne proﬁle measurements. Additional comparison at the lowest tropospheric layer was done against in-situ measurements from a 50-m-high mast. In general, the ground-based FTS and ACE-FTS proﬁles agreed within 10% below 20 km and within 30% in the stratosphere between 20 and 40 km. Our method was able to accurately capture reduced methane concentrations inside the polar vortex in the Arctic stratosphere. The method produced similar trend characteristics as the reference instruments even when a static prior proﬁle was used. Finally, we analyzed the time series of the CH 4 proﬁle datasets and estimated the trend using the dynamic linear model (DLM).


Introduction
After carbon dioxide (CO 2 ), methane (CH 4 ) is the second most important anthropogenic greenhouse gas in the atmosphere. While methane is not as abundant, it traps about 28 times more heat per unit mass than carbon dioxide [1]. In comparison to carbon dioxide, methane is relatively short lived with a lifetime of around nine years [2]. Despite this, the concentration of atmospheric methane has more than doubled since pre-industrial times, with global average surface concentration being 1863 ppb in May 2019 [3,4].
Atmospheric methane is produced on the surface by both human activities and natural processes. Major sources of methane include agriculture, waste management, production and use of fossil fuels and different processes in the biosphere involving microorganisms and plants, especially in the wetlands [5]. Increased methane concentrations introduce a potential positive climate feedback in the Arctic in the form of emission changes from the wetlands and potential release of methane trapped in the permafrost caused by a warmer climate [6,7]. crucial for vertical distribution research; most aircraft are unable to reach altitudes above roughly 13 km [26] and therefore cannot measure the stratospheric profile that is highly variable for methane, especially at high latitudes inside the polar vortex [14].
Another way to measure atmospheric profiles is using a satellite with limb viewing geometry, which allows the satellite to scan through different layers of the atmosphere as the sun (or other star) is ascending from behind the horizon or descending behind it. This method is used by the FTS instrument onboard the Atmospheric Chemistry Experiment (ACE) satellite to retrieve atmospheric profiles of CH 4 in the upper atmospheric layers. In this work, the ACE-FTS profile retrievals are used to evaluate our ground-based profile retrievals. A large ensemble of co-located ACE-FTS profiles is also used in the construction of the extra-tropospheric part of prior profiles.
In this work, we retrieved a time series of methane profiles from high-resolution short-wave infrared spectra measured by a ground-based Fourier Transform Spectrometer (FTS). Retrieval was done by modeling the temperature and pressure broadened absorption lines in short-wave infrared (SWIR) wavelength range using the Swirlab code that is publicly available [27]. Other remote sensing profile retrieval codes for ground-based FTS spectra were described and demonstrated by Hase et al. [28] and Zhou et al. [29]. The inverse problem related to the profile retrieval is ill-posed, thus, to retrieve physically meaningful profiles and to prevent artificial oscillations in the retrieved methane profile, the retrieval needs to be constrained. The other codes, PROFFIT ( [28]) and SFIT4TCCON ( [29]), use Tikhonov regularization [30] to obtain an appropriate prior for the retrieval. Swirlab on the other hand utilizes prior based dimension reduction [31], where the retrieved profile is constructed using only the largest singular vectors calculated from the prior covariance matrix. The prior itself is first built to reflect the mean and variability of the ACE-FTS measurements. An advantage of this method is that, in addition to regularization, the dimension of the problem is reduced, leading to a more computationally efficient retrieval algorithm [32].
The most significant advantage of the FTS measurements over the other profile measurements is the time resolution. The aircraft and balloon-borne measurements are superior in accuracy and vertical resolution. However, they are time consuming and resource intensive, and therefore not made on a daily or even weekly basis. Solar occultation measurements by ACE-FTS also have a higher vertical resolution than the ground-based remote sensing methods but the profiles cover only altitudes from roughly 10 km upwards. Due to the satellite's orbit and measurement geometry, the revisit time over Sodankylä varies between 1 and 70 days. Temporary data gaps in the ground-based FTS measurements are caused by clouds blocking the direct view of the sun. Excluding the polar night, data gaps longer than two days are very rare.
While the profile information is useful in itself, the profile retrieval method could also improve the total column measurements. The averaging kernels show better sensitivity to the variability in the stratosphere than in the case of a profile scaling method [32]. In addition, the solar zenith angle (SZA) dependence of the averaging kernels is lower, thus resulting in reduced SZA-induced biases when the true profile differs significantly from the prior profile used for the retrieval [32].
To evaluate the performance of Swirlab retrievals, the profiles were here compared to concurrent ACE-FTS profiles above 10 km altitude, AirCore balloon measurements below 30 km and to 50-m mast measurements for the lowest atmospheric level of the profile. In addition to one-to-one comparisons, a dynamic linear model (DLM) trend analysis [33] was carried out to each dataset, and the growth rates of the recent years were compared. The growth rates were also evaluated in the light of recently published Arctic methane growth rates [7,10,18].

Data Description
The study concentrated on data measured at or above Finnish Meteorological Institute's Arctic Space Centre in Sodankylä. The site is located at 67.37 • N, 26.63 • E (see Figure 1) and hosts a high-quality measurement infrastructure covering several research areas such as upper-air chemistry and physics, atmospheric column measurements, cryospheric and hydrological processes, biosphere-atmosphere interaction and satellite calibration and validation.

Ground-Based FTS Dataset
The Fourier Transform Spectrometer at Sodankylä became operational in 2009 [34]. The instrument participates in the Total Carbon Column Observing Network (TCCON). The TCCON facility at Sodankylä is equipped with a Bruker IFS 125HR spectrometer and a large solar tracker A547N (Bruker Optics, Ettlingen, Germany). The instrumentation is designed to record high-resolution solar spectra in the near-infrared spectral region. From the spectra, column-averaged abundances of CO 2 , CH 4 , N 2 O, HF, CO, H 2 O, HDO and other gases are derived.
The FTS instrument at Sodankylä has two room temperature detectors: an indium gallium arsenide (InGaAs, wavenumber range is 4000-11,000 cm −1 ) and a silicon diode (Si, range is 9000-15,000 cm −1 ). In addition, a liquid nitrogen cooled indium antimonide detector (InSb, detector covers 1800-6000 cm −1 ) was installed in July 2013. The additional detector expands the wavelength region covered by the instrument and provides a possibility to retrieve more atmospheric species. For stability and to keep water vapor low, the measurements are performed in vacuum. For the TCCON measurements, the spectral resolution is 0.02 cm −1 , with an optical path difference of 45 cm and the collection time for a single scan is 78 s. The instrument has been operated in an automated mode: after two InGaAs/Si scans, one InSb scan is taken. The instrument characterization has been performed using a hydrogen chloride gas cell.
For the TCCON retrieval, GGG2014 algorithm is used [35]. The algorithm retrieves total column dry air mole fraction of methane using three wavelength windows centered at 5938, 6002 and 6076 cm −1 (1684, 1666 and 1647 nm). The spectroscopy is based on a modified HITRAN2012 line list [36]. The National Centers for Environmental Protection and National Center for Atmospheric Research (NCEP/NCAR) reanalysis data are used to generate prior profiles of pressure, temperature, geopotential height and water vapor. Methane prior profiles are generated by empirical functions that include a secular increase, interhemispheric gradient, seasonal cycle and stratospheric decay based on the age of air. These functions are optimized to fit existing information from balloon-borne platforms, satellites and aircraft campaigns. The prior CH 4 profiles are iteratively scaled to generate forward-modeled spectra that best fit the measured spectra. To calculate the total column dry air mole fractions (XCH 4 ), the total dry air columns are calculated from retrieved O 2 columns, with a fair assumption of constant O 2 dry air mole fraction (20.95%). An empirical correction is applied to correct the solar zenith angle (SZA) dependence, determining any variations that are symmetric about local noon as SZA-dependent artefacts [21]. The TCCON data are publicly available at the TCCON data archive, hosted by Caltech-DATA [37]. The algorithm to retrieve profile information from the same spectra is described in Section 3.1.

AirCore Measurements
The AirCore was constructed as a 100-m-long coated stainless steel tube, with a volume of about 1400 mL. A data logger was built to record the temperature of the AirCore and the ambient pressure. An automatic valve was used to close the inlet valve of the AirCore system shortly after the landing of the payload. The AirCore was lifted by a meteorological balloon to the altitude of typically 25-30 km, and was filled during descent from the stratosphere down to the surface. The gas analysis was performed typically within 2-3 h after the landing of the AirCore using a cavity ring-down spectrometer (CRDS). Calibration of the CRDS was performed before and after analyzing the sample to ensure the traceability of the measurements to the World Meteorological Organization scale of CH 4 (X2004A).

ACE Data
Atmospheric Chemistry Experiment (ACE), also known as SCISAT, is a Canadian-led satellite mission for remote sensing of the Earth's atmosphere. It was launched in August 2003. The primary instrument onboard is the infrared Fourier transform spectrometer (ACE-FTS). It measures infrared spectrum in the wavenumber range of 750-4400 cm −1 with a resolution of 0.02 cm −1 . The spectra are measured during sunsets and sunrises in the limb viewing geometry with different slant paths and tangent heights from upper troposphere up to 150 km [38]. From the measured spectra, several tens of microwindows are analyzed to retrieve CH 4 profiles with the vertical resolution of about 4 km. In the comparison studies, the accuracy has been found to be within 10% in the upper troposphere-lower stratosphere region [39].
To gather a comprehensive dataset, the co-location criteria for ACE-FTS measurements was set quite loose. Measurements within 30 degrees in longitude and 3 degrees in latitude from Sodankylä were included in the analysis (see map in Figure 1). The comparison between ACE-FTS and Swirlab was done on the spatially closest ACE-FTS profile of the day and the daily average of SZA-corrected Swirlab profiles. We used the Level 2 ACE-FTS data version 3.5 and 3.6 interpolated to a 1 km grid with center altitudes from 0.5 to 149.5 km. The analysis was restricted to altitudes from 10 to 40 km as the ACE-FTS data below 10 km are sparse and the sensitivity of the Swirlab retrieval is very low above 40 km.

Continuous Mast Measurements
Sodankylä 50-m-high mast has gas line inlets at four different heights: 2, 10, 23 and 50 m a.g.l. Sample air from highest level inlet is measured for 50 min/h, and air from two other (2 and 23 m) levels for 5 min/h each. The 10 m level is not measured regularly. Gas lines are 8-mm OD stainless steel, and they are flushed with flow rate of ca. 2.5 L/min, which corresponds to less than 1-min residence time in the lines.
Sample air is measured with a CRDS instrument, which measures CO 2 , CH 4 , CO and water vapor concentrations at the same time. This instrument is automatically calibrated with WMO/GAW/CCL (World Meteorological Organization/Global Atmosphere Watch/Central Calibration Laboratory, i.e., NOAA ESRL GMD) traceable calibration gases ca. 2 times a month. In addition, two cylinders, "working standard" and "target", are measured every 11 and 34 h, respectively. The methane scale in use is WMO X2004A.
Humid air samples are measured directly, i.e., no dryer is used. This is compensated by determining and using individual water vapor correction factors for the instrument [40,41]. The raw data are first calculated as 1-min averages, and further to uncalibrated dry concentrations based on water vapor correction factors. Calibration is applied to these values, and hourly average concentrations are calculated from 1-min calibrated dry concentrations. For the time series analysis, a daily average value is calculated from the hourly values.

Profile Retrieval Method
Traditionally, ground-based direct solar radiation measurements at SWIR wavelengths have been used to retrieve methane total column. This is done by scaling the prior profile to produce a modeled spectrum that has the closest match to the measured spectrum (e.g., the method used in the TCCON retrievals [21]). However, the shapes of absorption lines include information about the vertical structure of the methane profile. The lines are broadened by temperature and pressure, resulting in wider absorption features at low altitudes and narrower at high altitudes. The algorithm described in this section finds the best fit between the modeled and measured spectra by allowing the methane profile shape to vary. This allows better fit between the measured and modeled spectra as also the width, in addition to the depth, of the absorption features in the modeled spectra is modified. The retrieved profiles, the vertically resolved concentrations and their changes can be studied. There are also advantages in using this method for total column retrievals: the averaging kernels show better sensitivity to the variability in the stratosphere than in the case of profile scaling method [32]. In addition, the SZA dependence of the averaging kernels is lower, thus resulting smaller SZA induced biases when there is a large difference between the shapes of the true profile and the prior profile used for the retrieval [32].
To retrieve vertical profiles of methane, we used the dimension reduction method described in detail in [32]. In this section, we give a brief description of the method. The goal has been to develop a method that would be only weakly dependent on the prior assumptions but would still provide information about the vertical structure and realistic profiles. The entire algorithm, hereafter referred to as Swirlab, is programmed to be run in Matlab and is publicly available [27].
The measured FTS data, described in Section 2.1, allows us to retrieve around 2-3 independent pieces of information about the vertical structure of CH 4 , as demonstrated in [32]. This means that we can actually extract more information from the measurement than the operational TCCON retrieval, which only has one degree of freedom. Thus, we construct the inverse problem so that there are 2-3 estimated parameters only. This kind of parameterization simplifies and stabilizes the inverse problem, and eases the computation in comparison to a layer-by-layer retrieval. The parameters need to be chosen so that they manifest the information content of the measurement, and we need a way to project the low-parameter profile back to the full space.
The profile of the atmospheric constituent of interest, in our cases CH 4 , is defined by first discretizing the atmosphere from 0 to 70 km into n layers and identifying the relative contribution of methane to the total atmospheric density (in parts per billion, ppb) per layer as x ∈ R n . The inverse problem of inferring methane concentrations from measured data y ∈ R m , also known as a retrieval, is derived from the relation where f is a non-linear atmospheric radiative transfer model and ε ∈ R m is a noise term describing the measurement and modeling error. The measured spectra used in this work are the same used for the TCCON retrieval but from a small microwindow with a wavenumber range from 6003 to 6005.5 cm −1 including five methane absorption lines. An example spectrum is given in Figure 2. We use the Bayesian analysis for solving the retrieval problem and assume Gaussian distributions for the prior and noise, x pr ∼ N ( x 0 , C), ε ∼ N (0, C y ), where x 0 ∈ R n is a priori mean profile, C ∈ R n×n positive definite prior covariance matrix and C y ∈ R m×m is the measurement error covariance matrix. The solution is obtained as posterior density of form We obtain the retrieved methane profile x from the retrieval method R as where R includes our prior assumptions, forward model and Bayesian formulation as well as additional parameters θ such as pressure, solar parameters, etc. These parameters can be retrieved as a part of the state vector, but to keep the notation compact we keep them fixed. The model atmosphere consists of a non-scattering atmosphere with 100 spherical layers from ground level up to 70 km. Absorption lines are calculated for each layer separately as Voigt profiles taking into account the temperature and pressure of the layer. Atmospheric model data from National Centers for Environmental Prediction (NCEP) are used as source for these atmospheric parameters.
As an update to description in Tukiainen et al. [32], the absorption coefficients are calculated using HITRAN2016 line database [42].
The dimension reduction approach that we adapted is based on those of Marzouk and Najm [43] and Solonen et al. [31]. The fundamental idea is to constrain the problem to a subspace that contains most of the variability allowed by the prior. This is achieved by using a smoothing prior covariance that has only a small amount of singular values significantly greater than zero, and thus most of this variability can be represented by a linear combination of the corresponding leading singular vectors. Thus, we start with the singular value decomposition of the prior covariance matrix where the diagonal matrix Λ contains the singular values of the prior covariance matrix and U contains the corresponding singular vectors u 1 , ..., u n as columns. The dimension of the problem is then reduced by using the k (1 ≤ k < n) largest singular values and the corresponding singular vectors and representing the unknown as a linear combination Here, the new unknown z k is a k-dimensional vector whose prior distribution is a k-dimensional Gaussian N (0, I k ) and P k ∈ R n×k , k < n, is a projection matrix from R k to R n . Inserting this new parameterization to the problem leads to a posterior distribution for the subspace parameter z k of the form The prior mean, x 0 , is estimated by combining the stratospheric part from the ACE-FTS satellite data and the tropospheric part from the TCCON prior data. The prior is calculated separately for three domains: summer (from May to August), polar vortex and the rest (denoted as winter). The location inside the polar vortex is defined using ERA-Interim atmospheric model data. The station is assumed to be inside the polar vortex when potential vorticity at 475 K potential temperature level exceeds 3×10 −5 K m 2 s −1 kg. The seasonal prior means are depicted in Figure 3. A trend of 5.57 ppb per year, determined from the TCCON prior profiles 2009-2016, is added to each layer. In addition, retrieval with trendless prior is made and the effect is discussed in Section 4.2. For the prior covariance construction, we follow the procedure from Tukiainen et al. [32] to produce a variability that matches well to the empirical covariance derived from an ensemble of ACE-FTS satellite based measurements. In short, this is done by first computing the diagonal elements of C as where h are the atmospheric layers and the parameters σ 1 , σ 2 , h 1 , h 1 , s 1 and s 2 are chosen appropriately (σ 1 = 300, σ 2 = 30, h 1 = 25 , h 2 = 5, s 1 = 10 and s 2 = 5). The off-diagonal terms are further given by where i and j are two atmospheric layers, dist(i, j) is the distance between these layers in kilometers and l is an appropriately chosen correlation length (l = 12 km). The full covariance matrix and the first four singular vectors are presented in Figure 4.

Solar Zenith Angle Dependence
Due to errors in radiative transfer modeling and spectroscopy, and because the measurements sensitivity to atmospheric composition at different altitudes changes when solar zenith angle (SZA) changes, there is a systematic difference between volume mixing ratios (VMR CH 4 ) measured at noon and at early morning or late evening (see Figure 5). This difference is compensated for by going through data for days when the SZA has had both small and large values. The measurement with the smallest SZA for the day is regarded as the reference value (REF CH 4 ) and the ratio (R) between all the daily measurements and the reference value is calculated.
The R is first calculated for days with the minimum SZA smaller than 44.5 degrees. Before calculating the ratios for the rest of the days, an exponential function is fitted to ratios R. For the days with the minimum SZA larger than 44.5 degrees, the reference value is corrected with this function and the following ratios are defined as where SZA ref is the smallest SZA of the evaluated day. The days are evaluated in order of minimum daily SZA. After each evaluated day, the f corr is updated and the evaluation ends when 10% of the time series (about 13,000 measurements) have been evaluated. Then, the final f corr is regarded as the zenith angle correction function for the time series. The correction function is calculated for each retrieved layer separately. The final corrected volume mixing ratio (cVMR CH 4 ) at layer i is then calculated as As an example, the SZA dependence of the volume mixing ratios at two different layers are depicted in Figure 5. The Figure also depicts the fits for all the layers and the spread of the difference between data points and the fit. The variability of the data points is described with interquartile range as there are some outliers skewing the distribution so that standard deviation would not describe the data well. Assuming that on average the methane concentration is stable and should not vary as a function of the solar zenith angle, it seems that at the layers close to the ground the mixing ratios are underestimated when the SZA increases, and in the stratosphere the values are overestimated when the sun is close to the horizon.

Averaging Kernel Correction to ACE Data
To compare the atmospheric profiles obtained from ACE-FTS measurements with Swirlab retrievals, the results must first be projected into the same space. This is because the ill-defined inverse problem is regulated by prior information and the retrieval is hence a biased estimator of true CH 4 . Moreover, the instruments naturally have different sensitivities at different altitudes, which needs to be accounted for as well. A widely used way to make two separate measurements by different instruments comparable is smoothing using the averaging kernel which measures the sensitivity of the retrieved valueˆ x to the actual true atmospheric state x. Averaging kernels for z parameters can be estimated from the known Swirlab measurements using Jacobians, i.e., the derivatives of forward model f of the form where A z ∈ R k×n , C y ∈ R m×m is the error covariance matrix, I k ∈ R k×k is the prior covariance of z k , and l ∈ R n contains the slant lengths of different atmospheric layers (for details, please refer to [32]). The full-state averaging kernel can then be computed as where A ∈ R n×n . The averaging kernels are SZA dependent, and thus can partly contribute to diurnal and also seasonal biases. A Swirlab averaging kernel matrix at SZA = 45 • is visualized in Figure 6. We see that the changes in the lowest altitudes have a large effect to the higher levels also. As the covariance structure allows the profile to vary the most at altitudes between 20 and 30 km, many of the changes at levels that are more restricted to follow the prior propagate to these altitudes.  Using the Swirlab averaging kernel, we follow [44] and get the smoothed ACE-FTS profile as

Time Series Analysis
In this study, our main interest is the altitude-dependent trend in the Sodankylä FTS time series. In general terms, trend can be defined as a change of distributional processes (e.g., mean) of the process that generates the observations. Typically, the interesting part is the more slowly varying background level after better-known variabilities, such as seasonality, have been accounted for. For time series trend analysis, we follow the approach by Laine et al. [45] and apply dynamical linear regression to the FTS data by using Dynamic Linear Model (DLM) toolbox for Matlab by Marko Laine [46]. In contrast to the more usual linear regression, dynamical regression allows the regression coefficient to evolve in time. This evolution is controlled by model variance parameters that are estimated from the data. DLM in extremely useful in atmospheric time series analysis since a lot of factors, such as external forcing, make the background mean trend nonlinear. On top of this, atmospheric processes are often a combination of short term rapid changes, such as seasonal cycles, and more long term background processes.
We present a short description of the DLM formulation used in our analysis. DLM with Gaussian errors is best described as a general linear state space model with an observation equation and an evolution equation as where y t contains the observations and u t the hidden and unobserved model state at time step t, the observations are gained via the observation operator H t , the state is evolving along the model evolution operator M t . In this formulation, the observation error v t and the model error w t are assumed to be zero mean Gaussian with covariance matrices given by R t and Q t , respectively. The hidden state vector u has components for each process that is fitted in the time series analysis. In our case, we get where µ t is the parameter for the local mean, α t is for local trend, β t,1 and β t,2 represent the seasonal components and η t is the autoregressive component of the state vector. In DLM, the background trend is modeled as a random walk process defined by local mean µ t and local trend α t . For these parameters, we can write the observation operator, model evolution operator and model error as where σ 2 trend is the error variance of the local trend. Similarly, the annual seasonal cycles are modeled using harmonic functions. For seasonal components β t,1 and β t,2 , we thus have M seas = cos(2π/s) sin(2π/s) − sin(2π/s) cos(2π/s) , where s = 365.242 is the number of days in a year (accounting also for leap years) and σ 2 seas is the error variance of the seasonal components. In this work, we also allow autocorrelations in the residuals by using a first-order autoregressive model. The corresponding observation operator, model evolution operator and model error are given by Combining the independent parts for trend, seasonality and autoregression, we get the full observation operator, model evolution operator and model error for the hidden state vector u as With the previous formulation, the DLM algorithm proceeds with Kalman Filter and Simulation Smoother to generate an ensemble of states that fit the observed data. This step is repeated with Markov chain Monte Carlo (MCMC), sampling the parameters σ trend , σ seas , σ AR and ρ. This procedure explores various choices of tuning parameters for the problem and gives the values that fit the observations the best, minimizing arbitrary fixing that is often done manually in other applications. For a detailed description of the Simulation Smoother and the MCMC, we refer the reader to the documentation in [45,46].

The Time Series
The Swirlab profile retrieval algorithm was used to retrieve methane profiles from all the SWIR spectra measured when non-obstructed sunlight was available in Sodankylä between May 2009 and November 2018. For further analysis and illustrative purposes, a time series of daily averages was calculated. Figure 7 shows a surface plot of the time series. In the figure, gaps shorter that seven days are filled and only few data gaps remain during the measurement season. The measurement season covers the time when the daily maximum solar elevation angle exceeds eight degrees, in Sodankylä roughly from mid-February to the end of October. On average there were 148 days with measurements per year between 2010 and 2018. The algorithm produces nicely the seasonal variability with summer minimum in the lower troposphere and really low values in the stratosphere throughout the spring when Sodankylä lies beneath the polar vortex. The seasonal behavior and trends at various altitudes are more closely studied in Section 4.6.

Sensitivity to Prior Profile
To demonstrate the profile information content in the ground-based FTS spectrum, we ran the Swirlab retrieval with different prior profiles. Even when using a single prior profile (Figure 8), the algorithm is able to detect the dip in stratospheric methane concentration inside the polar vortex. The use of a seasonal prior lets the retrievals inside the vortex reach a bit lower level. The introduction of the seasonal prior for summer also makes Swirlab retrieval agree better with the ACE-FTS mixing ratios (see (Figure 8) when the retrieval with a single prior leaves the summer values somewhat low compared to ACE-FTS.   Figure 9 shows the difference in total column mixing ratios calculated from Swirlab profiles with and without a trend in the prior. Swirlab is able to retrieve a positive trend in the methane total column even when the trend is excluded from the prior. However, the trend is roughly 0.3 ppb/year smaller than the one retrieved when a prior trend of 5.57 ppb/year is introduced. The prior trend is the average trend of total methane column calculated from TCCON prior profiles. The total column methane retrieved with prior trend agrees well with the TCCON data (see Section 4.5). Figure 9. Difference between the total columns calculated from the Swirlab profiles when the retrieval is made including or neglecting the trend in the prior profile.

Profile Comparison to AirCore Measurements
The concurrent Swirlab and AirCore profiles are shown in Figure 10. The AirCore profiles were obtained during 2013-2017. During multiyear flights, all seasons have been represented. In 2017, the flights were performed within ESA FRM4GHG Project ( [47,48]). There is, in general, a good agreement between AirCore and Swirlab posterior profile. For some retrievals, another prior profile could have resulted in a better agreement. For example, on 24 and 26 April 2017, the polar vortex prior would follow the AirCore profiles more closely but the potential vorticity in the ECMWF data suggests Sodankylä was not inside the polar vortex so the winter prior profile is used. Nevertheless, the methane deficit in the stratosphere is quite well detected by Swirlab even when using the winter prior profile.
For quantitative analysis, we calculated the root mean square error (RMSE) of Swirlab retrievals against the AirCore profiles. The RMSE is defined as where x is the Swirlab profile, z is the AirCore profile and N is the number of layers where we have co-occurring data points. We performed the same RMSE analysis also for another set of profiles, retrieved using the prior scaling method from [32]. The prior scaling retrieval is similar to the one used in TCCON and thus serves as a good reference point. The results are presented in Table 1. We can conclude that, in the sense of RMSE, Swirlab tends to give results that are on average closer to the AirCore profile. In addition, looking at the measurements where the profile scaling yields a smaller RMSE than Swirlab, the differences between the methods are still relatively small. In contrast, in the majority of cases where Swirlab provided a better retrieval, it was significantly better. We would also like to note that, as can be seen in Figure 10, the retrievals where profile scaling performs better are generally times when the measured AirCore profile is very close to the prior profile.

Profile Comparison with ACE-FTS Instrument
Swirlab retrievals were compared to both original and averaging kernel corrected ACE-FTS profiles ( Figure 11). The comparison was made between the closest ACE-FTS retrieval of the day within the co-location region (see Figure 1) and the concurrent Swirlab retrieval over Sodankylä. For the averaging kernel correction, the missing values in ACE-FTS profiles have been filled with the Swirlab prior profile. The gap-filling method might partly explain the small negative bias below 25 km and small positive bias above 25 km when comparing averaging kernel corrected ACE-FTS to Swirlab.
A large set of ACE-FTS profiles was used to create the prior means for Swirlab and also the prior covariance was built to resemble the empirical covariance of ACE-FTS profiles. Therefore, it was expected that on average ACE-FTS and Swirlab would agree well. However, the information from ACE-FTS for the priors is condensed to three seasonal prior means. The variability around these priors, even if constrained to a realistic and smooth solution by the covariance, is interpreted from the spectra. In this regard, using ACE-FTS as a reference is justified.
The raw, non-smoothed, ACE-FTS profiles and Swirlab profiles agree well below 25 km with Swirlab retrievals staying within 25% from the satellite data. Between 25 and 35 km the profiles agree well on average but the variability becomes larger. The distribution of differences against non-smoothed ACE-FTS profiles is slightly asymmetric with the median difference being close to zero. The median suggests that there is an equal amount of cases where Swirlab concentrations are lower than cases where they are higher than ACE-FTS values. The asymmetry suggests that, when ACE values are lower, the discrepancies between the two retrievals are larger. This could be caused by Swirlab being unable to sense large methane deficits occurring on a limited altitude range that can occur close to the polar vortex edge when air from outside and inside the vortex are mixed [14]. Close to 40 km altitude, the Swirlab retrievals get closer to prior profiles and on average compare well with ACE-FTS but the variability ACE-FTS can detect is much larger than what is possible with Swirlab.
The Swirlab profiles compare well also to the averaging kernel corrected ACE-FTS profiles, the distribution of the differences resembling the comparison to raw ACE-FTS profiles but with slightly smaller variability. However, the asymmetry above 25 km is even larger than in the comparison with raw ACE-FTS. This again could also, at least partly, be caused by the gap-filling with the Swirlab prior information that might not represent the true profile. At levels close to 40 km, the variability becomes smaller as the sensitivity of Swirlab retrieval is reduced and the averaging kernel smoothed ACE-FTS profiles are forced closer to the prior profiles as are the Swirlab retrievals.

Total Column Comparison to the TCCON Algorithm
We used the pressure profile and water vapor profile from NCEP/NCAR data to calculate the total column dry air mole fraction (XCH 4 ) from the Swirlab profiles. The daily averages of TCCON and Swirlab total columns are compared in Figure 12. The agreement is good with a mean bias of −2.8 ppb and standard deviation of 6.1 ppb. The trend of TCCON is followed well by Swirlab data, and based on the lack of seasonality in the residual time series, the SZA correction implemented to the Swirlab data is sufficient. A small step in the difference can be seen since the beginning of 2015 but its origin is unresolved.
The main difference of the two algorithms is the degree of freedom in finding the methane profile that would produce the measured spectra. TCCON scales the prior shape and Swirlab uses four singular vectors to modify also the shape of the profile. This improves the accuracy in cases where the prior profile significantly differs from the true profile. In addition, the residuals in the fit between the measured and modeled spectra become smaller as the absorption line shapes are allowed to vary [32].
There are also other differences. The wavelength range used by TCCON is wider as Swirlab uses one short microwindow (see Figure 2), whereas TCCON has three wavelength windows for its methane retrieval. The algorithms use slightly different absorption line data as Swirlab uses HITRAN2016 [42] and GGG2014 uses a modified HITRAN2012 [36]. Differences can also be partly attributed to the different way of defining the dry air column. As the Swirlab algorithm retrieves neither the O 2 column nor the water vapor column, the dry air column used to derive the dry air mole fraction (DMF) is calculated from the NCEP/NCAR model data. In TCCON retrieval, the air column is calculated using the retrieved O 2 column assuming that the DMF of oxygen is constant. The use of O 2 column for air column estimation improves the precision as it reduces the effects of instrumental or measurement errors that are common to both gases CH 4 and O 2 [49].

Ch 4 Trends
The CH 4 profiles retrieved with Swirlab allow us to quantify the growth rate of atmospheric CH 4 in every level during the 10-year time series. For this analysis, we drew samples from the posterior distribution of the level parameter µ t in our DLM analysis. For each sample, we computed the yearly trend at a timestep t as a change where s is the number of days considered. To assess the trend and the uncertainty, we drew 1000 samples and calculated the average trend and the standard deviation of the trend for every time step. The Swirlab CH 4 time series at four different levels (at the ground level, 11.5 km, 21.5 km and 26.5 km) and respective growth rates are shown in Figure 13. In the troposphere (upper panels in Figure 13; orange color), the results show a positive growth rate throughout the time series with moderate interannual variability, in particular at the ground level. The maximum growth rate of 8.1 ± 1.5 ppb/year occurs around winter 2013-2014 at the ground level, and also simultaneously higher in the troposphere. This is also clearly shown by Figure 14 (left panels), where the growth rate in all tropospheric levels increases until 2014, after which it begins to decrease but remains positive. It is notable that the growth rate has decreased significantly since 2017 at the ground level; however, this cannot be seen with confidence higher in the troposphere. The uncertainty estimates of the growth rate increase with altitude, as well as towards the edges of the time series. In the stratosphere (lower panels in Figure 13; orange color), the seasonal variability of CH 4 is stronger, as well as the interannual variability; thus, weak trends are more challenging to quantify reliably. The variability of the growth rate is mostly smaller-scale than the estimated uncertainty range. The average stratospheric trend is estimated to be weak, and above roughly 20 km altitude, systematically weakly negative; however, the uncertainty estimates are large. Contrary to the ground level and the troposphere, the stratospheric CH 4 growth rate based on Swirlab retrieval levels has been increasing since 2017 and even earlier in the lowest part of the stratosphere. In addition, the Swirlab data show a weak positive trend in the levels close to 40 km but at these high altitudes this is likely due to Swirlab retrieval closely following the prior with a predefined growth rate of 5.57 ppb/year.
We evaluated the Swirlab CH 4 trend signatures described above through comparisons with appropriate reference data at different atmospheric levels. Figure 13 (upper left panel) shows a comparison of Swirlab to the 50-metre mast measurements. In general, the time series from the Swirlab retrieval agrees well with the in-situ measurements, albeit the local variability in the concentration is stronger for the mast measurements. The large variability of the in-situ data results in large uncertainty estimates for the growth rate, exceeding those of Swirlab. Both measurements agree that the trend is lower, almost zero, at the end of the time series, although the uncertainty also increases. For the overlapping time period 2012-2018, the estimated growth rate is 6.3 ± 0.4 ppb/year (the uncertainty denoting the 95% confidence interval) for Swirlab and 6.8 ± 2.6 ppb/year for the mast measurements.  Surface plots of DLM trends at altitudes from 0 to 40 km based on two different methods: Swirlab retrieval using ground-based remote sensing and ACE-FTS retrieval using satellite-based solar occultation measurements. Swirlab trends are on the left and ACE-FTS trends on the right. The top row shows trends calculated for change over a one-year period, the middle row for a three-year range and the bottom row for a five-year range. The black line confines the areas where the trend differs from zero with a 95% confidence level.
At higher altitudes, the reference instrument is the ACE-FTS. As shown in Figure 13, the ACE-FTS and Swirlab time series agree well but the overall trends differ, although they are mostly within the uncertainty envelopes. As Swirlab is more constrained than ACE-FTS, the uncertainty estimates are larger for the latter. At 11.5 km, the eight-year trend from the beginning of 2010 until the end of 2017 (neglecting the edges of the DLM fit) shows a change of 4.9 ± 1.8 ppb/year in ACE-FTS measurements and 8.4 ± 1.2 ppb/year in Swirlab measurements. In the stratosphere, the ACE-FTS-derived growth rates are more subject to interannual variability. For the eight-year period of 2010-2017, overall, the trends at 21.5 km are 4.6 ± 4.8 ppb/year for ACE-FTS and 3.1 ± 3.4 ppb/year for Swirlab. At 26.5 km, the trends for the same period are 11.1 ± 11.0 ppb/year for ACE-FTS and −0.2 ± 4.4 ppb/year for Swirlab. The increasing growth rate towards the end of the ACE-FTS time series seems to be driven by the large increase from 2017 to 2018. This change is also noticeable in the Swirlab data albeit not as large as in the ACE-FTS data.
In Figure 14, the Swirlab growth rate is compared to that of the ACE-FTS at all levels and for different time scales. The ACE-FTS data are presented from 10 km upwards as the data are very sparse below that level. The analysis is also confined to the altitudes below 40 km as Swirlab is not very sensitive at altitudes above that and thus mostly follows the prior. Below 15 km, the instruments agree that there is a small but significant positive trend. However, between 15 and 20 km, ACE-FTS does not detect consistent growth while Swirlab does. From 20 to 35 km, Swirlab data show very stable methane concentrations with a weak yet insignificant negative trend until 2015 when the trend turns to weak positive, while ACE-FTS shows more variable trends in the beginning. In the latest years, Swirlab and ACE-FTS agree on a large but not statistically significant positive trend; however, in the ACE-FTS DLM fit, the trend is found significant in the three-and five-year time spans.
When comparing DLM trends for two datasets, it should be noted that the differences in concurrent measurements are not the only cause for differences in trends. The sampling of the data also has strong influence to the DLM fit. Random gaps in the data can create artifacts in the trends, especially if the measurement frequency is already sparse and the data gap hides a yearly maximum or minimum that is usually observed. Missing data in fall 2011 and early spring 2012 are partly the cause for the oscillations in the ACE-FTS stratospheric trends around winter 2011-2012. The larger amount of data in the Swirlab time series partially stabilizes the DLM fit but the missing values during the polar night make the fitting of the seasonal cycle challenging, as was also pointed out by Kivimäki et al. [50].

Discussion and Summary
We produced a ten-year time series of atmospheric methane profiles using a dimension reduction retrieval method to ground-based short wave infrared spectra measured with a Fourier Transform Spectrometer. After statistically correcting the solar zenith angle dependence, the retrieved profiles agree well with reference measurements: ACE-FTS above 10 km up to 40 km, AirCore from the surface to 25 km and 50-m mast air sampling measurements at the lowest altitude levels. The retrieval method succeeds in producing CH 4 profiles with realistic trends and seasonal variability even when using a single prior throughout the time series but the comparison to reference instruments performs better when a prior trend assumption and seasonal prior profiles are introduced.
The DLM approach in the paper results in realistic trend estimates. Missing data, especially in the ACE-FTS time series with less frequent measurements over Sodankylä, can induce short-term variability. The data gap due to the polar night mostly influences the seasonal cycle fit rather than the trend.
The tropospheric trends in the Swirlab time series agree with the trend in local 50-m mast measurements. In a wider perspective, the weaker tropospheric trends until 2011 and higher trends after that agree with other Arctic stations reported by Nisbet et al. [18]. Interestingly, the weaker trend in the end of the time series from 2017 to 2018 does not follow the global average CH 4 increase as reported by the Global Monitoring Division of NOAA's Earth System Research Laboratory [51]; however, this weaker trend is shown by both mast measurements and Swirlab. The slowdown of the growth rate could be interpreted as: (1) the decrease of surface fluxes in the influencing areas of the mast; (2) decreased biomass signals from e.g., Siberia; and (3) a transport pattern change. However, confirming the reason demands a thorough model study outside the scope of this work and thus can only be speculated here.
In the stratosphere, the CH 4 growth rates from ACE-FTS and Swirlab retrievals are generally in a good agreement. However, the confidence limits of the trend are larger than in the troposphere because of a larger uncertainty of the retrieved values, and a larger seasonal and interannual variability in the methane concentrations. Therefore, to be reliably detected from other variability, the growth rate would need to be consistent over a long time period. The significance of the estimated growth rate was evaluated by its deviation from zero with 95% confidence. As shown in Figure 14, both instruments detect mostly insignificant changes in the stratospheric CH 4 over Sodankylä. However, in the ACE-FTS time series, the positive stratospheric trends towards the end of the time series are found to be significant, contrary to the trends estimated for Swirlab. This could originate from higher vertical resolution of ACE-FTS but could also be caused by the differences in the DLM fitting due to less frequent data points in the ACE-FTS time series.
Albeit appearing towards the edge of the time series, the rising stratospheric trend from 2017 to 2018 draws attention. There is a slight suspicion that this could be an edge effect in the DLM fit, especially as there is a similar negative trend in the beginning of the time series. However, the late rise can be identified in both ACE-FTS and Swirlab data at the 26.5 km altitude as a high, sudden increase in CH 4 from 2017 to 2018, suggesting that the DLM-inferred positive trends are indeed described by the data. The stratospheric air at high latitudes mostly originates from the tropical troposphere. Therefore, the stratospheric methane concentrations in the Arctic are not expected to follow the changes in the local troposphere: the stratospheric air may have been transported several years after being in contact with methane sources at the surface. Hence, the recent sudden increase of CH 4 in the Arctic stratosphere could point to enhanced circulation from the tropics or to a decrease in stratospheric methane sink. However, continuing observations are needed to confirm and monitor this increase, and further studies involving atmospheric inverse modeling are necessary for pinpointing the reason for the observed growth rate.
Overall, the method discussed in this paper is a valuable addition to other methane profile measurements. While it does not reach the accuracy and vertical resolution of aircraft and balloon-borne measurements, the temporal coverage is much better. The information provided by the method will be a valuable addition to the existing continuous monitoring of the total column and the surface concentrations of methane. These continuous measurements can further help to develop models and reliably track the changes in the sources and sinks of methane. In this work, we considered a local scale but the method is publicly available and applicable to any ground-based, high-resolution spectrometer measuring in the short-wave infrared wavelength range.