The radiance and state of polarization of light at a given wavelength can be described by an intensity vector I which has the Stokes parameters I, Q, U, and V as its components (Chandrasekhar, 1960): I = I , Q , U , V T . (5)

Here, T indicates the transposed vector, and the Stokes parameters are defined with respect to a certain reference plane. We will use the local meridian plane as reference plane.

The DoLP can be computed from the Stokes parameters as: D o L P = Q 2 + U 2 I . (6)

So, in general, the forward model should simulate the Stokes parameters I, Q, and U for different spectral bands/pixels and viewing angles. The Stokes parameters are calculated from the optical properties by the SRON radiative transfer model LINTRAN v2.0 (Schepers et al., 2014). For a given spectral measurement in band/pixel i, the forward model should simulate E ̄ λ i = ∫ λ start λ end s e a r t h , i λ E λ d λ , (7)

where E represents either I, Q, or U, s _earth is the instrumental spectral response function (ISRF), and λ _i is the wavelength assigned to spectral measurement i. Application of Eq. 7 requires E(λ) to be computed at high spectral resolution. To avoid time-consuming line-by-line calculations, we use the linear-k method described by Hasekamp and Butz (2008).

If the dependence of E on λ is weak within the spectral range of the ISRF (e.g., at continuum wavelengths or for spectrally smooth absorbing trace gases), we can approximate Eq. 7 by E ̄ λ i = E λ i . (8)

In any case (Eq. 7 or Eq. 8), we need a radiative transfer (RT) model to compute I(λ), Q(λ), and U(λ) for a given model atmosphere.

For the forward model calculations, the model atmosphere is discretized into a number N _RT of homogeneous vertical layers, provided as input at a fixed number of altitude levels z _lev above the surface of the atmosphere. In the standard setup, we use 15 vertical layers (2 km spacing between 0 and 20 km, 4 km spacing between 20 and 36 km, and 1 layer above 36 km). First, we determine the pressure at these altitude levels by linear interpolation of the input pressure profile p _met as a function of altitude z _met from a meteorological data source (e.g., from the Copernicus Atmosphere Monitoring Services, CAMS), to the forward model levels z _lev + z _surf, where z _surf is the surface height from a digital elevation map (DEM).

Each layer k in the model atmosphere is characterized by the scattering optical thickness τ _sca,k, absorption optical thickness τ _abs,k, and scattering phase matrix P _k . τ a b s , k = τ a b s , m o l , k + τ a b s , a e r , k + τ a b s , c i r , k , (9) τ s c a , k = τ s c a , m o l , k + τ s c a , a e r , k + τ s c a , c i r , k , (10) P k Θ = τ s c a , m o l , k P m o l Θ + τ s c a , a e r , k P a e r Θ + τ s c a , c i r , k P c i r Θ τ s c a , k , (11) where the underscripts “mol,” “aer,” and “cir” refer to the contributions from molecules, aerosols, and cirrus, respectively, and Θ denotes scattering angle.

We compute τ _sca,mol,k using the Rayleigh cross-sections from Bucholtz (1995) and the sub-column of air in layer k. For the Rayleigh phase matrix P _mol, we refer to Hansen and Travis (1974). The molecular absorption optical depth is first calculated per absorber (j) and per discretized atmospheric layer (k). To account for the temperature and pressure dependence of the molecular absorptions within a model layer, we divide each model layer k in N _sub,k sublayers, such that the pressure thickness Δ P i sub of the sublayer < 10 hPa. The molecular absorption optical depth for layer k and absorber j is given by: τ abs,k, j λ = ∑ i sub = 1 N sub, k σ j p i sub , T i sub , λ N gas , k , j Δ p i sub Δ p k , (12) with the absorption cross-section σ j ( p i sub , T i sub , λ ) of molecule j at wavelength λ, pressure p i sub , and temperature T i sub at the center of sublayer i _sub of atmospheric layer k. Δ p i sub is the pressure thickness of sublayer i _sub and Δp _k of model layer k. Furthermore, N _gas,k,j is the sub-column of absorber j in layer k, given by N gas , k , j = r gas , k , j N a i r , k , (13) where r _gas,k,j is the dry-air mixing ratio of gas j in layer k, and N _air,k is the dry-air sub-column in layer k. The total molecular absorption optical depth for model layer k is then given by τ a b s , m o l , k λ = ∑ j τ abs,k, j λ . (14)

