The operational methane retrieval algorithm for TROPOMI

Abstract. This work presents the operational methane retrieval algorithm for the Sentinel 5 Precursor (S5P) satellite and its performance tested on realistic ensembles of simulated measurements. The target product is the column-averaged dry air volume mixing ratio of methane (XCH4), which will be retrieved simultaneously with scattering properties of the atmosphere. The algorithm attempts to fit spectra observed by the shortwave and near-infrared channels of the TROPOspheric Monitoring Instrument (TROPOMI) spectrometer aboard S5P. The sensitivity of the retrieval performance to atmospheric scattering properties, atmospheric input data and instrument calibration errors is evaluated. In addition, we investigate the effect of inhomogeneous slit illumination on the instrument spectral response function. Finally, we discuss the cloud filters to be used operationally and as backup. We show that the required accuracy and precision of


Introduction
Methane (CH 4 ) is the most important anthropogenic greenhouse gas after carbon dioxide (CO 2 ).While it occurs in smaller concentrations, it has a higher global warming potential per molecule than CO 2 .An accurate understanding of CH 4 sources and sinks is essential for a reliable prediction of climate change.Space-based measurements can provide continuous and global monitoring of CH 4 , leading to muchneeded improved constraints on the surface fluxes.
Several current and future satellite missions measure CH 4 abundances in the Earth's atmosphere.The SCanning Imaging Absorption spectroMeter for Atmospheric Cartogra-pHY (SCIAMACHY) aboard ENVISAT (Bovensmann et al., 1999) was the first space-based instrument to measure atmospheric CH 4 with sensitivity down to the Earth's surface.Frankenberg et al. (2005) were the first to use these measurements to constrain CH 4 surface fluxes.After loss of contact with ENVISAT in 2012, the Greenhouse gases Observing SATellite (GOSAT) is currently the only satellite measuring atmospheric CH 4 (Kuze et al., 2009).While GOSAT has a higher sensitivity and spatial resolution than SCIAMACHY, it has a fairly low spatial sampling.In 2017, the TROPOspheric Monitoring Instrument (TROPOMI) will be launched aboard the Sentinel 5 Precursor (S5P) satellite, and it will provide CH 4 measurements as one of its key products with unprecedented high precision, spatial resolution and global daily coverage.
The common goal of the above-mentioned missions is to provide atmospheric CH 4 concentrations with sufficient accuracy and spatiotemporal coverage to allow the assessment of CH 4 sources through inverse modelling.The observation strategy relies on measuring spectra of sunlight, backscattered by the Earth's surface and atmosphere, in the shortwave infrared (SWIR) spectral range.Absorption features of CH 4 molecules allow for retrieval of its atmospheric concentration with high sensitivity down to the Earth's surface where the main CH 4 sources are located.The applicability of such measurements for estimating source strengths, however, strongly depends on the precision and accuracy achieved.Residual systematic biases must be well below 1 % to facilitate inverse modelling (Meirink et al., 2006;Bergamaschi et al., 2007Bergamaschi et al., , 2009)).
Scattering by aerosols and cirrus is one of the major challenges for retrievals of CH 4 from space-based SWIR observations.While contamination by optically thick clouds can be filtered out reliably, optically thin scatterers are much harder to detect, yet they still modify the light path of the observed backscattered sunlight.This can lead to underestimation or overestimation of the true CH 4 column if not appropriately accounted for.The net light path effect strongly depends on the amount, size and height distribution of the scatterers as well as on the reflectance of the underlying surface (Aben et al., 2007;Gloudemans et al., 2008).Therefore, retrieval strategies rely on inferring the target gas concentration either simultaneously with atmospheric scattering properties or with a light path proxy.Frankenberg et al. (2005) introduced the "proxy" approach for CH 4 retrieval from SCIAMACHY measurements around 1600 nm, by using the simultaneously retrieved CO 2 column as a light path proxy.The proxy approach relies on the assumptions that scattering effects cancel in the ratio of the CH 4 column and the CO 2 column, and that a prior estimate of the CO 2 column is sufficiently accurate to recalculate the CH 4 column from the measured CH 4 / CO 2 ratio.In this case scattering is ignored in the forward modelling.Further applications of the proxy approach for CH 4 retrieval from SCIA-MACHY are described by Frankenberg et al. (2008) and Schneising et al. (2011).For GOSAT, the proxy approach has been successfully applied by Parker et al. (2011) and Schepers et al. (2012).
Alternatively, scattering-induced light path modification can be taken into account by simultaneously inferring the atmospheric CH 4 concentration and physical scattering properties of the atmosphere.Such "physics-based" methods have been developed for space-based CO 2 and CH 4 measurements from SCIAMACHY, GOSAT and the Orbiting Carbon Observatory (OCO); see e.g.Connor et al. (2008); Butz et al. (2009); Reuter et al. (2010); O'Dell et al. (2012).The physics-based methods make use of the oxygen (O 2 ) A-band in the near-infrared (NIR) around 760 nm and absorption bands of the target absorber in the SWIR spectral range.The advantage of physics-based methods for CH 4 retrieval compared to proxy methods is that they do not depend on prior information on the CO 2 column.On the other hand, the physics-based algorithms are more complex and computationally expensive.In addition, they may be limited by the information content of the measurement with respect to aerosol properties and related forward model errors in the description of aerosols.A detailed comparison between the two methods for GOSAT is provided by Schepers et al. (2012).
TROPOMI has four spectral channels in the ultraviolet (UV), visible (VIS), near-infrared (NIR) and shortwave infrared (SWIR), with spectral ranges of 270-320, 310-495, 675-775 and 2305-2385 nm, respectively.We use TROPOMI measurements in the NIR and SWIR for CH 4 retrievals.This spectral range does not allow for a light-pathproxy approach, and thus the effect of aerosols and cirrus needs to be accounted for by using a physics-based method as described above.The goal of this paper is to present the CH 4 retrieval algorithm for TROPOMI and investigate its sensitivity to algorithm assumptions, atmospheric input data and instrument calibration errors and filtering criteria.To this end, we simulated realistic TROPOMI measurements for aerosoland cirrus-loaded atmospheres under clear-sky and cloudy conditions.
The outline of this paper is as follows.We start with the methodology in Sect.2, giving an overview of the instrument, the retrieval algorithm and the methane data product.Then, a detailed error sensitivity study is presented in Sect.3, based on methane retrievals on a clear-sky global ensemble of simulated spectra.In Sect.4, a study of cloud filtering is performed on a TROPOMI orbit of simulated spectra that cover a realistic range of cloud parameters.Conclusions are presented in Sect. 5.

