Retrieved XCO₂ Accuracy Improvement by Reducing Aerosol-Induced Bias for China’s Future High-Precision Greenhouse Gases Monitoring Satellite Mission

Abstract

China is developing the High-precision Greenhouse gases Monitoring Satellite (HGMS), carrying a high-spectral-resolution lidar (HSRL) for aerosol vertical profiles and imaging grating spectrometers for CO₂ measurements at the same time. By providing simultaneous evaluation of the aerosol scattering effect, HGMS would reduce the bias of the XCO₂ retrievals from the passive sensor. In this work, we propose a method to reduce aerosol-induced bias in XCO₂ retrievals for the future HGMS mission based on the correlation analysis among simulated radiance, XCO₂ bias, and aerosol optical depth (AOD) ratio. We exercise the method with the Orbiting Carbon Observatory-2 (OCO-2) XCO₂ retrievals and AOD ratio inferred from the OCO-2 O₂ A-band aerosol parameters at 755 nm and the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) AOD at 532 nm at several Total Carbon Column Observing Network (TCCON) sites in Europe. The results showed that 80% of measurements from OCO-2 were improved, and data from six TCCON sites show an average of 2.6 ppm reduction in mean bias and a 68% improvement in accuracy. We demonstrate the advantage of fused active–passive observation of the HGMS for more accurate global XCO₂ measurements in the future.

1. Introduction

Based on the requirement of the National Civil Space Infrastructure and Long-term Development Plan (2015–2025), China will launch the High-precision Greenhouse gases Monitoring Satellite (HGMS) in 2023. HGMS was designed as an active–passive fusion satellite to carry a wide-swath greenhouse gases monitor (GGM) and an aerosol–cloud high-spectral-resolution lidar (ACHSRL) on a single platform in a 705 km solar synchronous orbit [1,2]. GGM covers a swath of more than 100 km and has four spectral bands: 0.753–0.768 μm, 1.595–1.625 μm, 2.04–2.08 μm, and 2.275–2.325 μm, which are used for the aerosol loads estimation and the column-averaged CO₂ and CH₄ mole fraction (XCO₂ and XCH₄) detection, respectively. ACHSRL consists of a 532 nm iodine-based HSRL and a 1064 nm Mie lidar. By separating the Mie backscatter return from the Rayleigh backscatter return, ACHSRL could directly provide the extinction profile of aerosols and clouds [3]. Furthermore, the optical properties at 532 nm and 1064 nm would help to better identify the aerosol types [4]. These two payloads make HGMS capable of simultaneous observation for CO₂ and aerosol information [5]. Since the influence of aerosols on XCO₂ retrieval is a recognized major issue [6,7,8,9,10], the simultaneous observation of XCO₂ and aerosols of HGMS will help estimate and correct the aerosol scattering effect on XCO₂ retrievals, which will significantly improve the accuracy of global XCO₂ measurements.

The mechanism and function of GGM’s first three bands are similar to those of the Orbiting Carbon Observatory-2 (OCO-2). The added fourth band at 2.276–2.325 µm is specially used for CH₄ retrievals, the second important greenhouse gas. Currently, the mainstream CO₂ retrieval method for passive satellite remote sensing is the full physical retrieval algorithm [9,11,12,13]. The aerosol information used in the algorithm is mainly from models or reanalysis data. However, aerosols have complex subtypes and relatively rapid temporal and spatial variability in the real atmosphere [14,15]. Aerosol observations from passive remote sensing can only retrieve the column aerosol optical depth (AOD) with restricted accuracy [16,17]. As a result, their ability to correct the aerosol scattering effect on XCO₂ retrievals is very limited. For example, it is easy to misestimate the influence on the upper atmospheric structure due to the lack of the ability to distinguish the upper layer (stratospheric) thin clouds and aerosols essentially. The data shows about a 20% deviation between cloud screen products from OCO-2 and cloud mask products from the Moderate Resolution Imaging Spectrometer (MODIS) [6]. The clouds and aerosols located at higher altitudes, such as stratospheric sulfate from volcanic eruptions, could cause larger XCO₂ bias (100 ppm) even with a small optical depth (such as AOD < 0.3) [7,18].

Studies of XCO₂ retrieval accuracy improvement considering the aerosol-induced bias have been performed. Before the launch of OCO-2, Wunch et al. proposed a method to reduce the XCO₂ systematic bias for the Greenhouse gases Observing SATellite (GOSAT) by constructing the correction equation between the XCO₂ and several parameters (albedo, pressure, and airmass) that may cause XCO₂ bias [10]. Furthermore, Guerlet et al. found the correlations between the XCO₂ bias and aerosol parameters from GOSAT and tried to reduce the interference from the aerosol scattering effect on account of the linear regression fit [7]. Based on Wunch and Guerlet’s research, O’Dell et al. constructed multilinear bias-correction fitting to the surface pressure, CO₂ gradient, and AOD to obtain the correlation with XCO₂, which helped to improve the OCO-2 XCO₂ retrieval and update the OCO-2 version 8 algorithm [18]. All the corrections aforementioned used the Total Carbon Column Observing Network (TCCON) measurements as the truth [19,20,21]. These XCO₂ bias-corrected algorithms almost rely on the predetermined aerosol models, which are not completely consistent with the real atmosphere.

Active remote sensing, such as lidar, can acquire the real-time and accurate vertical structure of the atmosphere [3,22,23]. Using lidar-observed aerosol products to improve XCO₂ accuracy has been studied and proved to be effective. Uchino et al. reduced GOSAT retrieval errors by introducing the vertical distribution of aerosols and cirrus observed by ground-based lidar to the retrieval algorithm [24]. Taylor et al. compared OCO-2 observation with collocated Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP) measurements to analyze the performance of the OCO-2 cloud screening algorithm and found the advantage of CALIOP in identifying high, thin clouds [25]. Merrelli et al. used the CALIOP-derived aerosol products as prior for the OCO-2 physical retrieval [26], which suggested that the OCO-2/CALIOP tight formation flying in A-train is excellent for studying active–passive lidar/OCO-2 retrievals.