The aerosol optical properties τ _sca,aer,k, τ _abs,aer,k, and P _aer in layer k are computed from their size distribution, refractive index, and shape. Herein, we assume that the aerosol size distributions consist of a number of size modes, where each mode, in principle, has its own refractive index, vertical distribution, and particle shape. The optical properties are computed per mode using the Mie- and T-matrix improved geometrical optics database by Dubovik et al. (2006) along with their proposed spheroid aspect ratio distribution for computing optical properties for a mixture of spheroids and spheres.

The size distribution n _aer(r) of each aerosol mode is either described by a log-normal or power-law function. The log-normal size distribution is written as: n a e r r = 1 2 π σ g r exp − ln ⁡ r − ln r g 2 / 2 σ g 2 , (15) where r is the radius (or radius of a volume equivalent sphere), r _g is the median radius, and σ _g is the standard deviation. Instead of r _g and σ _g , we use the effective radius r _eff and effective variance v _eff because they are less dependent on the actual shape of the size distribution (Hansen and Travis, 1974). r _eff and v _eff are related to r _g and σ _g by: σ g 2 = ln 1 + v eff (16) and r g = r eff / 1 + v eff 5 / 2 . (17)

The power-law size distribution is given by n a e r r = C for  r ≤ r 1 C r r 1 − p for  r 1 < r ≤ r 2 0 for  r > r 2 . (18) The cut-offs are r ₁ =0.1 μm and r ₂ = 10 μm, and the constant C is determined from normalization of the size distribution.

The refractive index of aerosols is described as a function of the MAP spectral sampling points and the central wavelengths for each spectrometer band but is assumed constant inside a spectrometer window. Each aerosol mode is assumed to be a mixture of a set of components, such as inorganic component (INORG), organic carbon (OC), black carbon (BC), dust (DU), and water, where the spectral dependence of the real part and imaginary part of each component are predefined based on the standard types of D’Almeida et al. (1991). Therefore, the complex refractive index as a function of wavelength for a given mode is calculated as follows: m λ = ∑ k = 1 n c c k m k λ , (19) where m _k (λ) represents prescribed refractive indices of the kth aerosol component and c _k is the coefficient of the component.

A normalized Gaussian altitude distribution with a center height z _aer and a full width half maximum (FWHM) w ₀ is used to calculate the number sub-column of aerosols N _{aer, k, i} for layer k and mode i from the total column of aerosol particles N _{aer, i} for mode i. Hence, for model layer k with a layer height z _k and a thickness Δz _k , the number sub-column yields: N aer, k, i = N aer, i B ⁡ exp − 4 ⁡ ln ⁡ 2 z k − z a e r 2 w 0 2 Δ z k , (20) where N _{aer, i} is the aerosol column number for mode i and B is the normalization constant of the Gaussian.

For the setup of the MAP-CO2I retrieval, we use a parametric 3-mode aerosol description, where the size distribution is described by 3 log-normal modes (i.e., N _modes = 3), with one fine mode and 2 coarse modes (soluble and insoluble). For the fine mode, the state vector includes r _eff, v _eff, N _aer, f _sph, and the refractive index coefficients c _k that correspond to the standard refractive index spectra of inorganic aerosol (real part), black carbon (imaginary part), organic carbon (imaginary part), and water. The coarse insoluble mode consists of non-spherical dust. For this mode, the state vector includes r _eff, N _aer, and a coefficient for the imaginary part of the dust refractive index. The fixed parameters are f _sph = 0, v _eff = 0.6, and c _k = 1 for the real part dust refractive index. One value for z _aer is included which is assumed to be the same for modes 1 and 2, while w ₀ is fixed at 2,000 m. The third mode is a coarse soluble mode. For this mode, the state vector includes r _eff, N _aer, and a coefficient c _k of the inorganic refractive index spectrum. The fixed parameters are f _sph = 1, v _eff = 0.6, and z _aer = 500 m. The aerosol parameters in the state vector for MAP-CO2I retrieval are shown in Table 4.