Methodology
The S5P satellite has a designed 7-year lifetime and will fly in a sun-synchronous orbit at 824 km altitude.Its single payload, TROPOMI, is a push-broom imaging spectrometer with a wide swath of 2600 km and a ground pixel of 7 × 7 km 2 in exact nadir.Approximately, TROPOMI observes a full swath per second, which results in ∼ 216 spectra per second.The instrument comprises two spectrometer modules, the first containing the UV, VIS and NIR spectral channels, and the second dedicated to the SWIR channel.We use the NIR and SWIR channels with spectral resolutions of 0.38 and 0.25 nm and spectral sampling ratios of 2.8 and 2.5, respectively (Veefkind et al., 2012).Since the NIR and SWIR detectors are incorporated in different instrument modules, the NIR spectra will be co-registered with the SWIR spectra before performing CH 4 retrievals.Examples of simulated TROPOMI NIR and SWIR spectra are shown in Fig. 1.

The retrieval algorithm
The S5P operational CH 4 retrieval algorithm is based on Re-moTeC, which was orginally developed for CO 2 and CH 4 retrievals from OCO and GOSAT obervations (Butz et al., 2009(Butz et al., , 2010)).A first performance study on TROPOMI mea- surements has been done by Butz et al. (2012), where a detailed description of the algorithm can be found.Here, we summarise the essentials and elaborate on the modifications made since then.

The forward model
The algorithm aims to infer the atmospheric state vector x from spectral measurements y in the NIR (757-774 nm) and SWIR (2305-2385 nm) ranges.This requires a forward model F that can accurately compute the measurement given the atmospheric state: where e y is the measurement noise error, and e F is the forward model error.The forward model incorporates the linearised radiative transfer model LINTRAN (Schepers et al., 2014).LINTRAN simulates the radiance at the top of the atmosphere I TOA on a fine internal spectral grid.The model needs a solar irradiance spectrum on the internal spectral grid as input, which is inferred from the daily solar measurements of TROPOMI using the deconvolution approach by van Deelen et al. (2007) and Wassmann et al. (2015).The simulated radiance measurement is obtained by a spectral convolution of I TOA with the instrument spectral response function (ISRF).
In addition to accuracy, an important requirement of the algorithm is computational speed.To be able to process TROPOMI's huge number of measurements, the CPU time per retrieval has to be on the order of seconds.Our algorithm achieves this by, among other things, using the linear kmethod for multiple scattering (Hasekamp and Butz, 2008).Single scattering is calculated line by line as it is computationally less expensive.To our knowledge, our algorithm is among the fastest in use for CH 4 full-physics retrievals, with an average CPU time of 7-10 s per retrieval (Hasekamp et al., 2016).
For a given model atmosphere, the forward model simulates spectra of backscattered sunlight, taking absorption and scattering by molecules and particles into account.The model atmosphere is defined for 36 pressure-equidistant vertical layers, with the top at 0.1 hPa and bottom at the surface pressure p surf .We calculate p surf by interpolating the meteorological surface pressure, from the European Centre for Medium-Range Weather Forecasts (ECMWF), on the surface elevation, from a digital elevation map (Danielson and Gesch, 2011;Farr et al., 2007).The absorbing trace gases of interest are O 2 in the NIR band, and CH 4 , H 2 O and CO in the SWIR band.The first guess layer subcolumns of these gases are calculated from input profiles of CH 4 and CO (from the global chemical transport model TM5; Houweling et al., 2014), and temperature, humidity and pressure profiles (from ECMWF forecast data).Molecular absorption features are calculated using appropriate spectroscopic databases (Tran et al., 2006;Rothman et al., 2009;Scheepmaker et al., 2012).Here, we evaluate the absorption cross sections on a 72-layer pressure-equidistant grid to account for the strong temperature dependence and pressure dependence of the cross sections.Molecular scattering properties are given by Rayleigh theory.Particulate absorption and scattering are computed with Mie theory using tabulated aerosol properties by Dubovik et al. (2006).
In our algorithm, the aerosol type is characterised by the refractive index and size distribution.The complex refractive index is fixed at n = 1.4 − 0.01i in the O 2 A-band and n = 1.47 − 0.008i in the SWIR band, which are averaged values derived from the ECHAM5-HAM model (Stier et al., 2005).