HGMS will carry active–passive fusion instruments, providing the real-time aerosol data to correct aerosol-induced bias in XCO₂ retrievals. To demonstrate the capability of HGMS to improve the XCO₂ accuracy, we used the OCO-2 and CALIOP data due to their tight flying formation in the following analysis [26].

In this study, a retrieved XCO₂ accuracy improvement method for HGMS was proposed based on the correlation analysis between XCO₂ bias and AOD ratio. Solar radiances at the weak CO₂ absorption band under different atmospheric scenes were firstly simulated for HGMS to quantify the aerosol scattering effect to the radiance bias. Then, the OCO-2 data and CALIOP AOD were matched to exercise our proposed method. Two-year (2015–2016) OCO-2/CALIOP measurements over land surface and six TCCON sites in Europe were analyzed. With reference to Guerlet et al. and O’Dell et al. [7,18], segmented linear correlation was calculated between XCO₂ bias and AOD ratio at four data sites, which accounts for the impact from different aerosol species, and used to correct XCO₂ at independent validation sites with close background at Bremen and Paris. Due to the complexity of aerosol layer height and subtype distribution, this paper mainly focuses on the seasonal and regional influence of CO₂ measurements caused by aerosol loadings in the troposphere.

After a brief introduction of HGMS, the construction method of HGMS’s solar radiance simulation and the substituted datasets of HGMS used in this work are described in Section 2. Section 3 first presents the simulated radiance distribution and analyzes the correlation between the simulated radiance bias and AOD ratio. Then, the seasonal and regional correlation between XCO₂ bias and AOD ratio are presented. The correlation analysis results are used to improve the OCO-2 official bias-corrected products for six TCCON sites in Europe. Section 4 provides the conclusion of this study.

2. Materials and Methods

2.1. Basic Configuration of HGMS

The detailed parameters of ACHSRL and GGM of HGMS can be found in

Table 1. Parameters of ACHSRL and GGM of HGMS.

ACHSRL		GGM
Parameter	Value	Parameter	Value
Laser wavelength	532/1064 nm	Spectral range	0.753–0.768/1.595–1.625/ 2.04–2.08/2.275–2.325 μm
Pulse energy	150/110 mJ
Frequency	20 Hz	Spectral resolution	<0.04/0.07/0.09/0.1 nm
Laser bandwidth	<100 MHz	Signal-to-noise ratio	>350/340/230/200
Filter bandwidth	<2 GHz	Spatial resolution	<3 km
Divergence angle	0.1 mrad	Swath	>100 km
Field of view	0.2 mrad

. GGM is essentially an imaging grating spectrometer, collecting high-resolution spectra of reflected sunlight in the O₂ A-band, the weak and strong CO₂ absorption band, and the CH₄ absorption band. Similar to OCO-2, the O₂ A-band spectrum has a high sensitivity to changes in the light path [27]. It will be primarily used to determine the surface pressure and screen for optically thick clouds [25]. Furthermore, it can be used to quantify the AOD and vertical distribution of aerosols for soundings that are sufficiently optically thin (AOD < 0.3) to yield accurate estimates of XCO₂. The other active lidar instrument, ACHSRL, consists of a 532 nm HSRL and a 1064 Mie scattering lidar, which is capable of aerosol and cloud optical properties profiling. The real-time aerosol and cloud information provided by ACHSRL could replace the model data during the GGM XCO₂ retrieval and help to accurately estimate the aerosol scattering effect.

The configurations and detection mechanism of GGM and ACHSRL are similar to OCO-2 and CALIOP. All these instruments are loaded on their corresponding satellite in a 705 km solar synchronous orbit (CALIOP lowered its orbit from 705 km to 688 km after September 2018) [18,28]. Their orbit inclination and repeat cycles are all 98.2° and 16 days. ACHSRL and CALIOP are both dual-wavelength lidar with Nd:YAG lasers that emit pulses at 532 and 1064 nm. Their pulse energy at 532 nm is 110 mJ and the repeat frequencies are both close to 20 Hz. Furthermore, their divergence angles are 0.1 mrad and the field of view of ACHSRL is 0.2 mrad, while the field of view of CALIOP is 0.13 mrad [28]. Considering the same orbit altitude, inclination angle, and laser repeat frequency, the global laser footprint coverage of ACHSRL and CALIOP could be assumed to be the same. Besides the instrument configuration, ACHSRL and CALIOP are both polarized lidar collecting the backscattering light to retrieve the optical properties of aerosols and clouds. The spectral ranges of OCO-2 are 0.76–0.77, 1.59–1.62, and 2.04–2.08 μm, respectively, which are nearly the same as GGM’s first three bands. The spatial resolution of OCO-2 is 1.29 × 2.25 km², while GGM’S resolution is close to 3 × 3 km² [18]. Both of their retrieval methods are the full physical retrieval algorithm. Though CALIOP and OCO-2 are loaded on different satellites, they could still acquire the collocated near-simultaneous measurements due to the A-train formation. Since HGMS has not been launched yet, simulated GGM and ACHSRL observations are substituted for OCO-2 and CALIOP datasets.

The product information can be found in

Table 2. Datasets used for the correlation analysis between XCO₂ bias and aerosol scattering effect.

	Instrument	Product Name	Parameters