The aerosol size distribution for CO2I-only retrieval is described by a single power-law mode which reflects the information content of the measurements (Butz et al., 2009). Here, N _aer and power p are included in the state vector. The refractive index is fixed by 1.4–0.003i for all spectral bands. Furthermore, z _aer is included in the state vector and w ₀ is fixed at 2,000 m.

Optical cirrus properties are computed using the parameterization scheme proposed by van Diedenhoven et al. (2012). It provides extinction cross-section σ _e , single scattering albedo ω, and scattering phase matrix which are tabulated in a look-up table that stores the optical properties as a function of effective radius r _eff, hexagonal prism aspect ratio α, degree of crystal distortion δ, and wavelength λ. Linear interpolation is used to obtain optical properties at the actual values of r _eff, α, δ, and λ. The model uses a description of ice crystals as a mixture of plates with α < 1 as proxies for complex ice crystals. Crystal distortion parameter δ parameterizes the optical effects of crystal complexity from microscales to macroscales (Macke et al., 1996; van Diedenhoven et al., 2014). A maximum of δ = 0.7 is used. The current look-up table contains optical properties at λ =1,589, 1,884, and 2,117 nm. Interpolation leads to moderate positive and negative absorption biases in SWIR-1 and SWIR-2 bands, respectively, which are not expected to substantially affect our sensitivity study. For the vertical distribution of cirrus number column N _cir, Gaussian distribution is also adopted with center height z _cir and FWHM w ₀. The state vector includes r _eff, α, δ, N _cir, and z _cir for cirrus if cirrus is fitted in the retrieval, and w ₀ is fixed to 1,000 m.

The sub-columns N _gas,k,j of the j different trace gases in layer k of the model atmosphere are given by N gas , k , j = s j N std , k , j , (21) where N _std,k,j is the sub-column for trace gas j corresponding to an a priori vertical profile and s _j is the corresponding scaling factor independent on k. The O₂ vertical profile follows the pressure profile. For MAP-only retrievals, we fix the scaling factors s _j to 1.0, that is, their atmospheric abundance is fixed to the a priori value throughout the retrieval, while for MAP-CO2I and CO2I-only retrievals, we include the scaling factors s _j corresponding to CO₂, CH₄, and H₂O in the retrieval state vector.

In the inversion procedure, we invert the linearized forward model of Eq. 3 for iteration step n to find the state vector x _n+1 for iteration step n + 1. Hereto, we minimize the following cost function (Tikhonov, 1963): x n + 1 = arg min x K x − y T S y − 1 K x − y + x − x a T γ 2 H − 1 x − x a , (27) which we transform to x ̃ n + 1 = arg min x ̃ K ̃ x ̃ − y ̃ T K ̃ x ̃ − y ̃ + γ 2 x ̃ − x ̃ a T x ̃ − x ̃ a , (28) where K ̃ = S y − 1 2 K H 1 2 , x ̃ = H − 1 2 x , x ̃ a = H − 1 2 x a , and y ̃ = S y − 1 2 ( y − F ( x n ) ) . x _a is the a priori state vector, S _y is the measurement error covariance matrix, γ is a regularization parameter, and H is a regularization matrix (assumed diagonal) which ensures that all state vector parameters of the same order of magnitude and determines the relative weight of parameters in the side constraint (Hasekamp et al., 2011). Note that if γ = 1 and H is the a priori error covariance matrix, Eq. 27 (and hence Eq. 28) reduces to the cost function of the optimal estimation method (Rodgers, 2000).