H. Hu et al.: TROPOMI methane retrieval algorithm
The size distribution is described by a power-law function with size parameter α (e.g.Mishchenko et al., 1999): where r 1 = 0.1 µm, r 2 = 10 µm, r is the particle radius and A is a normalisation constant.The amount of aerosol and its vertical distribution is provided by the vertically integrated column number density N aer and a normalised Gaussian height distribution h(z k ) with z k the height of layer k: where w = 2000 m, z aer is the central height and B is a normalisation constant.Thus, in layer k with thickness z k , the layer subcolumn n k is given by Butz et al. (2010) suggest that the exact choice of the fixedvalue parameters, such as the refractive indices or the width of the height distribution, does not affect the retrievals significantly.To support this claim, we have studied the effect of different values for the fixed-value aerosol parameters on our retrievals, and we arrive at the same conclusion; see Appendix A.
We have modified the retrieval forward model with respect to Butz et al. (2012) to account for chlorophyll fluorescence emission in a simplified manner as described by Frankenberg et al. (2012).In short, scattering of the fluorescence emission is ignored, and solely the spectral shape and absorption features by oxygen (O 2 ) are modelled.This allows fluorescence to be treated as a simple additive term to the radiance before convolution with the ISRF.The surface emission at the top of the atmosphere (TOA) is then modelled as where τ O 2 (λ) is the vertical optical thickness of oxygen, µ is the cosine of the viewing zenith angle and λ is the wavelength in nanometres.F surf s,755 and s represent the chlorophyll emission at 755 nm and its spectral slope, respectively.

The state vector
The state vector x consists of 25 elements; see Table 1.The target is a 12-layer pressure-equidistant vertical profile of CH 4 partial column number densities.The 12 retrieval layers are related to 36-layer model atmosphere by the shared interfaces; i.e. each retrieval layer is divided in three sublayers for the forward model calculations.Retrieved ancillary parameters are the total column number densities of H 2 O and CO, three scattering parameters N aer , α and z aer (related to amount, size and height), the surface albedo (constantterm and first-order spectral dependence1 and in the NIR and SWIR band), two terms to account for a spectral shift of the measurement (NIR and SWIR band) and two terms to account for chlorophyll fluorescence emission, F surf s,755 and s (see Eq. 3).An overview of the state vector elements and their a priori values is given in Table 1.We note that for the operational methane retrievals, the a priori value for F surf s,755 will come from a dedicated fluorescence retrieval scheme based on the Fraunhofer lines.In this study however, we simply took an a priori F surf s,755 = 0 photons s −1 m −2 nm −1 .Further, we would like to mention that our retrieved CO column should be regarded as a diagnostic for the CH 4 retrieval.The official TROPOMI CO product has a dedicated retrieval algorithm described by Landgraf et al. (2016).

The inversion procedure
The state vector is found by inverting Eq. ( 1), where the inverse method is based on a Philips-Tikhonov regularisation scheme (Phillips, 1962;Tikhonov, 1963).Regularisation is required because the inverse problem is ill-posed; i.e. the measurements y typically contain insufficient information to retrieve all state vector elements independently.Philips-Tikhonov regularisation aims to reduce contributions from measurement noise to the retrieved state vector while retaining valuable information.Because the forward model F(x) is non-linear in x, the inversion is performed iteratively by a step-size-controlled Gauss-Newton scheme, where at each iteration step the forward model is linearised.
The inverse algorithm finds x by minimising the cost function that is the sum of the least-squares cost function and a side constraint weighted by the regularisation parameter γ according to where S y is the diagonal measurement error covariance matrix, which contains the noise estimate.x a is an a priori state vector, and W is a diagonal weighting matrix that renders the side constraint dimensionless and ensures that only the CH 4 parameters and the scattering parameters contribute to its norm: W jj = 1/x a,j for the CH 4 column number densities and the three aerosol parameters, and W jj = 10 −7 /x a,j for all other state vector elements.The latter are thus retrieved in a least-squares sense.For determining γ , the L-curve criterion (Hansen, 1998) is applied in the baseline algorithm.
However, γ may be fine-tuned later using real observations.Table 1.State vector elements and their a priori values or references to source.

State vector element A priori value (units)
CH 4 subcolumns in 12 vertical layers TM5 (cm −2 ) CO total column TM5 (cm −2 ) H 2 O total column derived from ECMWF specific humidity (cm −2 ) Aerosol column N aer derived from aerosol optical thickness = 0.1 at 760 nm (cm −2 ) Aerosol size parameter α 3.5 (-) Aerosol height parameter z aer 5000 (m) Lambertian surface albedo in the NIR band maximum of measured reflectance in the NIR band (-) First-order spectral dependence surface albedo in the NIR band 0 (-) Lambertian surface albedo in the SWIR band maximum of measured reflectance in the SWIR band (-) First-order spectral dependence surface albedo in the SWIR band 0 (-) Spectral shift NIR 0 (nm) Spectral shift SWIR 0 (nm) Fluorescence emission at 755 nm F surf s,755 0 (photons s −1 m −2 nm −1 ) Fluorescence spectral slope s 0 (-)

The CH 4 data product
Although we retrieve methane in n = 12 sublayers, there is virtually no profile information in the measurement.The degree of freedom of signal of the retrieved methane profile is about 1.Therefore, the methane data product is given as a column-averaged dry air mixing ratio XCH 4 .This quantity is obtained from the methane entries of the retrieved state vector x i through where V air,dry is the dry air column (calculated from meteorological input surface pressure and the water vapour profile).
To interpret the retrieved XCH 4 correctly, one also needs the column averaging kernel A col that describes the sensitivity of the retrieved CH 4 column to changes to the true methane profile (see Rodgers, 2000, for details): As an example, in Fig. 2 we show the column averaging kernel corresponding to the methane retrieval performed on the simulated spectra of Fig. 1.We see that the column averaging kernel is around 1 in the lower atmosphere.From Eq. ( 5) it is clear that the closer A col is to 1, the more the retrieved column represents the true column (see also Eq. 6).This illustrates that the retrieval of methane columns from the SWIR has a nearly ideal sensitivity to methane in the troposphere and the tropospheric boundary layer.
For validation and interpretation purposes it is important to realise that the retrieved XCH 4 is related to the true methane profile x true and the a priori profile x a as where XCH 4,F is the bias caused by forward model errors, and XCH 4,y is the retrieval noise due to measurement noise.The standard deviation of the retrieval noise, i.e. the precision σ XCH 4 , follows from the error covariance matrix S x , which describes the effect of measurement noise on the retrieved parameters (see e.g.Rodgers, 2000, Sect.4.3): Together with XCH 4 and A col , the precision σ XCH 4 is given in the main data product.Note that for simulations, the bias XCH 4,F can be calculated from Eq. ( 6), since all other terms are either known from the simulations or calculated by the retrieval algorithm.In Sect.3, we evaluate XCH 4,F , XCH 4,y and their sensitivity to input errors.