XCO2	OCO-2	L2 Standard V10r	XCO2
		L2 Lite FP V10r	Bias corrected XCO2
	TCCON	GGG2014	XCO2
Aerosol	OCO-2	L2 Standard V10r	Aerosol parameters at 755 nm
	CALIOP	L2 5km A/CPro V4-20	Extinction coefficient at 532 nm
			CAD score

. The XCO₂ from OCO-2 Standard products were used to analyze the correlation between the bias of XCO₂ and aerosol loads, while the Lite products were used to validate the feasibility of the correction algorithm because of its better agreement with the ground base stations and model data, which is not in the Standard products [29]. The TCCON XCO₂ products are used as the truth [30,31]. The aerosol information analyzed was the aerosol parameters from the O₂ A-band at 755 nm provided by OCO-2 Standard products and the observed extinction coefficient profiles at 532 nm provided by CALIOP; the CALIOP AOD is assumed as the true value. The cloud–aerosol discrimination (CAD) scores of CALIOP products were used to filter data with cloud contamination and low confidence [28].

2.2. Radiance Simulation for HGMS

To better understand the HGMS XCO₂ bias caused by aerosol scattering effect, simulation for the solar radiance at the HGMS’s weak CO₂ absorption band under different aerosol optical properties was calculated. The radiative transfer model used for atmosphere absorption is the Line-By-Line Radiative Transfer Model (LBLRTM v12.9) [32,33]; multiple scattering simulation is calculated by library for radiative transfer (libRadtran v2.0.4) with the discrete ordinates radiative transfer solvers (DISORT) [34,35,36]; solar irradiance spectrum is simulated by the Kurucz compilation [37]. Other parameters chosen in this study mainly refer to Mao’s work, which can be found in

Table 3. Input parameters of radiance simulation for HGMS.

Parameters	Value
CO2 level	360 ppm
Atmosphere profiles	1976 U.S. Standard Atmosphere
Lambertian surface reflectance	0.2
Solar zenith angle	32°
Satellite viewing angle	nadir
Solar irradiance spectrum	Kurucz compilation
Spectral band range	6300–6400 cm−1

[8]. The vertical resolution of the atmosphere profiles (temperature, pressure, etc.) is 100 m and the surface reflectance is 0.2, considering the real situation of the European land surface [38].

Besides the molecular absorption and the surface reflection, the aerosol scattering effect is also considered. Different aerosol scenes, including marine, rural, and urban aerosols, were tested, which were the main aerosol types over European land [39,40,41]. All the aerosols are distributed below 2 km. All aerosol models are taken from LOWTRAN and the aerosol optical properties are calculated by Mie theory and the Henyey–Greenstein phase function according to Mao et al. [8], which can be found in

Table 4. Optical properties of different aerosols used in the radiance simulation.

Aerosol	Single Scattering Albedo	Asymmetry	Visibility/km
Marine	0.98	0.72	23
Rural	0.81	0.64	5
Urban	0.52	0.63	5

. All the models and parameters are input in the libRadtran software with Fortran program [35].

2.3. XCO₂ and AOD Datasets for HGMS

2.3.1. Data Screening and Fusion

To acquire matched comparison data pairs among OCO-2, CALIOP, and TCCON and eliminate the errors caused by temporal and spatial mismatch, a few steps were taken to screen the observations from these instruments. The dates when there are no data available for satellite footprints within a 200 km radius from the TCCON site, and the minimum distance between OCO-2 and CALIOP footprint is over 10 km, are removed. The number of available dates for different sites is counted in

Table A1. Summary of the available dates of different sites used for correlation analysis.

Site (Location)	Date
Bialystok (53.23° N, 23.025° E)	26 Feb 2015 23 Mar 2015 25 Mar 2015 30 Mar 2015 12 May 2015 11 Jun 2015 13 Jun 2015 6 Jul 2015	20 Jul 2015 31 Jul 2015 1 Sep 2015 22 Sep 2015 10 Oct 2015 25 Nov 2015 18 Mar 2016 25 Mar 2016	19 Apr 2016 26 Apr 2016 3 May 2016 21 May 2016 28 May 2016 4 Jun 2016 22 Jun 2016 29 Jun 2016	24 Jul 2016 31 Jul 2016 7 Aug 2016 1 Sep 2016 8 Sep 2016 19 Sep 2016 22 Nov 2016 29 Nov 2016
Karlsruhe (49.10° N, 8.44° E)	7 May 2015 17 Jun 2015 13 Aug 2015 21 Sep 2015 24 Nov 2015 19 Dec 2015	26 Dec 2015 20 Jan 2016 2 May 2016 10 Jun 2016 5 Jul 2016 7 Jul 2016	21 Jul 2016 30 Jul 2016 1 Aug 2016 6 Aug 2016 8 Aug 2016 31 Aug 2016	7 Sep 2016 14 Sep 2016 16 Oct 2016
Bremen (53.10° N, 8.85° E)	13 Aug 2015 7 Sep 2015	11 Apr 2016 2 May 2016	24 Aug 2016 31 Aug 2016	28 Nov 2016
Garmisch (47.48° N, 11.06° E)	18 Feb 2015 18 May 2015 10 Jun 2015 17 Jun 2015 5 Jul 2015	6 Aug 2015 13 Aug 2015 21 Sep 2015 10 Nov 2015 24 Nov 2015	20 May 2016 28 Jun 2016 5 Jul 2016 31 Aug 2016 7 Sep 2016	27 Oct 2016 3 Nov 2016
Orleans (47.97° N, 2.113° E)	21 Feb 2015 23 Mar 2015 25 Mar 2015 1 Apr 2015 10 May 2015 19 May 2015 21 May 2015 13 Jun 2015 15 Jul 2015	22 Jul 2015 29 Jul 2015 10 Sep 2015 24 Sep 2015 12 Oct 2015 19 Oct 2015 26 Oct 2015 4 Dec 2015 22 Dec 2015	20 Mar 2016 27 Mar 2016 10 Apr 2016 28 Apr 2016 5 May 2016 6 Jun 2016 24 Jun 2016 8 Jul 2016 15 Jul 2016	9 Aug 2016 27 Aug 2016 3 Sep 2016 10 Sep 2016 28 Sep 2016 19 Oct 2016 30 Oct 2016 6 Nov 2016
Paris (48.846° N, 2.356° E)	12 May 2015 21 May 2015 15 Jul 2015	10 Sep 2015 24 Sep 2015 26 Oct 2015	24 Jun 2016 10 Sep 2016 28 Sep 2016

