remote Obtaining Gradients of XCO2 in Atmosphere Using the Constrained Linear Least-Squares Technique and Multi-Wavelength IPDA LiDAR

: Integrated-path di ﬀ erential absorption (IPDA) LiDAR is a promising means of measuring the global distributions of the column weighted xCO 2 (dry-air mixing ratio of CO 2 ) with adequate accuracy and precision. Most IPDA LiDARs are incapable of discerning the vertical information of CO 2 di ﬀ usion, which is of great signiﬁcance for studies on the carbon cycle and climate change. Hence, we developed an inversion method using the constrained linear least-squares technique for a pulsed direct-detection multi-wavelength IPDA LiDAR to obtain sliced xCO 2 . In the proposed inversion method, the atmosphere is sliced into three di ﬀ erent layers, and the xCO 2 of those layers is then retrieved using the constrained linear least-squares technique. Assuming complete knowledge of the water vapor content, the accuracy of the retrieved sliced xCO 2 could be as high as 99.85% when the signal-to-noise ratio of central wavelength retrievals is higher than 25 (with a log scale). Further experiments demonstrated that di ﬀ erent carbon characteristics can be identiﬁed by the sign of the carbon gradient of the retrieved xCO 2 between the ABL (atmospheric boundary layer) and FT (free troposphere). These results highlight the potential applications of multiple wavelength IPDA LiDAR.