Data filtering
Our algorithm has been designed to be efficient and accurate under certain assumptions; e.g. the atmosphere is cloud-free and plane-parallel with optically thin scatterers.Therefore, it is essential to filter out those scenes for which these assumptions break down to ensure data quality and reach the required accuracy.Moreover, if filtering can be done a priori, for example using cloud detection, we can considerably reduce the number of performed retrievals and hence computation time.Column averaging kernel for a typical methane retrieval.
The retrieval was performed using the simulated NIR and SWIR spectra in Fig. 1.

A priori filtering
For operational cloud filtering, measurements from the Visible Infrared Imaging Radiometer Suite (VIIRS) aboard the Suomi-NPP satellite will be used.S5P is foreseen to operate in loose formation with Suomi-NPP, meaning that both missions observe the same ground scene with a time delay of about 5 min.In case VIIRS data are not available, we have developed a backup cloud filter using H 2 O retrievals from the weak and strong absorption band, assuming a non-scattering atmosphere.Here we make use of the fact that clouds and aerosols will modify the optical path length in the two bands differently due to the different absorption strengths (Taylor et al., 2016;Frankenberg, 2014).Thus, there will be a difference between the H 2 O column retrieved from the strong band and the H 2 O column retrieved from the weak band in the presence of clouds.
Based on tests with simulations, we chose 2329-2334 nm as the weak absorption band and 2367-2377 nm as the strong absorption band (see also Scheepmaker et al., 2016).In Fig. 3, we show the absorption features of H 2 O in the SWIR range and highlighted the proposed weak and strong absorption band.Because the weak-band and the strong-band retrievals are differently affected by scattering, the ratio |H2Oweak−H2Ostrong| H 2 O strong is strongly correlated with cloud contamination.Scenes in which this ratio exceeds a certain threshold can be flagged as cloudy and filtered out.Based on a realistic ensemble of synthetic measurements (see Sect. 4), we find that a threshold of 0.08 filters out most cloudy scenes and keeps most of the clear-sky scenes.However, a direct cloud filter based on VIIRS data is shown to be superior to this indirect two-band cloud approach.Therefore, the two-band cloud filter will only be used as backup when VIIRS data are unavailable.Further, before performing any CH 4 retrievals, we filter out cases with solar zenith angle (SZA) larger than 70 • , and viewing zenith angle (VZA) larger than 50 • .These thresholds have been derived from simulations and shall be finetuned after launch using real observations.

A posteriori filtering
Butz et al. ( 2012) identified an a posteriori filter based on retrieved scattering parameters: SWIR aerosol optical thickness τ swir at 2350 nm, size parameter α and height parameter z aer .Following this work, we filter out retrievals with This indicates that the retrieval algorithm has difficulties with scenes that have optically thick scattering layers of large particles at high altitude.Additionally, we filter out cases with retrieved albedo in the SWIR band smaller than 0.02.More a posteriori filters (e.g. based on goodness of fit) will be determined after launch using real observations.

Sensitivity studies
Here, we evaluate the sensitivity of the retrieved XCH 4 forward model error to instrument errors and auxiliary input to the retrieval algorithm (e.g.meteorological data).To put the errors in perspective, the S5P product accuracy of XCH 4 was envisioned to be within 2 % and the product precision to be within 0.6 % (Veefkind et al., 2012).Accuracy is defined as the mean deviation from the truth, and precision is the variation due to random processes such as instrument noise.More recently, the requirements have been slightly reformulated as 1 % bias and 1 % precision (Hasekamp et al., 2016).From the 1 % bias, 0.8 % is reserved for forward model errors and 0.6 % for instrument-related errors.