. CALIOP absolute CAD score higher than 70 indicates high confidence from its official data introduction (positive values signify clouds, whereas negative values signify aerosols) [4]. Since the cloud existence would lead to XCO₂ retrieval failure, CALIOP profiles with all the CAD scores lower than −70 are retained for further analysis. After orbit registration, cloud screening, and quality control, the remaining OCO-2 soundings and CALIOP footprints can be placed into pairs for the correlation comparison. The procedure is to choose the CALIOP footprints as the center point and draw a rectangle around the point along the orbit direction as the fusion range. The rectangle size is 5 km (the horizontal resolution of CALIOP L2 products) along the orbit direction and 10 km across the orbit direction, which is displayed in the dotted red circle in Figure 1 [42]. The soundings of OCO-2 located in the fusion range are averaged. The averaged AOD of OCO-2 and CALIOP are denoted as τ₇₅₅ and τ₅₃₂, respectively. The averaged XCO₂ of OCO-2 is denoted as X_{CO2, OCO2}. These fused data are used in Section 3.2 for data correlation analysis.

Six TCCON sites in Europe, including Bialystok, Bremen, Garmisch, Karlsruhe, Orleans, and Paris, are analyzed in this work [10,21]. Bremen and Paris are used to validate the XCO₂ accuracy improvement method, while the other four sites are used to analyze the correlation. Due to the high cloud fraction, large surface snow cover, and limited direct sunlight in winter (December–January–February, DJF) and spring (March–April–May, MAM), most of the TCCON and OCO-2 measurements are in summer (June–July–August, JJA) and fall (September–October–November, SON). Liang et al. compared the nadir mode OCO-2 retrievals and TCCON measurements by restricting the period within an hour and the area within ±2° latitude and ±2.5° longitude of the TCCON sites. The comparison found that the effect of averaging kernels changes could be neglected with the differences less than ±0.021 ppm on average [43]. In this work, a 200 km radius area from the TCCON sites and ±1 h period measurements were also chosen to obtain sufficient data and guarantee spacetime correlation (as shown in Figure 1). Synergistic data in 2015 and 2016 were analyzed because CALIOP has experienced lower energy laser shots within and near the South Atlantic Anomaly region since 2017 and it also lowered its orbit from 705 km to 688 km to resume formation with CloudSat in 2018. Furthermore, the TCCON data availability in Europe is higher in 2015 and 2016 [44]. Totally, 2625 comparison pairs were obtained, covering 111 days over two years.

2.3.2. Data Segmentation

For different TCCON sites, different seasons, and different true values of AOD are considered. The data are divided into four seasons (DJF, MAM, JJA, SON), and the data in the same season are divided into four regions based on the τ₅₃₂. When τ₅₃₂ = 0, it means the atmosphere should be clean without any aerosol. When τ₅₃₂ ≠ 0, it means there are aerosols in the atmosphere, which may lead to bias of XCO₂ measurements. Aben et al. found typical errors in the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) CO₂ column for nadir observations over land surfaces with AOD range from 0.1 to 0.3 [45]. Crisp et al. identified 0.3 as a reasonable AOD threshold below which to attempt XCO₂ retrievals [46]. Therefore, we divided the AOD into several regions based on the distribution of the values, which are 0 < τ₅₃₂ ≤ 0.1, 0.1 < τ₅₃₂ ≤ 0.3, τ₅₃₂ > 0.3. Here, we use XXX _{(aod1, aod2)} to describe a case in the correlation at a TCCON site. For instance, MAM _{(0.1, 0.3]} means the correlation is analyzed when 0.1 < τ₅₃₂ ≤ 0.3 and in spring.

The Ångström exponent law was introduced here to estimate the AOD ratio between different wavelengths:

τ/τ₀ = (λ/λ₀)^−k.

(1)

λ is the wavelength and τ is the corresponding AOD to be estimated; λ₀ is the reference wavelength and τ₀ is the corresponding reference AOD; k is the well-known Ångström exponent [47]. For example, Amiridis et al. used AErosol RObotic NETwork (AERONET) and CALIOP measurements to establish a multi-wavelength aerosol optical databases at 355, 532, 1064, 1570, and 2050 nm by employing the Ångström exponent law [47]. Holben et al. also used AERONET measurements to calculate the Ångström exponent at 440/870 nm [48]. Based on their calculation and statistics, the approximation of k_755/532 with the same size distribution and aerosol types could be roughly estimated as 1.6, which means the τ_755/532 could be estimated as 0.57. If τ_755/532 is far beyond or lower than 0.57, the τ₇₅₅ may have a large bias compared with its true value, assuming the τ₅₃₂ is accurate. The XCO₂ bias ∆X_CO2 can be expressed as [19]

∆X_CO2 = (X_{CO2, TCCON} − X_{CO2, OCO2})/X_{CO2, TCCON} × 100%,

(2)