The solution of Eq. 27 is given by: x ̃ n + 1 = x ̃ n + Λ K ̃ T K ̃ + γ 2 I − 1 K ̃ T y ̃ − γ 2 x ̃ n − x ̃ a . (29) Λ is a filter/damping factor between 0 and 1, which limits the step size for each iteration of the state vector. In this way, we use a Gauss–Newton scheme with reduced step size to avoid diverging retrievals. At each iteration, an optimal combination of γ and Λ is chosen that provides the best match between a simplified (fast) forward model and the measurements via χ ² assessment. χ ² is a measure of goodness of fit given by Eq. 30. χ 2 = 1 N meas ∑ i = 1 N meas y i − F i e i 2 , (30) where N _meas is the total number of measurements, y _i is the measurements, F _i is the forward model simulations with the retrieved parameters, and e _i is the measurement uncertainties.

The error covariance matrix S _x of the retrieved state vector is given by S x = S r + S e , (31) where S _r is the regularization error covariance matrix which describes the effect of the a priori error covariance matrix S _a on x, S r = I − A S a I − A T . (32) S _a takes the value of H in this study. S _e is the retrieval error covariance matrix that describes the effect of measurement and forward model errors on x, S e = D S y D T , (33) where D is the contribution or gain matrix D = K T S y − 1 K + γ 2 H − 1 − 1 K T S y − 1 , (34) and A is the averaging kernel A = D K . (35)

The iteration starts with a first guess of the state vector. The first guess of aerosol parameters and surface parameters for MAP bands are obtained using a multimode look-up table (LUT) if MAP measurements are used in the retrieval. Otherwise, the first guess of aerosol parameters is the a priori values. Details of the LUT retrieval are described by Fu and Hasekamp (2018). We use tropopause height (calculated from CAMS data) as the first guess of cirrus layer height if cirrus is retrieved. For the first guess of the BRDF scaling factor at each CO2I band, we use the approximation of a Lambertian reflector for an extinction-free atmosphere: A = I λ cont / F 0 λ cont π μ 0 . (36) Here, the first guess A is constant over one spectral window. λ _cont is the wavelength that gives the maximum reflectance I(λ)/F ₀(λ) over that spectral band. F ₀ is the incoming solar irradiance.

The full-physics retrieval process is demonstrated by the flowchart in Figure 1. The a priori values and weights of the state vector for MAP-CO2I and CO2I-only retrievals are given in Table 4.

The non-scattering retrieval uses a simplified forward model, ignoring scattering in the atmosphere, and also a simplified inversion procedure. The non-scattering forward model calculates the reflected top-of-atmosphere radiance I as a function of wavelength λ by neglecting scattering: I λ = F 0 λ A λ μ in π exp − τ abs λ / μ ̃ , (37) with the air mass factor: 1 μ ̃ = 1 μ in + 1 μ out , (38) where A is the Lambertian surface albedo with its spectral dependence in a polynomial form parameterized by Eq. 25. τ _abs is the absorption optical thickness integrated over vertical layers and molecular species. In addition to the scaling factors for O₂, H₂O, CO₂, and CH₄, the linear spectral dependence for Lambertian surface albedo and the wavelength shifts are included in the state vector.

Similar to the full-physics retrieval, the non-scattering inversion is performed iteratively. For each iteration step, we minimize a least-square cost function: x n + 1 − x n = arg min x K x − d n T S y − 1 K x − d n , (39) where d _n = y − F(x _n ). Since the non-scattering forward model is sufficiently linear in the state vector x, the cost function can be solved directly without step size reduction by x n + 1 = x n + K T S y − 1 K − 1 K T S y − 1 d n . (40)

The non-scattering retrieval is performed to each spectrometer window individually to obtain the total column of the target molecular species of that band according to Table 2. The state vector definition for the non-scattering retrieval is shown by Table 5.