Synthetic measurements: the global ensemble
We performed detailed sensitivity studies of the CH 4 algorithm on a global ensemble of simulated spectra consisting of land-only, clear-sky scenes on a 2.79 • × 2.8125 • latitude × longitude grid.This ensemble is, to a large extent, identical to the one used by Butz et al. (2012).It contains realistic aerosol-and cirrus-loaded scenes for 4 days, one per season.The treatment of aerosols and cirrus in the simulations is far more complex than in the retrieval forward model, where only one effective aerosol type is considered; see Sect.2.1.1.For the simulations, the aerosol physical properties and vertical distributions are derived from the global aerosol model ECHAM5-HAM (Stier et al., 2005) for five different chemical species and on a superposition of seven log-normal size distributions.The aerosol optical thickness is derived from MODIS observations (Remer et al., 2005).Furthermore, the simulations contains cirrus with optical thickness and vertical distribution based on CALIOP measurements (Winker et al., 2007).Finally, the surface albedo in the NIR is taken from the MODIS land albedo product in the 841-876 nm channel.For the albedo in the SWIR, the SCIA-MACHY surface albedo product at 2350 nm is used (Schrijver et al., 2009).For simplicity, we assume a constant surface albedo within the NIR and SWIR band in our simulations.We expect no difficulties to fit a more realistic spectral dependent surface albedo based on our experience with real GOSAT data.We refer to Figs. 2 and 3 of Butz et al. (2012) for the geographical distribution of the total optical thickness and surface albedo, respectively, used in the simulations.The measurements are simulated for the nadir viewing direction and a solar zenith angle that is representative of TROPOMI with an overpass time of 13:30 local time.We have 8633 simulated measurements in total.
While Butz et al. (2012) only investigated the scatteringinduced error, we increased the inconsistency between the simulation forward model and the retrieval forward model model.We have attempted to include the most important contributions to the forward model error, except for errors due to the underlying spectroscopic database, which have been investigated elsewhere; see Galli et al. (2012) and Checa-Garcia et al. (2015).The simulations were computed using a line-by-line radiative transfer model, whereas the retrieval method uses the linear k-method.In addition, the simulations have a higher vertical and spectral resolution than used in the retrieval.Furthermore, we have added chlorophyll fluores-cence emission.In the simulations, fluorescence is modelled to have a double Gaussian spectral shape (Guanter et al., 2010), which is different from the linear spectral shape assumed in the retrieval forward model (see Eq. 3).As in the retrieval scheme, we neglect scattering of fluorescence emission for simplicity.The fluorescence at the TOA then becomes For the parameters A 1 , A 2 , λ 1 , λ 2 , σ 1 and σ 2 , we use the same values as in Frankenberg et al. (2012).We only included fluorescence emission for scenes in the global ensemble with albedo NIR / albedo SWIR > 5 as a rough selection criterion for regions with vegetation.
After convolving the simulated TOA radiance with the ISRF, the spectra are superimposed with instrument noise from the TROPOMI noise model (Tol et al., 2011).For the NIR, the noise consists solely of shot noise, while for the SWIR, the noise is composed of both shot noise and a signalindependent term.The corresponding continuum signal-tonoise ratios are 500 in the NIR and 100 in the SWIR for a reference scene with surface albedo A s = 0.05, viewing zenith angle VZA = 0 • and solar zenith angle SZA = 70 • .
The baseline performance of the operational CH 4 algorithm is tested on the simulated global ensemble described above.In Fig. 4, we show a world map of the bias XCH 4,F before and after applying the a posteriori filters based on retrieved scattering parameters and albedo.Note that the error due to measurement noise XCH 4,y has been subtracted from the total XCH 4 error.This is a random error and is evaluated separately in Sect.3.3.1In Fig. 5, the cumulative probability distribution of the absolute XCH 4 retrieval error is shown (blue line).We get a convergence rate of 99 %, and 53 % of the converged retrievals pass the filters.Finally, 94 % of the valid retrievals have an absolute error < 1 %.In Fig. 5, we also plotted the cumulative probability distribution in case fluorescence emission is not fitted (red line).Then, we have a convergence rate of 95 and 92 % of the valid retrievals have an error < 1 %.Thus our retrieval results are improved by fitting fluorescence, which is why we have included this in the baseline.
We note that most of the rejected retrievals are found in the desert and dust regions in northern Africa and central Asia in April and July where our simulations capture dust storm events.For real TROPOMI measurements though, we will less likely filter out large regions because TROPOMI has higher spatial resolution than our simulations.Furthermore, our simulations only represent 1 day per season with a mean atmospheric state.Therefore our results merely indicate that certain regions will have fewer days with successful retrievals in certain seasons, not that entire regions will be rejected.

Sensitivity to atmospheric input data
We investigated the effect of imprecise atmospheric input from TM5 and ECMWF on the CH 4 retrievals.The results are summarised in Fig. 6.

A priori CH 4 profile
For the baseline performance test in Sect.3.1 we used the true profile of methane x true as the a priori profile for the retrievals.The bias XCH 4,F as defined in Eq. ( 6) should not depend on the choice of the a priori profile because this effect is accounted for by the averaging kernel.To test this, we take the zonal mean as the a priori profile, where we averaged over all longitudes within a 2.79 • latitude bin, with a column deviation up to ±2 %.To illustrate the effect on the global ensemble, we evaluate the root mean square (rms) of the XCH 4 bias (i.e. the total XCH 4 error minus the contribution due to noise) of all retrievals that pass the a posteriori filters.Figure 6 (first panel, blue line) shows that the a priori CH 4 does not influence the retrieval accuracy, in terms of the rms of the XCH 4 bias, nor the stability, in terms of the convergence rate or number of valid retrievals.

A priori H 2 O profile
To investigate the sensitivity to errors on the assumed H 2 O profile, we follow the same procedure as in Sect.3.2.1.The error on the prior H 2 O profile is established in the same way as for CH 4 , i.e. by taking a normalised zonal mean profile per latitude bin.Note that for H 2 O, there is an additional (minor) influence on the retrieval of XCH 4 through the dry air column.The H 2 O column error is varied up to ±10 %. Figure 6 (left panel, red line) shows that the prior H 2 O profile has neg-ligible influence on the rms of the XCH 4 bias, increasing it with < 0.01 %.There is a small effect on the convergence rate, reducing it from 99 to 97 % and leading to a reduction of valid retrievals from 53 to 50 %.We note that taking a zonal mean profile for H 2 O represents a worst case in terms of accuracy for the specific humidity of ECMWF.

Pressure
An erroneous pressure affects the retrieval of XCH 4 in two ways: first of all, through the pressure dependence of the cross sections and, secondly, through the retrieved air column that is used to convert the CH 4 total column to the dry air mixing ratio, XCH 4 .The latter will introduce a retrieval error of the same magnitude as the pressure error.To evaluate the net effect of a pressure error, the prior pressure profile is perturbed with a scaling factor up to ±0.3 %, corresponding to | P surf | ≈ 3 hPa.We expect a better accuracy from the ECMWF surface pressure together with the digital elevation map (Salstein et al., 2007;Danielson and Gesch, 2011;Farr et al., 2007).Figure 6 (right panel, blue line) shows that the increase in the rms of the XCH 4 bias is < 0.15 %.There is no visible effect on the stability of the algorithm.