where X_{CO2, TCCON} is the TCCON XCO₂ measurements averaged within ±1 h when OCO-2 passed through the site (assumed the amount of OCO-2 measurements is N and TCCON measurements is M). Meanwhile, the uncertainty of XCO₂ bias δ_XCO2 is calculated as [49]

(δ_XCO2/∆X_CO2)² = (δ_{XCO2, TCCON}/X_{CO2, TCCON})² + (δ_{XCO2, OCO2}/X_{CO2, OCO2})².

(3)

δ_{XCO2, TCCON} and δ_{XCO2, OCO2} are the uncertainty of TCCON and OCO-2 XCO₂ measurements, computed as [19]

$${({\delta }_{XCO}{{}_2}_{, OCO}{}_2)}^2=\frac{1}{N-1}{\displaystyle \sum _{i=1}^N\left(X_{CO2, OCO2}^i-{\overline{X}}_{CO2, OCO2}\right)},$$

(4)

$${\left({\delta }_{XCO}{}_{2, TCCON}\right)}^2=\frac{1}{M-1}{\displaystyle \sum _{i=1}^M\left(X_{CO2, TCCON}^i-{\overline{X}}_{CO2, TCCON}\right)}.$$

(5)

Further analysis between τ_755/532 and ∆X_CO2 is described in Section 3.2.

3. Results and Discussion

3.1. Radiance Simulation Analysis

The radiance distributions of top of atmosphere in HGMS’s weak CO₂ absorption band under different aerosol types and optical depths are illustrated in Figure 2. In Figure 2a, the radiance changes under different aerosol scenes with the same AOD. The radiation deviations caused by aerosols can be positive and negative compared with the clear sky, and the largest and smallest deviations occur under urban and rural aerosol scenes. Figure 2b shows the radiance distribution influenced by different AOD under urban aerosol scene. It is easy to see that the radiation deviations increased when the AOD increased.

Here, we defined R as the ratio between the online and offline radiance (r_on and r_off), R = r_on/r_off. The offline was set to 6382.7 cm⁻¹ and the online was from 6300 cm⁻¹ to 6340 cm⁻¹. The ratio change means the ratio variation when the XCO₂ increased by 1% under PBL, which is △R = (R’ − R)/R × 100%. R’ means the ratio when XCO₂ is 363.6 ppm. The ratio change distribution under different aerosol scenes can be found in Figure 3a–d. The AODs set in Figure 3a–c are all 0.5. The distribution indicates that the aerosol existence leads to positive ratio change, which means that the aerosol scattering effect could lengthen the optical path of CO₂ detection, and it resulted in an overestimated XCO₂ retrieval. Additionally, the largest △R was introduced under the marine aerosol scene while the smallest one occurred under urban aerosol scene with the same AOD of 0.5. Figure 3e–h shows the radiance change when the XCO₂ increased by 1% with different AOD under urban aerosol scene. In the clear sky (Figure 3h), the peak value (−0.4%) of the spectrum is the largest compared with the aerosol scenes, which means that the detection sensitivity would be the largest and the relative error would be the smallest where there is no aerosol load. The sensitivity is only half of that in the clear sky when the AOD is 0.5 in Figure 3g, which would demand higher signal-to-noise ratio (SNR). The bias in Figure 3e–h means the peak value bias of the spectrum compared with the clear sky.

The relationship between radiance bias and AOD ratio (τ_755/532) was also simulated to analyze the XCO₂ bias influenced by the inaccurate OCO-2 aerosol model products. The results can be found in Figure 4. As described in Section 2, the AOD at 532 nm was assumed to be the true value. The radiance bias ∆r_off in Figure 4a was calculated at the offline wavelength with different AOD input, ∆r_off = r_off,755 − r_off,532. The intercept on the x-axis is 0.57, which is assumed to be the accurate value of τ_755/532. When τ_755/532 = 0.57, it could be recognized that there is no AOD bias at 755 nm and the radiance bias is correspondingly zero. Figure 4a shows the radiance bias influenced by the τ_755/532 with different AOD values (0.05–0.4) at 532 nm. The higher τ_755/532 could lead to higher radiance bias. Figure 4b shows the relationship between median radiance changes ∆r_m and τ_755/532. ∆r_m are the medians of the spectrum radiance difference when the AOD inputs are set to 755 nm and 532 nm, respectively.

∆r_m = median|(r_band,755 − r_band,532)/r_band,532 × 100%|.

(6)

r_band means the radiance spectrum distributed from 6300 to 6400 cm⁻¹. The symbol of median || means calculating the median of radiance changes. The ∆r_m in Figure 4b are all positive, which means the XCO₂ retrieval would be overestimated. The larger the τ_755/532, the more the results are overestimated. Additionally, when the AOD at 532 nm is lower than 0.1, the radiance bias and changes have a linear correlation with the τ_755/532. With the increase of AOD, the correlation gradually turns into a nonlinear form, which is caused by complicated multiple scattering effect with higher aerosol loading. Consequently, in the high-AOD scenes, the quality of HGMS’s CO₂ retrieval is worrying. With limited dynamic detection range and spectral resolution, invalid retrieval would occur for HGMS if there is large bias for the AOD estimation.

3.2. Correlation Analysis

The simulation results in Figure 4 indicate that the correlation between AOD ratio and HGMS’s radiance bias could be linear or near-linear. As for the correlation between HGMS’s XCO₂ bias and AOD ratio, the correlation coefficients with the division of different cases are presented in

Table 5. The correlation coefficients between HGMS’s XCO₂ bias and AOD ratio at different seasons, different τ₅₃₂, and different TCCON Sites. LD means lack of data. LC means low correlation.