The proxy methane product uses the non-scattering output and the a priori estimate of XCO₂ ( X C O 2 apr ) to construct the estimate of XCH₄ proxy product ( X C H 4 proxy ) by X C H 4 proxy = C CH 4 C CO 2 X C O 2 apr , (41) where total columns C _CH4 and C _CO2 are retrieved under the non-scattering assumption from the measurements in the SWIR-1 band. Similarly, the proxy XCO₂ retrievals can be derived as follows: X C O 2 proxy = C CO 2 C CH 4 X C H 4 apr . (42) The a priori estimates X C O 2 apr and X C H 4 apr should be sufficiently close to the true column-average dry-air mole fractions.

For the filtering of cirrus clouds, we will use the result from non-scattering retrievals in the different spectral bands (NIR, SWIR-1, and SWIR-2) to filter for clouds, following Taylor et al. (2016). Here, we retrieve independent estimates of the CO₂ and H₂O column abundances using observations taken at 1.61 μm (SWIR-1 weak CO₂ band) and 2.06 μm (SWIR-2 strong CO₂ band), while neglecting atmospheric scattering (see above). The CO₂ and H₂O column abundances retrieved in these two spectral regions differ significantly in the presence of cloud and scattering aerosols. Ratios of the retrieved CO₂ (R _CO2) and H₂O (R _H2O) column abundances are computed as follows: R gas = C gas W C gas S , (43) where C gas W and C gas S represent the vertical column density of the retrieved gas (CO₂ or H₂O) in the weak and strong absorption bands, respectively. For CO2M cirrus filtering, we use R CO 2 = C SWIR 1 CO 2 / C SWIR 2 CO 2 and R H 2 O = C SWIR 1 H 2 O / C SWIR 2 H 2 O , where C SWIR 1 CO 2 and C SWIR 2 CO 2 are the columns of CO₂ retrieved at SWIR-1 and SWIR-2 bands, and C SWIR 1 H 2 O and C SWIR 2 H 2 O are the columns of H₂O retrieved at SWIR-1 and SWIR-2 bands, respectively.

Cirrus clouds and aerosols modify the optical path length in the two bands differently, producing column abundance ratios significantly different from unity (Taylor et al., 2016). There are two fundamental reasons why the ratio deviates from unity in the presence of scattering. First, for most terrestrial surfaces, albedos in the 1.6 μm band are most often higher than at 2.0 μm. This yields a variable fractional contribution of scattered light to the measured radiances. Second, the 1.6 and 2.06 μm band strengths are very different, resulting in different sensitivities to atmospheric scattering. If no scattering is assumed in the retrieval, a deviation from unity in the ratio thus indicates a substantial variation in the photon path length (PPL) distribution between the two bands, while in the absence of scattering, this ratio approaches unity.

The upper panel of Figure 4 illustrates ratios R _CO2 and R _H2O as a function of true COT values, showing that most cirrus-free pixels (with COT = 0) have ratio values close to 1, that is, R _CO2 < 1.04 and R _H2O < 1.06. Moreover, for cirrus-contaminated pixels, R _CO2 and R _H2O have a positive correlation with COT, whereas the correlation is stronger for R _CO2. This suggests the ratios can be applied for the filtering of cirrus.