Temperature
An error in the temperature will propagate to the XCH 4 retrievals though the temperature dependence of the cross sections.To investigate this effect, the temperature profile is offset up to ±2 K. Figure 6 (right panel, blue line) shows that the increase in the rms of the XCH4 bias is < 0.15 %.There is a small effect on the stability of the algorithm, reducing the convergence rate to 97 % and the number of valid retrievals to 47 %.

Sensitivity to instrument errors
We investigated the effect of different possible instrument and calibration errors on the CH 4 retrievals.The results are summarised in Table 2. Below we discuss each effect separately.

Signal-to-noise ratio
The simulated spectra include instrument noise as described in Sect.3.1.The precision is given by the standard deviation of the retrieval noise σ XCH 4 ; see Eq. ( 7).The world map in Fig. 7 shows the precision relative to the retrieved XCH 4 .Typically the precision is better than the accuracy.The signal-to-noise ratio only becomes a limiting factor for scenarios with snow-covered ground and large SZA, which is why we filter for SZA < 70 • and albedo > 0.02 to keep this error relatively small.

Instrument spectral response function
The synthetic measurements were created by convolving the underlying line-by-line spectra with a Gaussian ISRF with a full width half maximum (FWHM) of 0.38 nm in the NIR band and 0.25 nm in the SWIR band.The mission requirements state that the ISRF shall be known within 1 % of its maximum (Langen et al., 2011).This is achieved approximately by varying the FWHM by 1 %.Table 2 give the results for retrievals with an assumed error in the FWHM.We find that the retrievals are mostly sensitive to the accuracy of the ISRF in the SWIR band, leading to an increase of the rms XCH 4 bias from 0.56 % (baseline) to 0.67 %.We note that of all the instrument errors investigated here, the ISRF gives the largest error contribution to the methane retrievals.Thus, our study indicates that accurate calibration of the ISRF should have high priority.

Spectral calibration
According to the instrument requirements, the centre wavelengths of spectral channels are known within 2 pm (picometres).For our global ensemble, a spectral shift of 2 pm has negligible effect on the error characteristics.Here, we evaluate the effect of a wavelength shift of 1/20 of the spectral sampling distance, i.e. 10 and 5 pm for the NIR and SWIR band, respectively.The reference retrieval fits a spectral shift.
To test this fitting option, the synthetic spectra were shifted with a constant wavelength: where λ k is the real wavelength, λ k,meas is the measured wavelength at pixel k and λ the spectral shift.In Table 2, results are given for the case that a spectral shift is not fitted and is fitted (between brackets).When fitted, the performance is as good as for the reference retrieval, i.e. simulations with-out an error in the spectral position.This is as expected and indicates that the spectral shift fitting is robust.Optionally a wavelength-dependent shift, λ squeeze , can be fitted.To test this fitting option, the synthetic wavelength grid was "squeezed": where λ end and λ mid are the wavelengths at the end and middle of the band, respectively.Table 2 shows the performance of retrievals with this assumed error on the measured wavelength grid.Since a spectral squeeze is a second-order effect and has a much smaller impact on the XCH 4 retrievals than a spectral shift, it is not fitted in the baseline.However, our results show that, if needed, the option to fit a spectral squeeze can be used reliably.

Radiometric offset: additive factor
The effect of an unknown systematic offset in the Earth radiance is investigated.The offsets in the NIR and SWIR bands are independently varied with ±0.1 % of the continuum.Table 2 shows the effect of a radiometric offset to the XCH 4 retrievals.We note that a radiometric offset in the SWIR band causes a larger XCH 4 error than an offset in the NIR band.The latter is partly compensated by the retrieved fluorescence.

Radiometric gain: multiplicative factor
The absolute radiometric accuracy of the measurement of the Earth spectral radiance shall be better than 2 % according to the system requirements.To investigate the effect of such an error, the synthetic spectra were multiplied by a scaling factor G. Table 2 shows that there is negligible effect of an error of 2 % in radiometric gain.This error is largely compensated by the retrieved surface albedo.It follows that interaction be- Table 2. Effect of instrument calibration errors on convergence rate, fraction of valid retrievals after filtering and rms values of XCH 4 bias and precision for the global ensemble.Note that all sensitivities include the baseline forward model error, caused mainly by aerosol and cirrus scattering.The terms between brackets are for the cases where the relevant quantity is also retrieved.For each instrument calibration error, multiple simulation runs were performed with all combinations of errors in NIR and SWIR channels.The results shown here correspond to the runs with poorest performance in terms of the rms error.tween surface albedo and aerosols has a negligible impact for gain errors < 2 %.

Heterogenous slit illumination
For the two-dimensional TROPOMI push-broom spectrometer, light in across-slit dimension is dispersed by the instrument grating in order to spectrally resolve the received signal.The along-slit direction is aligned across flight direction to achieve the desired spatial resolution.For a homogeneous illumination of the instrument slit, the spectral instrument response is characterised extensively during the pre-flight calibration of the TROPOMI instrument, and it is used as baseline to simulate the TROPOMI radiometric measurement in our retrieval.In space, however, the instrument slit will be illuminated inhomogeneously due to ground scene heterogeneities on scales smaller than the instrument's field of view.Inhomogeneous illumination across the slit leads to a distortion of the ISRF as described by Noël et al. (2012), Caron et al. (2014) and Landgraf et al. (2016), and it can affect the retrieval accuracy of the TROPOMI methane product.Small-scale heterogeneities of the ground scene are generally caused by spatial variations of surface reflection and by broken clouds.Because of our strict cloud filtering, spatial variations in surface reflection are the only cause of methane retrieval biases due to inhomogeneous slit illumination.To evaluate the effect of surface scene heterogeneity on our XCH 4 product, we employ the instrument model described by Landgraf et al. (2016) for both the NIR and SWIR bands of our retrieval.Furthermore, we use a high spatial resolution MODIS albedo map for the 50 × 50 km 2 marsh region in central Siberia with structures in the surface reflection due to ponds, shown in Fig. 8. Depending on the scene heterogeneity in the flight direction, the XCH 4 error shows an oscillation structure with a maximum amplitude ≤ 0.4 %, a standard deviation of 0.12 % and a mean error of −0.01 %.For a particular temporal and spatial sampling of the scene, a pseudo-random scatter is introduced to the XCH 4 product.This means that overall the effect can be considered small.One may consider this error as a limitation when interpreting very localised sources in surroundings of heterogenous surface reflection, but for most applications some averaging either in time or space will be done, which reduces this error.