Site	Season	τ532 = 0	(0, 0.1]	(0.1, 0.3]	(0.3, ∞)	NHC/NDA	Data Utilization
Bialystok	DJF	LD	LD	0.62	LD	17/18	94.4%
	MAM	LC	LC	0.56	0.58	95/199	47.7%
	JJA	0.56	0.62	LC	0.64	221/285	77.5%
	SON	LC	0.61	LC	LC	122/275	44.4%
Garmisch	DJF	0.63	LD	LD	LD	50/50	100%
	MAM	LD	0.57	LC	LD	29/39	74.5%
	JJA	0.58	LC	0.56	0.69	192/267	71.9%
	SON	LC	LC	LC	LC	0/281	0
Karlsruhe	DJF	0.58	LC	LD	LD	42/94	44.7%
	MAM	LC	0.64	LC	LD	45/93	48.3%
	JJA	LC	0.60	LC	LC	70/187	37.4%
	SON	LC	LC	0.57	LD	106/253	41.9%
Orleans	DJF	LD	LD	0.79	LD	10/17	58.8%
	MAM	LC	LD	0.64	LC	35/129	27.1%
	JJA	0.62	LC	LC	0.70	120/196	61.2%
	SON	0.59	LC	0.64	LD	143/242	38.5%
TOTAL				1297/2625			49.4%

, where LD means lack of data and LC means low correlation in terms of a linear relationship (correlation coefficient R < 0.5). N_DA means the amount of data available after the screening method described in Section 2. N_HC means the amount of data of high correlation (R > 0.5).

The results indicate that about half of the cases (22/48, regardless of the LD cases) have a certain correlation between ∆X_CO2 and τ_755/532 (τ₇₅₅). For the 22 cases with a certain correlation, we assume the linear regression equation is expressed as y × 100 = b + a × x, where x means τ_755/532 (τ₇₅₅) and y means ∆X_CO2. The results are displayed in

Table 6. Lookup table of linear regression equation coefficients of 22 cases.

Site	Case	b	a	Site	Case	b	a
Bialystok	JJAτ532 = 0	0.09	4.71	Garmisch	DJFτ532 = 0	0.37	11.8
	JJA(0, 0.1]	0.49	0.53		JJAτ532 = 0	0.18	5.44
	SON(0, 0.1]	0.64	0.47		MAM(0, 0.1]	0.057	0.39
	DJF(0.1, 0.3]	0.76	0.27		JJA(0.1, 0.3]	0.14	0.34
	MAM(0.1, 0.3]	−0.20	0.59		JJA(0.3, ∞)	0.062	1.58
	MAM(0.3, ∞)	0.62	0.90	Orleans	JJAτ532 = 0	0.18	5.44
	JJA(0.3, ∞)	0.28	1.89		SONτ532 = 0	0.42	7.91
Karlsruhe	DJFτ532 = 0	0.15	9.98		DJF(0.1, 0.3]	0.41	1.27
	MAM(0, 0.1]	0.15	0.48		MAM(0.1, 0.3]	0.27	0.35
	JJA(0, 0.1]	0.18	0.26		SON(0.1, 0.3]	0.13	0.52
	SON(0.1, 0.3]	0.48	0.51		JJA(0.3, ∞)	0.30	1.18

.

Correlation plots of the six cases are displayed in Figure 5 and Figure 6, with τ_{532 = 0} and τ_{532 ≠ 0}, respectively.

3.3. Validation

The difference in surface albedo may lead to systematic errors in CO₂ retrievals [18]. To make the correlation analysis as independent of surface albedo as possible, the weekly averaged surface albedo at the six TCCON sites from April to December (no data available for validation in January, February, or March at Bremen and Paris) were calculated based on the MODIS MCD43C3 products at band 6 (1.628–1.652 μm) [50]. The results can be seen in Figure A1. The weekly averaged albedo at Bialystok is 0.42 ± 0.01, while the albedo at Bremen is 0.44 ± 0.01. The weekly averaged albedo at Orleans is 0.25 ± 0.02, and the albedo is 0.23 ± 0.01 at Paris. We decide to use the correlation analysis results at Bialystok to improve the Bremen observation, and we used the Orleans results to improve the Paris observation due to their similar weekly albedo change.

According to Equation (2) and the linear regression between ∆X_CO2 and τ_755/532, it can be derived that

∆X_CO2 = b + a × τ_755/532.

(7)

The optimization OCO-2 observations are denoted as X_{CO2, NEW}, which is used to replace the X_{CO2, TCCON} in Equation (2) because it is assumed that the TCCON measurements are the true values, and then we can obtain the optimization X_{CO2, NEW} after approximate linear regression correction by calculating

(X_{CO2, NEW} − X_{CO2, OCO2})/X_{CO2, NEW} × 100% = b + a × τ_755/532.

(8)

The X_{CO2, NEW} can be expressed as

X_{CO2, NEW} = X_{CO2, OCO2}/(1 − b − a×τ_755/532) × 100%.

(9)

The optimization X_{CO2, NEW} is obtained based on the linear regression coefficient a and b selected from Table 6. The OCO-2 Lite products are used to validate the feasibility of the correction algorithm. Figure 7 shows the comparison results between OCO-2 and TCCON measurements for Bremen and Paris sites.

The OCO-2 measurements include official Standard, Lite products, and new results after linear regression correction, noted as X_{CO2, Std}, X_{CO2, Lite}, and X_{CO2, NEW}. The ∆XCO₂ in Figure 7 equals X_{CO2, i} minus X_{CO2, TCCON} (i = Std, Lite, and New), while the legends in Figure 7 represent the statistical results of the absolute value of the XCO₂ difference. The ∆XCO₂ in Figure 7 are divided into several regions, and the heights of different colored histograms represent the counts distributed in the corresponding region.