To investigate the capability of R _CO2 and R _H2O for cirrus filtering, we analyze the fraction of ‘correct positives’ f _cp, that is, the fraction of cirrus-contaminated pixels that were correctly identified as such, and the fraction of ‘correct negatives’ f _cn, that is, the fraction of cirrus-free pixels that were correctly identified as such. The lower panels of Figure 4 show f _cp (blue line) and f _cn (red line) as a function of R _CO2 and R _H2O, respectively. Obviously, the goal is to have both f _cp and f _cn as large as possible. As expected, with larger R _CO2 and R _H2O, fewer pixels no matter cirrus-contaminated or cirrus-free will be identified as cirrus-contaminated by the filter. Therefore, f _cp decreases with the value of the ratios. On the other hand, f _cn increases with the larger value of the ratios, but more cirrus-free cases (e.g., with large aerosol load) are incorrectly identified as cirrus-contaminated. So, the optimal value for the ratios is a compromise between the two effects. By using a threshold of R _CO2 =1.034, 90% of cirrus-free pixels are preserved and 90% of cirrus-contaminated pixels are filtered out. We consider this a good compromise. A filter based on R _H2O is not able to achieve a 90% success rate, so we consider it inferior compared to the filter based on R _CO2. Therefore, we apply R _CO2 < 1.034 as a cirrus filter to obtain cirrus-free pixels for RemoTAP retrievals in its baseline setup. The cirrus filter will be applied after filtering for thick clouds. In our ensemble, 18.5% of the pixels remain after this first filter step. Then, of the pixels that are free of (thick) clouds, 80% remains after applying the cirrus filter (∼15% of the total).

After applying the cirrus filter, we perform the MAP-CO2I and CO2I-only retrievals for the remaining pixels using RemoTAP in its baseline setup (not fitting cirrus). We apply an a posteriori filter based on goodness-of-fit keeping pixels with χ ² < 1.5 for MAP-CO2I retrievals and χ ² < 1.2 for CO2I-only retrievals.

The upper panel in Figure 5 shows the performance of the MAP-CO2I retrievals. After applying the χ ² < 1.5 filter, about 85% of the processed pixels remain, and 95% of cirrus-contaminated pixels are filtered out. The root mean square error (RMSE) and bias for XCO₂ are 0.6303 and 0.0741 ppm, respectively. These small errors are possible because the aerosol properties are accurately retrieved. The RMSE and bias of AOT at 765 nm are 0.0277 and 0.0032, respectively. Table 6 summarizes the retrieval accuracy of the most important aerosol properties. Overall, the listed aerosol properties are accurately retrieved from the MAP-CO2I combination and are mostly within the aerosol requirements formulated for climate research (Mishchenko et al., 2004; Hasekamp O. P. et al., 2019).

The results of the baseline CO2I-only retrieval are shown in the middle panel of Figure 5. Here, the XCO₂ product has an RMSE of 1.4742 ppm and a bias of 0.1945 ppm. So, the accuracy is significantly degraded compared to the MAP-CO2I retrievals, where the RMSE on XCO₂ is by a factor of 2.3 larger than MAP-CO2I retrievals. The larger errors are caused by a poor capability to retrieve the atmospheric aerosol, with an RMSE and bias on AOT at 765 nm of 0.12 and −0.08, respectively.

In order to illustrate the effect of light path modification by aerosol, we also show the non-scattering retrievals of XCO₂ using measurements of SWIR-1 and SWIR-2 bands in Figures 5E and F, respectively. By applying the filters of χ ² < 1.05 for SWIR-1 and χ ² < 2.0 for the SWIR-2 band, the retrievals can achieve similar passing rates compared to the two cases of full-physics retrievals. Without the light path correction of fitting aerosol, the XCO₂ retrievals at the SWIR-1 band have a large positive bias of 1.513 ppm and an RMSE of 3.6157 ppm, while the retrievals at the SWIR-2 band have a larger negative bias of −3.3184 ppm and an RMSE of 5.4022 ppm. By fitting aerosol in the CO2I-only retrieval case, the accuracy of XCO2 retrievals is improved by a factor of 2 and 4, compared to the non-scattering results for SWIR-1 and SWIR-2, respectively, in terms of RMSE, and the magnitudes of biases are reduced by a factor of > 8. For the combined MAP-CO2I retrievals, the improvements in RMSE are a factor of 5.7 and 8.5, respectively. This demonstrates the large effectiveness of light path correction for the combined MAP-CO2I retrievals.