Cloud filtering
The global ensemble from Butz et al. (2012) as used in Sect. 3 cannot be used to evaluate the cloud filters, because it consists purely of cloud-free scenes.Therefore, the cloud filters are tested using a new ensemble of synthetic measurements, which also include scenes with water clouds.

Synthetic measurements: the TROPOMI test orbit
To test the performance of the proposed backup cloud filter, we have constructed synthetic TROPOMI L1B radiance spectra for an entire orbit that passes over Africa.We used realistic viewing geometries from a TROPOMI orbit simulator provided by the Royal Netherlands Meteorological Institute (KNMI).The meteorological data are from ECMWF.We used CO and CH 4 model profiles from TM5 (Houweling et al., 2014).The surface elevation comes from a digital elevation map constructed by KNMI based on USGS (Danielson and Gesch, 2011) and NASA data (Farr et al., 2007).The aerosol and cirrus properties and surface albedos are taken from the global ensemble, which means that these are the same for all TROPOMI ground pixels within a ∼ 3 • × 3 • latitude × longitude box.While the global ensemble used for the sensitivity studies is land-only and clear-sky, the test orbit contains cloudy scenes over land and ocean.Cloud fraction, cloud optical thickness and cloud top pressure are obtained from MODIS Aqua measurements at 5 × 5 km 2 .Here, the cloud top height is derived from the cloud top pressure and  the surface pressure, also provided by MODIS.The cloud properties are collocated in time and space to the TROPOMI orbit and used in the measurement simulation.For fractional clouds, the independent-pixel approximation is used to combine the cloudy and clear-sky parts of the scene.For reference, the cloud fraction used in the simulations is shown in Fig. 9, upper left.We note that 28 % of the TROPOMI ground pixels are fully cloud-free for this test orbit.Globally, one would expect, on average, 20 % cloud-free pixels (Krijger et al., 2007).Resampling of the auxiliary data on the TROPOMI ground pixels has been performed using the Multi-Instrument Preprocessor (MIPrep) developed at SRON.

Performance of cloud filters
First, we show the performance of the XCH 4 retrievals on the test orbit when no a posteriori data filtering is applied, only a priori filtering of ocean pixels and SZA > 70  cloud-contaminated measurements lead to large XCH 4 errors (> 2 %).
Assuming that we have cloud data from VIIRS, we would then be able to filter out the cloudy pixels almost perfectly.To illustrate the effect, we filtered out pixels with cloud fraction > 0.02 in Fig. 9, lower left.Note that in this case we have also applied a posteriori filtering based on retrieved scattering parameters and albedo.One is then left with valid retrievals of ∼ 3 % of all simulations in the test orbit.
In comparison, the performance of the backup cloud filter based on the difference between the H 2 O column retrieved from strong and weak bands (see Sect. 2.3) is shown in Fig. 9, lower right.The backup cloud filter removes most cloudy pixels, but some remain.In Table 3 and Fig. 10, the statistics of XCH 4 retrievals on the orbit are summarised.After cloud filtering with MODIS data (representative of operational VI-IRS data), the results for the test orbit are comparable to the clear-sky global ensemble.However, the backup cloud filter is less effective.The rms of the XCH 4 bias is then 0.71 % instead of the 0.56 % that is expected for the operational VIIRS cloud mask.

Conclusions
This paper describes the algorithm baseline of the operational methane retrievals from the S5P measurements.The level 2 product includes the column-averaged dry air mixing ratio XCH 4 , the column averaging kernel and the noise standard deviation.In order to account for the effect of aerosols and cirrus, the algorithm developed retrieves the methane column simultaneously with effective scattering parameters related to particle amount, size and height distribution.The choice of scattering parameters reflects the information content of the measurements as closely as possible.The retrieval algorithm uses the radiance and irradiance measurements in the SWIR 2305-2385 nm band and additionally in the NIR band between 757 and 774 nm (O 2 A-band).The forward model of the retrieval algorithm uses online radiative transfer calculations, fully including multiple scattering in an efficient manner.Absorption cross sections of the relevant atmospheric trace gases and optical properties of aerosols are calculated from lookup tables.The inversion is performed using Phillips-Tikhonov regularisation in combination with a reduced-step-size Gauss-Newton iteration scheme.
To test the developed algorithm we generated two ensembles of simulated measurements that cover the range of scenes that will likely be encountered by the S5P instrument: one clear-sky global ensemble and one test orbit containing cloud-contaminated measurements.Overall, the algorithm developed performs well in correcting for the effect of aerosols and cirrus clouds on the retrieved XCH 4 .For both ensembles, ∼ 80 % of the cases have an XCH 4 error < 0.5 % and ∼ 95 % have an error < 1 %.To achieve this, it is necessary to apply a priori filtering of cloud-contaminated scenes and a posteriori filtering based on retrieved parameters.We illustrated the performance of the proposed backup cloud filter based on retrievals of H 2 O from weak and strong absorption bands in the SWIR under the assumption of a nonscattering atmosphere.It should be noted that the cloud filter based on S5P measurements itself is less efficient than the VIIRS cloud mask for water clouds.
Apart from forward model errors induced by aerosols, we also studied effects of model errors in temperature, pressure and water vapour profiles.We expect to stay within product requirements for errors in input profiles of water, pressure and temperature below 10 %, 0.3 % and 2 K, respectively.Another relevant source of errors to the CH 4 data product could be spectroscopic errors.This has been studied in detail by Galli et al. (2012) andCheca-Garcia et al. (2015).Note that a study is ongoing to improve the spectroscopic data for TROPOMI's SWIR spectral range (Loos et al., 2015).Concerning instrument errors, we found that the most critical error source is an error in the ISRF in the SWIR band.To conclude, we have shown that for a compliant instrument, our algorithm provides a methane product that meets the requirements.