In Figure 7a–c, the difference between X_{CO2, OCO2} and X_{CO2, TCCON} is decreased from 1.68 ± 0.91 ppm to 0.99 ± 1.03 ppm on September 7 2015, with 0 < τ₅₃₂ ≤ 0.1, while the difference between X_{CO2, Lite} and X_{CO2, TCCON} is 0.77 ± 0.7 ppm. In addition, the difference between X_{CO2, OCO2} and X_{CO2, TCCON} is decreased from 2.59 ± 1.36 ppm to 0.95 ± 1.12 ppm on November 28 2016, with 0.1 < τ₅₃₂ ≤ 0.3, while the difference between X_{CO2, Lite} and is 1.12 ± 1.16 ppm. In Figure 7g–i, the difference between X_{CO2, OCO2} and X_{CO2, TCCON} is decreased from 7.98 ± 3.23 ppm to 4.99 ± 3.63 ppm on May 21 2015, with 0.1 < τ₅₃₂ ≤ 0.3, while the difference between X_{CO2, Lite} and X_{CO2, TCCON} is 4.78 ± 2.75 ppm. The difference between X_{CO2, OCO2} and X_{CO2, TCCON} is decreased from 3.81 ± 0.81 ppm to 1.64 ± 1.58 ppm on September 10 2015, with τ₅₃₂ = 0, while the difference between X_{CO2, Lite} and X_{CO2, TCCON} is 2.26 ± 0.77 ppm.

A total of 16 days (7 days at Bremen and 9 days at Paris) of linear regression correction results are presented in Figure 8. Here, the mean bias (MB) of the XCO₂ is introduced to further indicate the overall bias between the OCO-2 and TCCON measurements. MB can be expressed as

$$\mathrm{MB}=\frac{1}{N}{\displaystyle \sum _N}|X_{CO2, \mathrm{TCCON}}-X_{CO2, OCO2}|.$$

(10)

Both X_{CO2, OCO2} and X_{CO2, TCCON} are daily averaged results. The correlation coefficients R are improved from 0.96 to 0.98 for Bremen and from 0.65 to 0.68 for Paris, while the MB is decreased from 2.39 ppm to 0.89 ppm for Bremen and from 5.86 ppm to 1.07 ppm for Paris.

Table 7. Summary of the XCO₂ among TCCON, OCO-2, and optimization results for Bremen and Paris during 2015–2016. SD_CO2 is the standard deviation of XCO₂.

Site	Season	TCCON	OCO-2 Data				Optimization OCO-2 Data
		XCO2	Count	XCO2	SDCO2	R	Count	XCO2	SDCO2	R
Bremen	MAM	403.82	17	400.46	1.67	0.96	15	401.99	1.22	0.98
	JJA	398.53	155	397.12	1.90		106	398.35	2.09
	SON	399.30	152	396.90	3.86		140	398.74	3.92
Paris	MAM	403.07	58	396.30	4.10	0.65	23	399.96	3.17	0.68
	JJA	399.51	53	394.69	4.64		43	400.52	4.28
	SON	398.67	97	392.69	2.36		97	397.56	2.02

gives more detailed information about the comparisons among TCCON, OCO-2, and optimization results of XCO₂ for Bremen and Paris during 2015–2016.

Linear regression correction results of all six TCCON sites are analyzed and presented in Figure 9. The correlation coefficient R is improved from 0.76 to 0.81, and the MB is decreased from 3.81 ppm to 1.21 ppm.

3.4. Discussion

In Table 5, the overall data utilization is 49.4%, which means that about half of the data is correlated based on the division. The statistical histogram is presented in Figure 10. The rectangular bars with different line types in the same season mean different TCCON sites. The single bar is the data amount for one site in a specific season. The amount of the green pattern is equal to the N_HC in Table 5, while the orange pattern represents data with a lower correlation (N_DA-N_HC). For different seasons, it is clear that the total amounts (N_DA) in DJF and MAM are lower than the other seasons, which contribute only 24.3% of the N_DA due to the high cloud fraction and lack of reflected solar radiation from the sun in winter and spring in the Northern Hemisphere. As for the cloud fraction over land surface in Europe in 2016, Chen et al. found that the highest clear-sky fraction occurs in summer (55.1%), while the lowest occurs in winter (29.5%). The spring and fall have close values, which are 39.3% and 40.0%, respectively [6]. In winter, the Bialystok site has the lowest N_DA because of its lower solar radiance located at high latitude. Additionally, N_DA would be different since the TCCON data availability was different, considering the instrument failure and weather conditions. Based on the official TCCON datasets, there are 394 days of available data for Bialystok, 278 days for Garmisch, 267 days for Karlsruhe, and 415 days for Orleans in 2015–2016, respectively [51].

As for the cases with lower correlation in Table 6, the possible reasons are the following: a) there is limited N_DA in this case, such as the MAM_{(0, 0.1]} in Bialystok, the MAM_{(0.1, 0.3]} in Garmisch, and the MAM_{(0.3, ∞)} in Orleans, the N_DA of which are, respectively, 12, 10, and 12. Lower N_DA may lead to lower correlation because a single datum can significantly influence the correlation; b) though the N_DA is enough, the data are from a single date, which made the results reflect just one day’s situation. For instance, the data in MAM_{(0.1, 0.3]} in Karlsruhe came from only 7 May 2015 and 2 May 2016. The N_DA is 44 on May 2 while the N_DA is 4 on May 7, which represents the case analysis for the single day’s correlation and cannot reflect the seasonal characterization; c) the same aerosol load in the same season can lead to different errors to △X_CO2 if the aerosol subtype and vertical distribution are quite different, which is consistent with the finding in Section 3.1.