Introduction
At present, nearly half of the total anthropogenic CO 2 emission has been offset by unidentified CO 2 sinks [1]. Dense observations of column xCO 2 (dry-air mixing ratio of CO 2 ), supplementing the existing ground-based CO 2 monitoring networks, are urgently needed to provide stronger constraints on inversions of carbon fluxes using the atmospheric inversion technique [2]. Since 2000, several satellite-based sensors, such as AIRS, SCIAMACHY, and IASI, have been sent into orbit and are able to obtain data of XCO 2 (column-average xCO 2 ). In recent decades, some dedicated greenhouse gas monitoring missions, including GOSAT, OCO-2, and TanSat, have successfully obtained XCO 2 products with higher accuracy and precision than before. It is widely witnessed that these data are useful in deepening our understanding of the carbon cycle process, despite their susceptibility to aerosols and dependency on sunlight as well as insufficient coverage over high-latitude zones [3,4].
Supplementing the passive remote sensing of greenhouse gases, differential absorption LiDAR (DIAL) plays an increasingly critical role in atmospheric CO 2 sensing, owing to the remarkable developments made in both lasers and detectors [5,6]. Integrated path differential absorption (IPDA) LiDAR, a special type of DIAL, relies on backscattered signals by hard-targets, thus providing a higher signal-to-noise ratio (SNR) [7]. Several groups from different countries are dedicated to developing equipment for the active remote sensing of CO 2 equipment based on the integrated-path differential absorption (IPDA) technique, which is expected to obtain CO 2 measurements from space with high accuracy and better coverage at night and at high latitudes [8][9][10]. In recent years, airborne CO 2 -IPDA LiDARs have been extensively tested by different research groups funded by the National Aeronautics and Space Administration (NASA) [5,[11][12][13]. In these experiments, results with errors lower than 1% were obtained and researchers are anticipating the realization of spaceborne CO 2 -IPDA LiDAR in the near future.
To understand the transportation process of CO 2 in the atmosphere, it is important to identify vertical gradients of xCO 2 , especially in the lower troposphere [14,15]. Gradients of xCO 2 between the ABL (atmospheric boundary layer) and FT (free troposphere) can be used to estimate the surface net ecosystem exchange [16]. Vertical gradients of CO 2 would also help us to gain insight into the feedback between the carbon cycle and climate change [15,[17][18][19]. The ground-based dual-wavelength DIAL can measure the range-resolved xCO 2 within 3 km [18,20,21]. However, it relies on aerosol backscattered signals with a relatively lower SNR and thus is unlikely to provide products with accuracies of up to 1 ppm. Moreover, ground-based equipment is difficult to deploy and operate in remote places. The spaceborne CO 2 -IPDA LiDAR is a suitable and effective means to obtain measurements of CO 2 with a high accuracy globally. However, most dual-wavelength IPDA systems can only retrieve column-weighted xCO 2 and are incapable of measuring gradients of xCO 2 .
To tackle the above-stated problem, we propose a novel multi-wavelength inversion method, which is based on the constrained linear least-squares technique. In the field of IPDA/DIAL, the multiple-wavelength strategy was proposed soon after the dual-wavelength differential absorption inversion algorithm. In 1985, R.E. Warren first presented a generalized DIAL signal model and then proposed an inversion method of the concentration of a target material by using maximum likelihood estimation [22]. The main goal of R.E. Warren's method is to retrieve multiple-target materials by means of multiple-wavelength absorption LiDAR. In 1999, Fukuchi presented a multiple wavelength DIAL for measuring range-resolved SO2 and demonstrated that the accuracy can be improved by using averaged signals of a group of on-line/off-line wavelengths instead of a single pair [23]. Later, in 2002, their group successfully developed a five-wavelength SO2-DIAL and obtained range-resolved SO2 concentrations by means of a curve-fitting technique [24]. In 2014, Abshire et al. demonstrated airborne measurements of column-weighted xCO 2 using a 30-wavelength CO 2 -IPDA LiDAR, which is the first multiple-wavelength CO 2 -IPDA LiDAR [5]. In their inversion method, a layered atmospheric model was used to compute an absorption spectrum for a given path, CO 2 concentrations, and some other parameters. On this basis, the CO 2 concentration was solved by minimizing the error between the measured spectrum and the modeled one. J.R. Chen presented a novel multiple-wavelength IPDA inversion method using symmetrically measured LiDAR soundings, substantially reducing the error from spectral distortion and laser frequency drift in 2014 [25]. In 2016, Xiang et al. proposed another multiple-wavelength DIAL inversion method which is based on the weighted averaging of observation pairs [26], and they verified the method using a range-resolved, 11-wavelength CO 2 -DIAL in a laboratory environment. This method was further improved via using a better fitting kernel in 2017 [27]. We can conclude from these previous works that the technical feasibility of multiple-wavelength IPDA LiDAR to measure the CO 2 vertical gradient has been established [28]. In this work, we aim to discuss the feasibility of obtaining gradients of xCO 2 through our recommended method by multiple-wavelength IPDA LiDAR [13].
In the proposed framework, we sliced the atmosphere into several layers. Then, we constructed simultaneous equations using IPDA equations of multiple wavelengths. On this basis, the simultaneous Remote Sens. 2020, 12, 2395 3 of 20 equations were solved using the constrained linear least-squares technique, obtaining sliced xCO 2 . The linear least-squares technique is widely used as a standard approach to approximate the solution of over-determined problems [29]. However, it searches for a mathematically optimal solution without consideration of physical constraints [30]. The spatial distribution of CO 2 has its own special features, which can impose physical constraints on inversions [31]. The ordinary linear least-squares technique would obtain results with a high accuracy in terms of mathematics but with evident errors in terms of physics. We thus consider using a constrained linear least-squares method to obtain results with physical meanings. Such a technique has been used to resolve problems regarding remote sensing [32] but has not yet been used to deal with data of IPDA LiDAR. The sign of the xCO 2 gradient between the ABL and FT would be auxiliary data to help identify the characteristic of carbon flux (source, neutral, and sink). Moreover, one can retrieve the regional carbon flux using xCO 2 gradients based on the equilibrium boundary layer theory [33].
The remaining sections are as follows: the principle of the proposed method is presented in Section 2, as well as relevant methodologies used in the experiments. The performance of the proposed method is evaluated in Section 3 in terms of the SNR, the number of wavelengths, the number of layers, and the ABLH (the height of the atmospheric boundary layer). Moreover, we tested the ability of the proposed framework to retrieve gradients of xCO 2 . In Section 4, we discuss possible effects of water vapor on retrieving sliced xCO 2 using the proposed method. Finally, the major findings of this work are summarized in Section 5.

IPDA Inversion Algorithm
The first task is to introduce the hard-target LiDAR equation, as described in (1) [34]. Here, E λ is the received energy for a specified wavelength λ and for a given distance. Q is the surface reflectance. ξ i is the total instrument efficiency for wavelength i. In addition, r is the detection range (distance from target to receiver). E 0 is the output energy of a single laser pulse. A is the receiver area; N CO2 is the density value of CO 2 . T atm is the optical transmission due to other atmospheric constituents. σ(λ) is the absorption cross-section of CO 2 at wavelength λ, which is related to the atmospheric pressure and temperature [34,35].
On the basis of the LiDAR equation, the dual-wavelength IPDA technique utilizes a differential process to calculate absorption optical depth due to CO 2 . XCO 2 is retrieved using Equation (2): XCO 2 represents CO 2 concentration, and p surface and p toa represent pressure at the surface and top of the atmosphere, respectively. E off is the received energy of λ off , E on is the received energy of λ on , while E 0 on represents the emitted energy of λ on , E 0 o f f represents the emitted energy of λ off . The wavelength located in the strong absorption position of the CO 2 molecule is called the on-line wavelength (λ on ). The wavelength located in the weak absorption spectrum is considered the off-line wavelength (λ off ). WF (λ, p, T) is the weighting function of CO 2 , which is related to pressure and temperature. WF could be determined by Equation (3). where g denotes the intensity of the local gravitational field of the Earth, and m dry and m H2O represent the average masses of dry air and water vapor molecules, respectively, and q H2O is the dry-air volume mixing ratio of water vapor.

Multi-Wavelength Inversion Framework
The classic IPDA algorithm utilizes the observations of two wavelengths for the differential process. Given n+1 wavelengths across a single absorption line of CO 2 , we would have n differential equations. For each wavelength except the off-line wavelength, a classic dual-wavelength IPDA equation, as expressed in Equation (4), can be established. A series of such equations can be then expressed by Equation (5). ln represents xCO 2 of the i th layer, i = 1,2,...m. WF j is the total WF calculated by wavelengths λ on(j) and λ off , WF i j is i th layer WF calculated by the wavelengths of λ on(j) , and λ off, E 0 on( j) represents the emitted energy of λ on(j) . E on( j) represents the received energy of λ on(j) , and ξ j represents the error of total DAOD, calculated by the wavelengths of λ on(j) and λ off , j = 1,2 . . . ,n.
∆p i represents the difference of pressure in layer i, i=1, 2, . . . , m. Consequently, the parameters of W and obs would be expressed by Equations (6)- (8).
, ..., ln where ξ is the detection random error, m is the number of layers in the atmosphere, and n is the number of on-line wavelengths. In Equation (9), obs is an n-dimensional observation vector, W is an n × m-dimensional coefficient matrix, x is a one-dimensional unknown parameter vector, and ξ is a residual vector of observations. On this basis, the constraint linear least-squares technique can be used to minimize ξ. We need to solve Equations (6)-(9) to determine x i CO2 . To solve Equation (9), an implied condition is n > m.
Remote Sens. 2020, 12, 2395 5 of 20 In Equation (10), A is an m×m-dimensional coefficient matrix and b is a one-dimensional vector. The inequation A·x ≤ b can be used to describe the correlation of xCO 2 of different layers. For example, if we had a-priori CO 2 profiles describing the vertical structure of xCO 2 , this information can be imposed on our inversion via A·x ≤ b. The equation Aeq·x = beq is used to inverse slice xCO 2 , with a higher accuracy when the total column concentration and WF are well known. Aeq represents the WF of CO 2 in layers, x is layer xCO 2 , and beq is the DAOD of the total layer of atmosphere. l ≤ x ≤ u is used to set lower/upper boundaries to x by specifying l and u. Varieties of auxiliary data can be also integrated into the proposed frameworks, such as in-situ surface xCO 2 measurements and column average xCO 2 retrievals from ground-based and spaceborne Fourier transform infrared spectrometers. The principle of constrained Linear Least-Squares Technique could be found in Appendices A and B. If there is no a-priori knowledge of the inversion process, Equation (7) would be solved using the general least-squares method.
In the experiments of this work, we imposed the following constraints on the inversion process: (a) xCO 2 of each layer should be less than 425 ppm (b) xCO 2 of each layer should be greater than 370 ppm (c) A · x ≤ b A and b in Equation (10) were determined by Equation (11).
x is the sliced xCO 2 in each layer, and XCO 2 max is the upper boundary of XCO 2 . We utilized Xiang's method [26] to obtain the first guess of XCO 2 . It is worth noting that we only set an initial value for the column weighted xCO 2 as a constraint in the inversion process. All sliced xCO 2 are not assigned to certain initial values. The value of b in Equation (11) is based on the principle of double standard deviation limitation of probability distributions. Xiang's method used weighted averaging of XCO 2 retrieved by different wavelengths to obtain a more accurate XCO 2 , but this method cannot obtain sliced xCO2 using IPDA LiDAR. In this work, we used the first guess of XCO 2 as one of the constraints in solving equations using the linear least-squares technique. This means that this is never a compulsory procedure. One can use XCO2 provided by other equipment or retrievals of a traditional dual-wavelength IPDA inversion instead. N air (layer i ), (i = 1,2,3 . . . ,m) represents the number density of air molecules in each layer; N air (total) represents the number density of air molecules in the total atmosphere. SNR (T, λcentral) represents the signal-to-noise ratio of the DAOD measured by λ central , which is further described in Section 3.1. A simple explanation of constraint c herein is that differences between retrieved XCO2 using traditional IPDA algorithms and the proposed sliced xCO 2 inversion algorithm are limited to a certain extent, which is related to the SNR of signals.
Regarding settings of the layer number, m is set to 2, 3, and 4 in this work. Here, n is the number of wavelengths of λ on(j) . A key problem is to determine a rational range for each layer. From a physical point of view, the sensitivity to the vertical profile of CO 2 mainly comes from the differences in weighting functions between channels, which come from the fact that the absorption cross-section depends on the pressure. Therefore, a larger m means more unknown parameters to be solved and thus calls for more independent information. As the number of wavelengths increases, the wavelength spacing decreases accordingly, which exacerbates the correlation between the different detection channels. Thus, an increase in m will result in an extremely high increase in n, which in turn will result in a significant decrease in the energy that can be allocated to each channel. Therefore, we do not think it is wise to set an oversized m. A detailed analysis of the effect of m can be seen in the Discussion. Referring to measurements of the vertical profiles of CO 2 [36,37], there is an evident gradient of xCO 2 between ABL and FT, which could be of great significance for estimating CO 2 fluxes. To obtain the xCO 2 gradient there, we need to at least know xCO 2 in ABL and FT. Consequently, we set m to 3 in this work and sliced the atmosphere into ABL, FT, and STA (top of the atmosphere) as the default configuration. In Section 3.3, we further discuss the effect of the number of layers. Settings of xCO 2 are demonstrated in Table 1 as a part of the methodology, referring to the NASA LaRC Airborne Science Data for Atmospheric Composition (https://www-air.larc.nasa.gov/cgi-bin/ArcView/ascends.2017). Atmospheric profiles of pressure and temperature were preset using the US standard atmosphere model. The surface temperature and pressure are 297 K and 101.32 hPa, respectively. In the future of real applications, the number of layers and ranges of layers can be set according to location.

Effect of Signal-to-Noise Ratio
SNRs of received signals of wavelength λ can be expressed by Equation (12).
where E λ represents the received energy of λ, M is the internal gain factor of the detector, R denotes the response of the detector, B is the electrical bandwidth, e is the elementary charge, F is the excess noise factor of the detector (which accounts for additional noise due to the internal amplification statistics), and i 2 D is the dark current noise density. E back is the energy of solar background radiance. One can conclude from Equation (12) that SNR λ is highly related to the hardware configuration. Moreover, some environmental factors can also affect SNR λ , such as the aerosol optical depth and the surface reflectance which determines P λ . Moreover, the laser wavelength is also a critical parameter determining SNR λ because P λ varies significantly with λ due to CO 2 absorption.
δ (λ,CO2) is the measured optical depth due to CO 2 absorption on a certain wavelength of λ, which can be calculated by Equation (13). The relationship between SNR λ and the relative random error of δ (λ,CO2) can be expressed by Equation (14). N is the average number of the shots, which depends on the average time and emitted frequency of the future satellite IPDA. We further defined SNR (T,λ) using Equation (15), where ∆δ (λ,CO2) is the random error of δ (λ,CO2) . SNR (T,λ) is very different from SNR λ , The former is used to describe the SNR of measured δ (λ,CO2) , which is proportional to XCO 2 , while the latter is used to describe the SNR of received laser signals. Given the specific configurations of hardware and detection environments, both SNR λ and δ (λ,CO2) are negatively correlated with random errors of retrievals, as shown in Equation (14). However, SNR λ decreases with the increasing δ (λ,CO2) . Therefore, it is unwise to establish a relationship between the detection error and SNR λ . Compared with SNR λ , SNR (T,λ) can reflect the accuracy of retrievals more intuitively. Hence, we define SNR (T,λ) herein and correlate it with the detection error in the following experiments.
SNR (T,λcentral) , where λ central is the central wavelength of on-line wavelengths, would be better than 25 based on the technology of the current hardware [38]. In this group of simulations, the preset vertical profile of xCO 2 is shown in Table 1. The number of wavelengths used in this group of simulations is 22. If SNR (T,λcentral) is certain, the SNRs (T,λi) (i = 1,2,3 . . . m, m = 21) would be calculated by Equation (14).
In order to evaluate the influence of SNR (T,λcentral) on the final results, we performed 10,000 simulations to evaluate the performance of our recommended method under different SNR (T,λcentral) (20)(21)(22)(23)(24)(25)(26) values. The number of wavelengths is 22 and the assumed ABLH is 1.5 km. Then, the xCO 2 of three sliced layers was retrieved, as shown in Figure 1. Atmospheric profiles of pressure and temperature were preset using the US standard atmosphere model. The temperature and pressure of the surface are 297 K and 101.32 hPa, respectively. In the subsequent experiments of other subsections, we used the same settings of the atmosphere, unless we have declared otherwise.
Remote Sens. 2020, 11, x FOR PEER REVIEW 7 of 20 In order to evaluate the influence of SNR(T,λcentral) on the final results, we performed 10,000 simulations to evaluate the performance of our recommended method under different SNR(T,λcentral) (20)(21)(22)(23)(24)(25)(26) values. The number of wavelengths is 22 and the assumed ABLH is 1.5 km. Then, the xCO2 of three sliced layers was retrieved, as shown in Figure 1. Atmospheric profiles of pressure and temperature were preset using the US standard atmosphere model. The temperature and pressure of the surface are 297 K and 101.32 hPa, respectively. In the subsequent experiments of other subsections, we used the same settings of the atmosphere, unless we have declared otherwise.  As is shown in Figure 1, the mean error (ME) decreases with the increasing SNR(T,λcentral). It is worth nothing that MEs would be −0.98 ppm, 0.27 ppm, and −0.15 ppm in ABL, FT, and STA, respectively, when SNR(T,λcentral) is 25 dB. STD has a negative relationship with the SNR(T,λcentral) in three layers. However, both ME and STD performance slightly change when SNR(T,λcentral) is larger than 25 db. Hence, we set SNR(T,λcentral) to 25 db in the following simulations. Moreover, it is necessary to discuss the feasibility of an SNR(T,λcentral) of 25 db. We have used the performance simulator and configurations of China's CO2-IPDALiDAR described in our previous work [38] to calculate the SNR and found that it is higher than 25 db. Moreover, according to the flight test reported by Zhu et al., an SNR of 25 db is feasible in real conditions [39].

Effect of the Number of Wavelengths and Selecting Modes
A series of experiments was carried out to explore the effect of the number of wavelengths, in which the central wavelength was set to 1572.335 nm. Except for the number of wavelengths, two types of sampling mode were tested to determine whether the strategy of wavelength sampling would exert influences on outcomes. In two different sampling modes, we set the same wavelength step. With such a configuration, we can determine what the dominant factor is, the wavelength number or the wavelength step, in determining the performance of the proposed multi-wavelength method. Figure 2 demonstrates an example of symmetric sampling with 22 wavelengths (off-line wavelength is not shown), using both solid and dotted lines. An example of unilateral sampling with 12 wavelengths is also shown in this figure, using only the solid lines.
In the next experiment, we wanted to investigate whether different sampling modes would cause differences in the performance of the proposed method. The center wavelength of a CO2 spectrum would move towards the shorter wavelength with the decreasing pressure. Such a phenomenon is called pressure-induced shift [40]. Therefore, symmetric sampling in the absolute sense is not possible. This effect has been taken into consideration in the simulations to include possible errors. In our experiments, two sample modes were used: symmetric distribution sampling and unilateral distribution sampling. The SNR(T,λcentral) was set to 25 dB and the ABLH is 1.5 km. For the simulations in this section, the setting of the number of wavelengths is shown in Table 2 and the detailed setting of the layers and concentration of CO2 in each layer is shown in Table 1. As is shown in Figure 1, the mean error (ME) decreases with the increasing SNR (T,λcentral) . It is worth nothing that MEs would be −0.98 ppm, 0.27 ppm, and −0.15 ppm in ABL, FT, and STA, respectively, when SNR (T,λcentral) is 25 dB. STD has a negative relationship with the SNR (T,λcentral) in three layers. However, both ME and STD performance slightly change when SNR (T,λcentral) is larger than 25 db. Hence, we set SNR (T,λcentral) to 25 db in the following simulations. Moreover, it is necessary to discuss the feasibility of an SNR (T,λcentral) of 25 db. We have used the performance simulator and configurations of China's CO 2 -IPDALiDAR described in our previous work [38] to calculate the SNR and found that it is higher than 25 db. Moreover, according to the flight test reported by Zhu et al., an SNR of 25 db is feasible in real conditions [39].

Effect of the Number of Wavelengths and Selecting Modes
A series of experiments was carried out to explore the effect of the number of wavelengths, in which the central wavelength was set to 1572.335 nm. Except for the number of wavelengths, two types of sampling mode were tested to determine whether the strategy of wavelength sampling would exert influences on outcomes. In two different sampling modes, we set the same wavelength step. With such a configuration, we can determine what the dominant factor is, the wavelength number or the wavelength step, in determining the performance of the proposed multi-wavelength method. Figure 2 demonstrates an example of symmetric sampling with 22 wavelengths (off-line wavelength is not shown), using both solid and dotted lines. An example of unilateral sampling with 12 wavelengths is also shown in this figure, using only the solid lines.
In the next experiment, we wanted to investigate whether different sampling modes would cause differences in the performance of the proposed method. The center wavelength of a CO 2 spectrum would move towards the shorter wavelength with the decreasing pressure. Such a phenomenon is called pressure-induced shift [40]. Therefore, symmetric sampling in the absolute sense is not possible. This effect has been taken into consideration in the simulations to include possible errors. In our experiments, two sample modes were used: symmetric distribution sampling and unilateral distribution sampling. The SNR (T,λcentral) was set to 25 dB and the ABLH is 1.5 km. For the simulations in this section, the setting of the number of wavelengths is shown in Table 2 and the detailed setting of the layers and concentration of CO 2 in each layer is shown in Table 1.      Figure 3 shows that both ME and STD decrease with the increasing number of wavelengths in three layers. We also found that the number of wavelengths, rather than the wavelength step, is the dominant factor in determining the performance of inversions. Two sampling modes yield very similar performance with the same number of wavelengths. With 22 wavelengths, MEs of retrievals are −0.98 ppm, 0.27 ppm, and −0.15 ppm in ABL, FT, and STA, respectively. We think that 22 wavelengths are adequate and appropriate to obtain retrievals of each layer with an error of less than 1 ppm. Consequently, we set the number of wavelengths to 22 in the following simulations.
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 20    Figure 3 shows that both ME and STD decrease with the increasing number of wavelengths in three layers. We also found that the number of wavelengths, rather than the wavelength step, is the dominant factor in determining the performance of inversions. Two sampling modes yield very similar performance with the same number of wavelengths. With 22 wavelengths, MEs of retrievals are −0.98 ppm, 0.27 ppm, and −0.15 ppm in ABL, FT, and STA, respectively. We think that 22 wavelengths are adequate and appropriate to obtain retrievals of each layer with an error of less than 1 ppm. Consequently, we set the number of wavelengths to 22 in the following simulations.
(a)  Figure 3a-c, red circles and squares represent the ME and STD obtained using the symmetric sampling mode. Black circles and squares represent the ME and STD obtained using the unilateral sampling mode.

Effect of the Number of Layers
In order to determine the optimal number of sliced layers, we discussed different options for the number of slicing layers, including two, three, and four layers, into the inversion. In this experiment, we changed the default setting of the CO2 profile into a four-layer profile, as shown in Table 3. The simulated signals were then generated on this basis and were used to obtain sliced xCO2 of two-layer and three-layer. We also calculated the pseudo-truth of sliced xCO2 using weighted averaging in terms of sums of atmospheric pressures in different layers for two-layer and three-layer inversions. Then, the mean error was calculated as the difference between retrieved XCO2 and the corresponding pseudo-truths. The results of the simulations are shown in Table 3. In this section, the SNR(T,λcentral) was set to 25 dB, the setting of the number of wavelengths was 22, and the detailed setting of the layers and concentration of CO2 in each layer shown in Table 1.  Figure 3a-c, red circles and squares represent the ME and STD obtained using the symmetric sampling mode. Black circles and squares represent the ME and STD obtained using the unilateral sampling mode.

Effect of the Number of Layers
In order to determine the optimal number of sliced layers, we discussed different options for the number of slicing layers, including two, three, and four layers, into the inversion. In this experiment, we changed the default setting of the CO 2 profile into a four-layer profile, as shown in Table 3. The simulated signals were then generated on this basis and were used to obtain sliced xCO 2 of two-layer and three-layer. We also calculated the pseudo-truth of sliced xCO 2 using weighted averaging in terms of sums of atmospheric pressures in different layers for two-layer and three-layer inversions. Then, the mean error was calculated as the difference between retrieved XCO 2 and the corresponding pseudo-truths. The results of the simulations are shown in Table 3. In this section, the SNR (T,λcentral) was set to 25 dB, the setting of the number of wavelengths was 22, and the detailed setting of the layers and concentration of CO 2 in each layer shown in Table 1. Table 3. Performance evaluation of different number of sliced atmospheric layers. ABL B represents the layer below ABLH, ABL U represents the layer upper ABLH, LFT represents the lower free troposphere, and UFT represents the upper free troposphere. The column, xCO 2 , represents preset xCO 2 in each layer. Table 3 shows that the number of sliced layers have evident influence on the performance of our retrieval model. Both the accuracy and the precision of retrievals decrease with the increasing number of sliced layers. The scheme of two layers yields results with the highest accuracy and precision, in which MEs are −0.46 ppm and 0.11 ppm, respectively, whose corresponding STDs are also the lowest among the three schemes. Moreover, these results exhibit that both the accuracy and the precision are related to the range of each layer. The proposed method can obtain xCO 2 with errors of less than 1 ppm for each layer when the number of sliced layers is less than three.

Number of Layers Range (km) Pseudo-Truth of xCO 2 (ppm) Mean Error (ppm) STD (ppm)
More sliced layers help us to gain more information on the vertical mixing of CO 2 but lead to lower accuracy and precision. Hence, there is a compromise in terms of the selection of an appropriate number of sliced layers. Near-surface CO 2 has a strong relationship with human activity and ecosystems. The CO 2 gradient around PBL is meaningful for us to quantitatively calculate carbon fluxes. If we sliced the atmosphere into two layers, the concentration of CO 2 in ABL would be acquired with high accuracy, but the range of ABL U layer is too large compared with the range of ABL B . It is thus not reasonable to calculate the CO 2 gradient around the ABL. If we sliced the atmosphere into four layers, the accuracy of retrievals in the ABL and FT would be insufficient to obtain accurate estimates of CO 2 gradients between ABL and FT. Therefore, three sliced layers would be the most appropriate.

Effect of the Atmospheric Boundary Layer Height (ABLH)
The ABL is the lowest layer of the atmosphere, which is directly affected by the Earth's surface and responds to radiative forcing. LiDAR is capable of measuring the ABLH with high accuracy. ABLH has a relationship with the location, the landform, and meteorological conditions. To evaluate the performance of the proposed retrieval method with respect to different heights of PBL, we carried out inversion simulations under different settings of the PBLH, namely 1.0 km, 1.5 km, and 2.0 km. In this section, the SNR (T,λcentral) was set to 25 dB, and the setting of the number of wavelengths was 22. Table 4 shows that the accuracy of the retrieved xCO 2 in ABL would increase with the ascending ABLH. Especially for the xCO 2 in ABL, the mean error of xCO 2 in ABL was −1.76 ppm when ABLH was set to 1 km. Along with the ME, STD also falls with the increasing ABLH, implying that the precision of xCO 2 is also affected by the ABLH. We believe that the mechanism by which the ABLH affects the performance of the proposed method can be explained by Equation (9). A thicker ABL means a longer optical path, which leads to a higher DAOD. This would produce better adjustment in the ABL layer to promise the accuracy for xCO 2 in ABL during the process of our retrieval model. In this subsection, experiments are designed to explore the ability of the proposed multi-wavelength method to retrieve gradients of CO 2 under different circumstances. The ultimate goal of any space-based CO 2 sensor is to map the dynamic distribution of carbon sources and sinks. The vertical gradient of CO 2 is a critical parameter for distinguishing carbon sources and sinks. With enough accuracy, the vertical gradient of CO 2 can provide the potential ability to estimate carbon fluxes quantitatively. Three typical CO 2 profiles are assumed, as shown in Table 5, to simulate the possible characteristics of different carbon fluxes. In this series of experiments, we assumed three scenarios according to differences between X CO 2 in ABL and FT and named them as shown in Table 5. It is worth noting that the vertical profiles of CO 2 are determined by many factors. In some extreme cases, a lower xCO 2 in ABL with respect to FT could be witnessed for CO 2 source regions and vice versa. Herein, we utilized the potential source, sink, and neutral to name the three scenarios, purely for simplification. In this group of experiments, the number of wavelengths was set to 22, the SNR (T,λcentral) was set to 25 dB, and the wavelength range was set to 1572.305-1572.365 nm.  Figure 4 shows the ME and STD of retrievals in different layers for three scenarios. Each subfigure shows the performance of the proposed method in a specific layer and different colors are used to represent different characteristics of CO 2 vertical profiles, namely a potential source, a potential sink, and a potential neutral area. The MEs of xCO 2 retrievals in ABL are −0.98 ppm, −0.31 ppm, and −0.11 ppm for a potential carbon source, a potential sink, and a potential neutral region, respectively. Compared with retrievals in ABL, retrievals in the other two layers are much more accurate, yielding an ME of no more than 0.34 ppm regardless of carbon characteristics. The STD never exceeds 0.73 ppm in any case. These results indicate that the proposed method can obtain accurate and precise results regardless of carbon characteristics.  The mean errors of xCO2 in ABL and xCO2 in FT would be anti-correlated (one is positive while the other negative) for the source and sink cases because, in the retrieval model, it would adjust the value of DAOD in each layer according to the weight matrix W in the manuscript. Additionally, it would promise that the total DAOD is certain. The DAOD in STA is larger than that of ABL and FT. Hence, the main adjustments of DAOD occur in ABL and FT, and mean errors in DAODs in ABL and FT would have a neutralized relationship. Meanwhile, it is also worth noting that the performance of the proposed method varies with different characteristics of CO2 vertical profiles. Figure 5 shows the ME and STD of retrieved xCO2 gradients, including the xCO2 gradient between ABL and FT and between FT and STA. For a potential carbon source, the ME and STD of the The mean errors of xCO 2 in ABL and xCO 2 in FT would be anti-correlated (one is positive while the other negative) for the source and sink cases because, in the retrieval model, it would adjust the value of DAOD in each layer according to the weight matrix W in the manuscript. Additionally, it would promise that the total DAOD is certain. The DAOD in STA is larger than that of ABL and FT. Hence, the main adjustments of DAOD occur in ABL and FT, and mean errors in DAODs in ABL and FT would have a neutralized relationship. Meanwhile, it is also worth noting that the performance of the proposed method varies with different characteristics of CO 2 vertical profiles. Figure 5 shows the ME and STD of retrieved xCO 2 gradients, including the xCO 2 gradient between ABL and FT and between FT and STA. For a potential carbon source, the ME and STD of the retrieved xCO 2 gradient (xCO 2ABL -xCO 2FT ) are −1.25 ppm and 1.05 ppm, xCO 2ABL is xCO 2 in ABL, and xCO 2FT is xCO 2 in FT. For a potential carbon sink, the ME and STD of the retrieved xCO 2 gradient (xCO 2ABL − xCO 2FT ) are −1.03 ppm and 1.02 ppm. For a potential neutral region, the ME and STD of the retrieved xCO 2 gradient (xCO 2ABL -xCO 2FT ) are −0.17 ppm and 0.99 ppm. Consequently, the accuracy and precision of the retrieved xCO 2 gradients are also influenced by the characteristics of the vertical CO 2 profiles. Considering that the true carbon gradient (xCO 2ABL -xCO 2FT ) is set to 8 and −12 ppm, its accuracy would be acceptable for its largest uncertainty of only around 16%. Figure 4. The retrievals of sliced xCO2 with settings of the potential areas of carbon source (red), carbon sink (green), and neutral (blue) characteristics, corresponding to ABL, FT, and STA layers, respectively.
The mean errors of xCO2 in ABL and xCO2 in FT would be anti-correlated (one is positive while the other negative) for the source and sink cases because, in the retrieval model, it would adjust the value of DAOD in each layer according to the weight matrix W in the manuscript. Additionally, it would promise that the total DAOD is certain. The DAOD in STA is larger than that of ABL and FT. Hence, the main adjustments of DAOD occur in ABL and FT, and mean errors in DAODs in ABL and FT would have a neutralized relationship. Meanwhile, it is also worth noting that the performance of the proposed method varies with different characteristics of CO2 vertical profiles. Figure 5 shows the ME and STD of retrieved xCO2 gradients, including the xCO2 gradient between ABL and FT and between FT and STA. For a potential carbon source, the ME and STD of the retrieved xCO2 gradient (xCO2ABL-xCO2FT) are −1.25 ppm and 1.05 ppm, xCO2ABL is xCO2 in ABL, and xCO2FT is xCO2 in FT. For a potential carbon sink, the ME and STD of the retrieved xCO2 gradient (xCO2ABL-xCO2FT) are -1.03 ppm and 1.02 ppm. For a potential neutral region, the ME and STD of the retrieved xCO2 gradient (xCO2ABL-xCO2FT) are −0.17 ppm and 0.99 ppm. Consequently, the accuracy and precision of the retrieved xCO2 gradients are also influenced by the characteristics of the vertical CO2 profiles. Considering that the true carbon gradient (xCO2ABL-xCO2FT) is set to 8 and −12 ppm, its accuracy would be acceptable for its largest uncertainty of only around 16%. Figure 5. The retrievals of XCO2 gradient with settings of the potential areas of carbon source (red), carbon sink (green), and neutral characteristics (blue) corresponding to the carbon gradient between ABL and FT and the carbon gradient between FT and STA. Figure 5. The retrievals of X CO 2 gradient with settings of the potential areas of carbon source (red), carbon sink (green), and neutral characteristics (blue) corresponding to the carbon gradient between ABL and FT and the carbon gradient between FT and STA.
As for the xCO 2 gradient between FT and STA, the uncertainty of the retrieved xCO 2 gradients is less than 7.5% in all three scenarios. Moreover, both the accuracy and precision of the retrieved xCO 2 gradients between FT and STA are higher than the carbon gradients between ABL and FT. Referring to Figure 4, we believe that the lower accuracy and precision of the xCO 2 retrievals in ABL are responsible for such a phenomenon.

Discussion
The water vapor is widely regarded as the most significant interference in measuring CO 2 by using the IPDA LiDAR. In the previous experiments, we ignored interferences from water vapor, which may result in misleading conclusions. Hence, we try to determine the influence on CO 2 concentration derived from the errors in water vapor information in this subsection.

Interference of Water Vapor in the Absorb Spectrum
For a CO 2 -IPDA LiDAR working near 1.6 µm, there are three commonly used candidate lines, namely R14, R16, and R18 in the 30012←00001 band of CO 2 . We select all lines of water between 6358 and 6362 cm −1 and then set a cutoff of 1E-28 cm/molecule to the line intensity to screen useful  [41]. Figure 6 shows cross-sections of CO 2 and water vapor in the wave number span of 6358 cm −1 to 6362 cm −1 . Since line intensities of water vapor are smaller than those of CO 2 by several orders, cross-sections of two molecules mark different vertical axes. Moreover, line intensities of water vapor are diverse in terms of absorption lines. Therefore, cross-sections of water vapor are shown in a logarithmic scale. Different isotopologues of the water vapor are shown in different colors, helping readers to understand the effect of isotopologues on the water vapor. The natural abundances of HDO, H 2 18 O, and H 2 O are set to 0.00031, 0.002, and 0.998, respectively. and 6362 cm −1 and then set a cutoff of 1E-28 cm/molecule to the line intensity to screen useful lines. After this, there are 24 lines left, in which there are 9, 12, and 3 lines of H2O, HDO, and H2 18 O, respectively. Spectroscopic parameters of CO2 and water vapor used in this experiment would be acquired from HITRAN 2016 [41]. Figure 6 shows cross-sections of CO2 and water vapor in the wave number span of 6358 cm -1 to 6362 cm -1 . Since line intensities of water vapor are smaller than those of CO2 by several orders, crosssections of two molecules mark different vertical axes. Moreover, line intensities of water vapor are diverse in terms of absorption lines. Therefore, cross-sections of water vapor are shown in a logarithmic scale. Different isotopologues of the water vapor are shown in different colors, helping readers to understand the effect of isotopologues on the water vapor. The natural abundances of HDO, H2 18 O, and H2O are set to 0.00031, 0.002, and 0.998, respectively. Figure 6 shows that the majority of water interferences are related to HDO. One line of H2O is located at the right wing of the R18 line of CO2, whose central wave number is 6361.250 cm −1 . From the perspective of the absorption of a single molecule, the water vapor would cause evident errors because the cross-section of water vapor is about 0.1-1% of that of CO2. The total error budget for a space-borne XCO2 measurement is merely 0.3% [42]. Therefore, the interference from water has been regarded as a critical error source for a long time.   Figure 6 shows that the majority of water interferences are related to HDO. One line of H 2 O is located at the right wing of the R18 line of CO 2 , whose central wave number is 6361.250 cm −1 . From the perspective of the absorption of a single molecule, the water vapor would cause evident errors because the cross-section of water vapor is about 0.1-1% of that of CO 2 . The total error budget for a space-borne XCO 2 measurement is merely 0.3% [42]. Therefore, the interference from water has been regarded as a critical error source for a long time.
However, in the application of measuring atmospheric CO 2 in natural conditions, we should focus on the measured absorption optical depth due to trace gases but not cross-sections. Figure 7 shows the simulated absorption optical depth of carbon dioxide and water vapor from the top of the atmosphere to the surface. In this experiment, we used US standard atmosphere 1976 to provide atmospheric parameters such as T, P, air density, and g. Concentrations of carbon dioxide and water vapor are set to constants 400 ppm and 15 g/kg, respectively. It is worth noting that the concentration of water vapor would decrease quickly with the altitude. Consequently, we have already overestimated significantly the effect of water vapor in this experiment. Figure 7a also shows that interferences from water vapor are negligible for R14 and R16 in terms of the measured absorption optical depth. For R18, there are evident interferences at the right wing. However, interferences are negligible at the left wing of R18. Consequently, one can select appropriate on-line and off-line wavelengths to avoid interference from the water vapor. For the multi-wavelength inversion framework, we need to scan a whole spectrum. Figure 7b shows a detailed comparison of absorption optical depth (AOD) in R16 line. The ratio of AOD_water and AOD_CO 2 fluctuates between 0.00001% and 0.02%. If we further confine the range of wavelengths to 1572.305-1572.365 nm, namely 6359.846-6360.089 cm −1 , the ratio would fluctuate between 0.0013% and 0.0019%, which accounts for merely less than 1% of the total error budget. Hence, the influence of water vapor on the absorption spectrum would be neglected. optical depth. For R18, there are evident interferences at the right wing. However, interferences are negligible at the left wing of R18. Consequently, one can select appropriate on-line and off-line wavelengths to avoid interference from the water vapor. For the multi-wavelength inversion framework, we need to scan a whole spectrum. Figure 7b shows a detailed comparison of absorption optical depth (AOD) in R16 line. The ratio of AOD_water and AOD_CO2 fluctuates between 0.00001% and 0.02%. If we further confine the range of wavelengths to 1572.305-1572.365 nm, namely 6359.846-6360.089 cm -1 , the ratio would fluctuate between 0.0013% and 0.0019%, which accounts for merely less than 1% of the total error budget. Hence, the influence of water vapor on the absorption spectrum would be neglected.   AOD_water_vapor (unitless) Figure 7. Absorption optical depths of water vapor and carbon dioxide. In both subfigures, the absorption optical depth (AOD) of CO 2 marks the right vertical axis while the AOD of water vapor marks the left axis. Concentration of CO 2 and water vapor are set as 400 ppm and 15 g/kg, respectively. Figure 7 (b) is the enlargement of Figure 7 (a) at the R16 line.

The Influence of Water Vapor on the Value of Dry-Air Mixing Ratio of CO 2
Though the last subsection has demonstrated that interference from water vapor can be ignored in terms of measuring DAOD, there is another mechanism by which additional errors would be introduced owing to the neglect of water vapor. The DAOD is proportional to the number density of CO 2 molecules. Water vapor would influence the value of DAOD of CO 2 . After obtaining the number density of CO 2 molecules, we need to transform it into the dry-air mixing ratio. That is to say, we need to exclude water vapor molecules from the wet air molecules. To consider such an effect, we simulated the vertical profile of water vapor using the concept of water vapor saturation pressure, as shown in Equation (16).
where W is the content of water in the atmosphere (g/kg), M (H2O) is the molar mass of water vapor, M (air) is the molar mass of air molecular, P is the atmospheric pressure, and p w represents the water vapor saturation pressure, which could be calculated using the Goff-Gratch method, recommended by WMO (World Meteorological Organization) [43]. However, in actual experiments, it would exclude errors in water vapor information regardless of whether they were obtained by reanalysis meteorological data or by actual measurements. Hence, we discussed the influence of different accuracies of water vapor on the final retrieval results, including 1%, −1%, 5%, and −5%, which could be calculated by Equations (17) to (18). The results are presented in Table 6.
where N air (i) is the number density of air molecules in layer i, i = 1,2,3. ratio (H2O) (i) is the dry-air mixing ratio of water vapor in layer i (i = 1,2 . . . n, n = 4). N CO2 (i) is the number density of CO 2 molecules in layer i, r i is the range interval of layer i. DAOD i CO2 and AOD i H20 represent the DAOD of CO 2 and DAOD of water vapor in layer i (i = 1,2 . . . n, n = 4). As is shown in Table 6, water vapor has a larger influence on xCO 2 in ABL than in FT and STA. Since the concentration of water vapor mostly exists in the lower atmosphere, there is nearly no water vapor when the altitude is higher than 20 km. However, this amount would be only 0.06 ppm even though the uncertainty of water vapor information is 5%. In conclusion, the influence of water vapor on sliced xCO 2 would be neglected based on our recommended retrieval model.

Conclusions
A multi-wavelength inversion framework for CO 2 -IPDA LiDAR is proposed to obtain sliced xCO 2 . In this framework, we slice the atmosphere into three layers, ABL, FT, and STA. Then, we construct multiple differential absorption equations using different on-line wavelengths. On this basis, the constrained linear least-squares technique can be used to obtain the optimum solution. Our method achieves a reliable performance when the SNR(T, λcentral) value is 25 and the wavelength number of 22. Further experiments confirmed the potential ability to distinguish different carbon characteristics by our recommended method. Even the largest error in the sliced xCO 2 would be less than 0.03%, while the accuracy of the xCO 2 gradient is within 1.25 ppm. Moreover, the results of sliced xCO 2 would be reasonable to neglect the influence of the present uncertainty of water vapor information. Therefore, if a multi-wavelength IPDA launched successfully, it would produce meaningful data to help us learn more about the vertical diffusion of CO 2 globally.
It could be proven that KKT conditions are sufficient and necessary conditions for the objective function to obtain the minimum value, based on Equation (6).

(A9)
Then, we could calculate the KKT conditions to acquire the extreme minimum value of f (x). 3. Both equality and inequality constraints like Equation (7) are calculated as follows: The solution of x should satisfy the following conditions, as in Equation (8): Then, the solution of x could be satisfied by the KKT conditions and minimize the Lagrangian function to obtain a feasible solution under these constraints.