Data availability
The underlying research data for the simulated atmospheric profiles and TROPOMI measurement simulation and retrievals are available upon request from Haili Hu (h.hu@sron.nl).HITRAN spectroscopic line parameters (Rothman et al., 2009) are available through HITRAN online (http://hitran.org),and the line parameters from Scheepmaker et al. (2012) are available in their Supplement.Appendix A: Sensitivity of retrievals to fixed-value aerosol parameters Our retrieval algorithm retrieves three effective aerosol parameters, namely total column, size parameter and central height.Other aerosol parameters are assumed to be fixed, such as the refractive indices n and width w of the aerosol layer (i.e.FWHM of Gaussian height distribution).The values in the baseline retrieval are fixed at n NIR = 1.4 − 0.01i, n SWIR = 1.47 − 0.008i and w = 2 km, respectively.To justify this assumption, we studied the effect of the fixed-value aerosol parameters on the methane retrievals for our synthetic global ensemble as described Sect.3.1.We have performed four retrieval test runs with refractive indices of 1.37 − 0.0005i and 1.55 − 0.02i, representative of the majority of aerosol types (Dubovik et al., 2002), and with widths of 1 and 3 km.For each test run the other values are kept to the baseline values.For simplicity, we used the same refractive index for NIR and SWIR in the test runs.The results are summarised in Table A1 from which we conclude that, within a reasonable range, the exact value of the fixed aerosol parameters does not affect the methane retrievals significantly.We do see a slight deterioration of our results when a relatively low refractive index is used.

Figure 1 .
Figure 1.Simulated TROPOMI spectra in the NIR and SWIR channels for a scenario from the global ensemble; see Sect.3.1.The scenario, located in a rural area in Utah (USA), is observed at nadir and SZA = 50 • .

Figure 2 .
Figure 2. Column averaging kernel for a typical methane retrieval.The retrieval was performed using the simulated NIR and SWIR spectra in Fig. 1.

Figure 3 .
Figure 3. Simulated reflectance spectrum showing the absorption features of H 2 O in the SWIR spectral range for the same scenario as in Fig. 1.The red shaded area represents the weak absorption band, while the blue area represents the strong absorption band of H 2 O.The difference between the water column retrieved from the weak and strong band is used for the backup cloud filter.

Figure 4 .Figure 5 .
Figure 4. Residual forward model XCH 4 errors for the baseline retrieval method for the clear-sky global TROPOMI measurement ensemble before (a) and after (b) a posteriori filtering.The ensemble covers scenes for a day in January (JAN), April (APR), July (JUL) and October (OCT) as described in Sect.3.1.

Figure 6 .
Figure6.Influence of errors in atmospheric input on the accuracy (panels 1 and 2) and stability (panels 3 and 4) of XCH 4 retrievals.The upper and lower x axes refer to perturbations.Accuracy refers to the XCH 4 bias, and stability refers to the fraction of converged and valid retrievals.The profiles of methane and water have been perturbed in panels 1 and 3; the pressure and temperature profiles have been perturbed in panels 2 and 4.

Figure 7 .
Figure 7. Relative precision of XCH 4 due to the instrument noise for the a posteriori filtered data set for the clear-sky ensemble in Fig. 4.

Figure 8 .
Figure 8. Methane bias due to heterogenous slit illumination for spatially varying surface reflection over a marsh scene at Siberia close to the river Ob at latitude 62.8 • N and longitude 72.1 • E. Measurement simulations are performed with the instrument model by Landgraf et al. (2016) for an instantaneous field of view of 3.4 km across the slit and 7.0 km along the slit.

Figure 9 .
Figure 9. Simulations of a TROPOMI orbit.Panel (a) shows MODIS cloud fraction resampled on the orbit's ground pixels.Panel (b) gives the XCH 4 bias of all processed pixels that converged.Panel (c) and (d) give XCH 4 bias of the valid retrievals after cloud filtering with MODIS data and the backup cloud filter, respectively.Here, we also applied the a posteriori filter based on retrieved scattering parameters and albedo.

Figure 10 .
Figure 10.Cumulative probability distribution of the absolute XCH 4 bias for the simulated level 1b orbit.Cloud-contaminated measurements are filtered using MODIS data (blue line) or the backup cloud filter (red line).

Table 3 .
Error statistics of XCH 4 retrievals from the L1B orbit using different cloud filters in comparison to the global clear-sky orbit.The performance of the MODIS filter is expected to be comparable with the operational cloud filter using VIIRS data.

Table A1 .
Effect of fixed-value aerosol parameters on convergence rate, fraction of valid retrievals after filtering and rms values of XCH 4 bias and precision for the global ensemble.See text for values used in the baseline run.Convergence Valid retrievals rms of XCH 4 bias rms of XCH 4 precision