Figure 5 indicates that most of the τ₇₅₅ is less than 0.4, which is a reasonable range of the AOD from the OCO-2 aerosol model. The bigger the τ₇₅₅ is, the bigger the △X_CO2 is for most data, revealing that the increase of τ₇₅₅ can lead to the increase of △X_CO2 when τ₅₃₂ = 0 in summer at these two sites. In Figure 6, most of the τ_755/532 is larger than 0.57, which shows that the τ₇₅₅ is overestimated, assuming the τ₅₃₂ is the true value. The correlation is consistent with the simulation results in Figure 4, which means that the linear regression fitting equation applied here is reasonable. Similar to the trends in Figure 5, it can also indicate that the increase of τ_755/532 can lead to the increase of △X_CO2 when τ₅₃₂ ≠ 0 in specific seasons and at specific sites in Figure 6.

As for the validation section, the distribution of stacked histograms in Figure 7 shows significant trends where most Standard products are smaller than the assumed true value of TCCON measurements. In comparison, the Lite products and the optimization results are closer. Though there are some cases where the optimization results are not as improved as the Lite products, the deviations are usually less than 0.2 ppm, which is far less than the uncertainty of each product. It is worth mentioning that the uncertainty of new optimization results is occasionally larger than the Standard products. This is mainly related to the linear regression equation. Not all the residuals are small enough in the fitting regions, which means that if the τ_755/532 near the Bremen site is close to zero or greater than 3, the optimization results may not be satisfied. The 16 days’ validation results in Figure 8 and Table 7 show that 79.7% (424/532) data pairs are successfully improved using this method. For the Bremen site, the differences between the average XCO₂ from TCCON sites and OCO-2 measurements passing nearby regions is decreased from 2.39 ± 2.48 ppm to 0.89 ± 2.41 ppm, with a 1.5 ppm mean bias decrease. For the Paris site, the result is decreased from 5.86 ± 3.70 ppm to 1.07 ± 3.16 ppm, with a 4.79 ppm mean bias decrease. The accuracy is improved by 51.9% and 81.7%, respectively. Though the SD_CO2 of optimization data is not decreased in JJA and SON at the Bremen site from Table 7, it is still acceptable because the increments are less than 0.2 ppm (0.19 and 0.06 ppm, respectively). As for the total six TCCON site measurements in Figure 9, the differences between the averaged XCO₂ from TCCON sites are decreased from 3.81 ± 2.62 ppm to 1.21 ± 1.89 ppm, with a 2.6 ppm mean bias decrease and 68.2% accuracy improvement. In general, the statistical results verify the feasibility of the proposed method to reduce aerosol-induced bias in XCO₂ retrievals for the future HGMS.

4. Conclusions

In this study, an XCO₂ accuracy improvement method was proposed based on the radiative transfer simulation and two-year measurements from OCO-2, CALIOP, and six TCCON sites to demonstrate the advantage of active–passive fusion observation. The method was designed for HGMS XCO₂ observations in the future.

Simulations for the solar radiance at the HGMS’s weak CO₂ absorption band under different aerosol parameters (AOD, subtype) were calculated. For AOD, the radiance change was −0.4% when XCO₂ increased by 1% in the clear sky, but the change would be half of that when the AOD is 0.5; for aerosol subtypes, the largest radiance ratio change was introduced under the marine aerosol scene while the smallest one occurred under the urban aerosol. The lower radiance change represents lower detection sensitivity, which means that higher SNR was needed for the detector and a larger XCO₂ bias would occur in the urban scenes with large AOD. In addition, a near-linear correlation was found between the radiance bias and AOD ratio (τ_755/532), which provided the main idea of our XCO₂ accuracy improvement method.

Under the guidance of radiance simulation, the correlations between XCO₂ bias and AOD ratio of HGMS were analyzed. OCO-2 and CALIOP measurements were substituted for the corresponding datasets of HGMS. First, a screening and fusion method was used to acquire proper data for correlation analysis. Under specific TCCON site, season, and AOD, the XCO₂ bias and AOD ratio have a linear correlation. Among 2625 comparison pairs, 49.4% of them show a certain linear correlation (R > 0.56). For different seasons, the minimum N_DA and N_HC occur in winter due to the high cloud fraction and lack of solar radiation in the Northern Hemisphere. A lookup table of the linear regression equation coefficient was summarized for 22 different cases at four TCCON sites. An XCO₂ accuracy improvement method was proposed based on the lookup table and linear regression equations. The method was then applied to the other two TCCON sites (Bremen and Paris). Under the comparison of OCO-2 official Lite products, 80.0% (424/532) of the data used to validate the linear regression correction method are successfully improved. Similar to the four experimental TCCON sites, there is no retrieved XCO₂ that is bias-corrected in winter because of the cloud existence and solar radiance issues for the two validation sites. For the rest of data in spring, summer, and fall, the XCO₂ bias is decreased from 2.39 ± 2.48 ppm to 0.89 ± 2.41 ppm for Bremen and from 5.86 ± 3.70 ppm to 1.07 ± 3.16 ppm for Paris. All six TCCON data show a 2.6 ppm mean bias decrease and 68.2% accuracy improvement. The correlation coefficient was also improved from 0.76 to 0.81. The feasibility of the XCO₂ accuracy improvement method by reducing aerosol-induced bias was demonstrated.

For future work, more measured data with different aerosol layer heights and subtypes could be considered, as aerosol profiling and classification are the advantage of lidar. Since the XCO₂ data used in this study are from 2015–2016 over land surface in Europe, more data in different years and near different TCCON sites around the world could be analyzed for further verification, which will be useful to update the optimization method applied to HGMS in the future and will help to obtain more accurate XCO₂ observations globally.