Figure 6 shows the dependence of XCO₂ errors (ΔXCO₂) as a function of SZA, AOT at 1,600 nm, COT at 1,600 nm, and the BRDF scaling factor at 2,000 nm for both retrieval methods (upper panels for MAP-CO2I retrievals and lower panels for CO2I-only retrievals). The fraction of retrievals with an error larger than 3 ppm (the plotting range) is denoted as f _out. The MAP-CO2I retrievals show virtually no dependence of ΔXCO₂ on SZA and AOT, and only a slight dependence on the surface albedo at 2,000 nm (somewhat larger error at low albedo than at large albedo). Regarding the dependence on surface albedo, it should be noted that aerosols have different effects on the light path for low and high albedo, respectively. For low albedo, aerosols cause a shortening of the light path because the light is reflected back to space at higher altitudes, whereas for high surface albedo, aerosols cause a longer light path because of multiple scattering of light between the surface and the aerosol layer. For the dependence on COT (for the few pixels that were missed by the cirrus filter), we see an underestimation of XCO₂ for all COT values. The CO2I-only retrievals show much more pixels with an error larger than 3 ppm. Also, we see a clear dependence of ΔXCO₂ on AOT and a much stronger dependence of surface albedo than for the MAP-CO2I retrievals. The dependence of ΔXCO₂ on AOT and surface albedo may manifest itself in regional and seasonal dependent biases which hamper accurate quantification of CO₂ emissions on a regional/city scale. For the dependence on COT, we also see an underestimation at large COT but more overestimation at smaller COT than for the MAP-CO2I case.

Figure 7A shows the RMSE, bias, and fraction of successful retrievals (of processed pixels and of all pixels without a priori filtering) for different RemoTAP setups: MAP-CO2I retrievals in baseline setup on cirrus-filtered (cf) data (“MAP_CO2I_cf”), CO2I-only in baseline setup on cirrus-filtered data (“CO2I_cf”), MAP-CO2I retrievals in baseline setup on all data where a cirrus filter is not applied (“MAP_CO2I_all”), MAP-CO2I retrievals where cirrus properties are fitted (Cir) on all data (“MAP_CO2I_Cir_all”), CO2I-only retrievals on all data (“CO2I_all”), and proxy retrievals on all data (“Proxy_all”). We found that the CO2I-only measurements do not contain sufficient information to include additional cirrus-fitting parameters, so these retrievals are not included in the analysis.

We have already discussed the results for the first two retrieval setups. It can be seen that for MAP-CO2I baseline retrieval, applying the a priori cirrus filter and the a posteriori χ ² filter leads to a fraction of successful retrievals of 85% of the processed pixels, which yields 68% of all pixels (MAP_CO2I_cf). For the MAP-CO2I baseline retrievals on all pixels (MAP_CO2I_all), we see that the fraction of successful retrievals increases to 74% of all pixels, and the RMSE increases to 0.82 ppm. It is important to note that here a lot of cirrus-contaminated pixels are filtered out by the a posteriori χ ² filter. If we include cirrus properties in the retrieval state vector (MAP_CO2I_Cir_all), we see the fraction of successful retrievals further increasing to 78% (of all pixels) with an RMSE of 0.80 ppm. So, compared to the baseline setup (MAP_CO2I_cf), the fraction of successful retrieval is increased by 10 percent points (of all pixels), but also the RMSE on XCO₂ is increased by 0.12 ppm. The retrieval accuracy of the cirrus parameters is given in Table 6, where we can see that COT is well retrieved with an RMSE of 0.0161 and a bias of -0.0059. For CO2I-only retrievals, the RMSE increases to 1.82 ppm when applied to all pixels (CO2I_all), and the fraction of successful retrievals is 75% of all pixels.

Next, we consider only the true cirrus-contaminated pixels (cc) (Figure 7B). Here, the MAP-CO2I retrieval in its baseline setup (MAP_CO2I_cc) keeps 12.6% of the cirrus-contaminated pixels (1.6% of all pixels) after χ ² filtering. The RMSE and bias for these specific pixels are large: 2.3 and −1.7 ppm, respectively. When we fit cirrus (MAP_CO2I_Cir_cc), the RMSE and bias for these pixels reduce significantly to 1.27 and −0.22 ppm, respectively, where 65% of the cirrus-contaminated pixels pass the χ ² filter (8% of all pixels). For CO2I-only retrievals on cirrus-contaminated pixels (CO2I_cc), we see that the RMSE and bias are −3.17 and −1.66 ppm, respectively, and 30% of the cirrus-contaminated pixels remain after χ ² filtering (3.7% of all pixels).

We conclude that we can effectively filter out cirrus and get very accurate retrievals when we fit aerosol properties in the MAP-CO2I retrieval. We can also include cirrus properties in the retrieval which will lead to a better coverage but reduced accuracy. In the end, the decision whether to fit cirrus properties or not, that is, the choice between better coverage or better accuracy, should be made based on the evaluation of real retrievals.

For the proxy retrievals, cirrus has little effect on the performance, therefore, neither a cirrus filter nor a χ ² filter is applied. By using the truth for X C H 4 apr , the retrievals on all pixels have an RMSE 0.9 ppm and bias 0.309 ppm. It should be noted that these numbers represent an overly optimistic view of the CO₂ proxy retrievals because it is assumed that the a priori XCH₄ is perfectly known. Still, the MAP-CO2I error performance is slightly better than that of the proxy product which demonstrates the strength of this approach. The proxy retrievals are mostly interesting for quantifying point sources of CO₂, which do not simultaneously retrieve methane (so the a priori XCH₄ error shows no spatial pattern) (Krings et al., 2011). The proxy-retrieved XCO₂ can then be compared against background values so that the a priori XCH₄ error cancels out. Proxy retrievals for point sources are only applicable when either a priori XCO₂ or a priori XCH₄ is well known, which is usually not the case for wildfires and coal mines.

Figure 8 shows the performance for XCH₄ for MAP-CO2I and CO2I-only retrievals with the baseline setup, that is, not retrieving cirrus and applying both cirrus filter and χ ² filter. XCH₄ retrievals by the MAP-CO2I method have an RMSE of 3.8395 ppb and a bias of 0.9222 ppb. The RMSE is a factor of 1.6 smaller than for CO2I-only retrievals. This, on a relative scale, agrees very well with the XCO₂ performance estimates.

The upper panel and lower panel of Figure 9 show, respectively, the RMSE and bias of the MAP-CO2I and CO2I-only retrievals as a function of the month for different regions: China (CN), India (IN), US, Europe (EU), Sahara, and Southern Africa. Also, the monthly mean of AOT is indicated in the upper panel of Figure 9. For China and India, CO2I-only retrievals have a large RMSE for XCO₂ (up to 2.8 ppm), and the RMSE roughly follows the variation in AOT. The RMSE exceeds the CO2M precision requirement of 0.7 ppm by a factor of 4. For other regions, the RMSE is smaller and mostly in the range of 1.4–2.0 ppm but still well above the CO2M requirement. MAP-CO2I retrievals show RMSE smaller than the required 0.7 ppm for most regions and months (but not all), and there is no clear correlation with AOT.

For China and India, the bias of the CO2I-only retrieval exceeds the 0.5 ppm bias requirement for many months with a maximum bias of about 2 ppm in India for April/May. For the other regions, the bias for CO2I-only retrievals is significantly smaller but still exceeds the required 0.5 ppm in a number of cases. There is a large spatial variability in the bias distribution for CO2I-only retrievals, for instance, the bias over China is on average positive, but over the Sahara it is negative. This is mostly because the aerosol compositions and loadings are different in different regions, which leads to different optical properties. The difference in aerosol properties in combination with the difference in surface reflectance leads to inter-region differences. In contrast, the bias for the MAP-CO2I retrievals is below 0.5 ppm for virtually all cases. The regional analysis demonstrates even more convincingly than the global analysis the need for a MAP on the CO2M mission.