Averaging bias correction for the future space-borne methane IPDA lidar mission MERLIN

The CNES (French Space Agency) and DLR (German Space Agency) project MERLIN is a future integrated path differential absorption (IPDA) lidar satellite mission that aims at measuring methane dry-air mixing ratio columns ($X_{{\mathrm{CH}}_{\mathrm{4}}})$ in order to improve surface flux estimates of this key greenhouse gas. To reach a 1 % relative random error on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ measurements, MERLIN signal processing performs an averaging of data over 50 km along the satellite trajectory. This article discusses how to process this horizontal averaging in order to avoid the bias caused by the non-linearity of the measurement equation and measurements affected by random noise and horizontal geophysical variability. Three averaging schemes are presented: averaging of columns of $X_{{\mathrm{CH}}_{\mathrm{4}}}$, averaging of columns of differential absorption optical depth (DAOD) and averaging of signals. The three schemes are affected both by statistical and geophysical biases that are discussed and compared, and correction algorithms are developed for the three schemes. These algorithms are tested and their biases are compared on modelled scenes from real satellite data. To achieve the accuracy requirements that are limited to 0.2 % relative systematic error (for a reference value of 1780 ppb), we recommend performing the averaging of signals corrected from the statistical bias due to the measurement noise and from the geophysical bias mainly due to variations of methane optical depth and surface reflectivity along the averaging track. The proposed method is compliant with the mission relative systematic error requirements dedicated to averaging algorithms of 0.06 % (±1 ppb for $X_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{1780}\mathrm{ppb}$) for all tested scenes and all tested ground reflectivity values.

1 Introduction

Methane (CH₄) is the second most important anthropogenic greenhouse gas after carbon dioxide (CO₂) (IPCC, 2013). Despite its key role in global warming, there are still uncertainties in the cause of the observed large fluctuations in the growth rate of atmospheric methane. Measuring atmospheric CH₄ concentration on a global scale with both high precision and accuracy is necessary to improve the surface flux estimate and thus develop the knowledge of the global methane cycle (Kirschke et al., 2013; Saunois et al., 2016).

The Methane Remote Sensing Lidar Mission (MERLIN – website: https://merlin.cnes.fr/, last access: 22 October 2018) is a joint French and German space mission with a launch scheduled for 2024 (Ehret et al., 2017). This mission is dedicated to the measurement of the integrated methane dry-air volume mixing ratio ($X_{{\mathrm{CH}}_{\mathrm{4}}})$. The German Space Agency (DLR) is responsible for the payload while the French Space Agency (CNES) is responsible for the platform (MYRIADE Evolution product line). The payload data processing centre is under CNES responsibility with significant contributions from DLR.

MERLIN's active measurement is based on a space-borne integrated path differential absorption (IPDA) lidar. Just like a differential absorption lidar (DIAL), MERLIN's IPDA lidar uses the difference in transmission between an online pulse with a frequency accurately set in the trough of several CH₄ absorption lines and an offline pulse whose wavelength has a negligible CH₄ absorption (Ehret et al., 2008). Furthermore, the two wavelengths are set close enough in such a way that the differential effects of any other interaction, excluding CH₄ absorption, are minimized. Figure 1 shows the positioning of the two wavelengths. However, unlike a DIAL, MERLIN's IPDA lidar provides the column content of a specific trace gas along the line of sight rather than the range-resolved profile of CH₄. This column-integrated methane mixing ratio can be retrieved from the return signals after they are backscattered on a hard target such as the surface of the Earth or dense clouds. The much higher backscatter signal from these targets allows for a system with a relatively small power-aperture product as compared to a DIAL, which has to rely on atmospheric backscatter.

Figure 1 Laser frequency positioning of the online and offline laser beams. The online frequency is positioned in the trough of one of the methane absorption line multiplets. The offline frequency is positioned so that the methane absorption is negligible.

The MERLIN measurements require a well-defined processing chain that ensures the final performance of the mission. The processing chain is divided into four levels. Level 0 (L0) consists of raw data (backscattered signals and auxiliary data), and level 1 (L1) processes the vertically resolved products and the differential absorption optical depth (DAOD) values for both individual calibrated signal shot pairs and for a horizontal averaging window. Level 2 (L2) computes the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ for both individual calibrated signal shot pairs and for a horizontal averaging window, additionally using operational analyses from numerical weather prediction (NWP) centres. Note that in the presence of clouds, two products are provided; the first one computes the average for clear-sky shots only, and the other one averages all shots. Finally, level 3 (L3) produces $X_{{\mathrm{CH}}_{\mathrm{4}}}$ maps using a Kalman filter approach (Chevallier et al., 2017).

To reach a usable precision, space-borne IPDA lidar missions often require an averaging of measurements along the orbit's ground track (Grant et al., 1988). This process of averaging data horizontally is a general concern for IPDA lidar missions. The data processing of the NASA Active Sensing of CO₂ Emissions over Nights, Days, and Seasons (ASCENDS) mission considers the averaging of multiple lidar measurements along track over 10 s (70 km with no gaps) to reduce the random error on the carbon dioxide mixing ratio: $X_{{\mathrm{CO}}_{\mathrm{2}}}$ (Jucks et al., 2015). Likewise, MERLIN's averaging process is included into L1 and L2 algorithms in order to reduce the relative random error (RRE) of DAOD and $X_{{\mathrm{CH}}_{\mathrm{4}}}$ (Fig. 2). For the MERLIN mission, measurements are averaged over a nominal window length of 50 km corresponding to about 150 shot pairs to reach an RRE of approximately 20 ppb.

Figure 2 Principle schematics of the MERLIN IPDA lidar measurement. The lidar emits two laser beams with slightly different wavelengths (λ_on and λ_off). Every measurement corresponds to the small fraction of the two laser beams – called online and offline signals – that are reflected by a “hard” target (Earth's surface, top of dense clouds) to the satellite receiver telescope. For clarity, the three averaging windows are represented with four measurements instead of 150. On every averaging window, geophysical parameters such as altitude (or scattering surface elevation when there are clouds) or reflectivity vary.

The non-linearity of the equation relating calibrated signals and DAOD in combination with both the statistical noise inherent to any measurement and the varying geophysical quantities (altitude, pressure, reflectivity) of the sounded scene increases the relative systematic error (RSE or bias) and impairs measurement accuracy. Werle et al. (1993) describe RRE reduction when averaging signals using the concept of Allan variance. Up to an optimal integration time, measurement variance reduces because the measurement is dominated by white noise. For greater integration times, the estimation is biased due to drifts inherent to the measurement systems. The aim of the present article is not to correct biases caused by real system drift but to correct biases that are caused by the non-linearity of the IPDA lidar measurement equation.

MERLIN must reach an unprecedented precision and accuracy on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ with a targeted RRE of 1 % (18 ppb). The targeted RSE must remain under 0.2 % (±3 ppb) in 68 % of cases; a limited budget of 0.06 % (±1 ppb for a $X_{{\mathrm{CH}}_{\mathrm{4}}}$ of 1780 ppb) is allocated to biases introduced by averaging algorithms with algorithms to correct these averaging biases. To reach the RRE target, levels 1 and 2 of MERLIN's signal processing requires a horizontal averaging of data over 50 km along track (Kiemle et al., 2011). Thus, the single shot online and offline random error is reduced by a factor of $\sqrt{\mathrm{150}}\approx \mathrm{12}$. For instance, for the typical target reflectivity (0.1), the online and offline signal-to-noise ratios (SNRs) are of the order of 6.1 and 16.5, respectively, and the equivalent SNRs for the averaged signals are, respectively, 79.6 and 197.2. This process greatly decreases the RRE of the $X_{{\mathrm{CH}}_{\mathrm{4}}}$.

Section 2 gives an overview of the IPDA measurement and MERLIN data processing. Section 3 defines and compares biases of several averaging schemes (described below) and suggests correction algorithms. Section 4 presents a comparative evaluation of these averaging schemes and associated bias correction procedures using modelled scenes based on real satellite data. And finally, in Sect. 5, the results of the simulation are described and a “best approach” algorithm (i.e. the least biased on tested scenes) is proposed for the MERLIN processing chain.

2 Overview of IPDA measurement and the MERLIN processing chain

2.1 IPDA measurement

MERLIN active measurement is based on a short-pulse IPDA lidar. The column content of methane between the satellite and a “hard” target (ground, vegetation, clouds, etc.) is retrieved by measuring the light that is reflected by the scattering surface, which is illuminated by two laser pulses with a slight wavelength difference. Figure 2 schematically shows the principle of the nadir-viewing space-borne lidar MERLIN. The pulse-pair repetition rate is 20 Hz, and the sampling distance is 350 m considering a ground spot velocity of about 7 km s⁻¹. The online and offline ground spots are separated by about 2 m, which is negligible compared to the ground diameter of the spots of about 100 m (90 % encircled energy). Shot pairs will be averaged over a 50 km window (about 150 shots pairs). The online wavelength λ_on (1645.552 nm; 6076.998 cm⁻¹) is positioned in the trough of one of the methane absorption line multiplets, whereas the offline wavelength λ_off (1645.846 nm; 6075.903 cm⁻¹), which serves as reference, is positioned such that the methane absorption is negligible (Fig. 1). Both wavelengths are close enough so that interactions with the ground and the atmosphere and instrumental response can be considered identical, notably for reflectivity, which is defined as the ratio of the power reflected toward the satellite receiver to that incident on the hard target. The difference is thus mostly sensitive to the difference in methane absorption.

2.2 MERLIN processing chain

When the offline and online radiation reach the photodetector (Avalanche Photo Diode), it is converted to photoelectrons and to an electrical current. The measured raw signal obtained is the sum of the lidar signal and a background signal that is produced by background light, detector dark current and electronic offset. This background signal must be estimated to be removed from the raw signal. In the presence of measurement noise, when the SNR is low, this process of background signal removal can lead to a negative estimated lidar signal.

For the sake of conciseness, we introduce for any variable X the notation X^on,off that interchangeably represents the online or offline variables X^on or X^off. By measuring the online and offline pulse energies denoted P^on,off (P^on or P^off, respectively), it is possible to compute the DAOD of methane and then retrieve $X_{{\mathrm{CH}}_{\mathrm{4}}}$ for the sounded column. We denote Q^on,off the measurements after normalization by the laser pulse energies, denoted E^on,off, and range r, which is the distance from the satellite to the reflective target:

$$\begin{array}{}\text{(1)}& Q^{\mathrm{on},\mathrm{off}}={\displaystyle \frac{P^{\mathrm{on},\mathrm{off}}\cdot r^{\mathrm{2}}}{E^{\mathrm{on},\mathrm{off}}}}.\end{array}$$

The quantity Q^on,off will be referred to as calibrated signals in the following sections of the present article. The DAOD used in this study, in which the contributions of other gases are neglected, is denoted δ and is computed as Eq. (2):

$$\begin{array}{}\text{(2)}& \mathit{\delta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{Q^{\mathrm{off}}}{Q^{\mathrm{on}}}}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\mathit{\tau }}^{\mathrm{2}}\right),\end{array}$$

where τ² is the relative two-way transmission. From δ, we can derive $X_{{\mathrm{CH}}_{\mathrm{4}}}$ from Eq. (3) (Ehret et al., 2008; Kiemle et al., 2011):

$$\begin{array}{}\text{(3)}& X_{{\mathrm{CH}}_{\mathrm{4}}}={\displaystyle \frac{\mathit{\delta }}{\mathrm{IWF}}}={\displaystyle \frac{{\int }_{p_{\mathrm{surf}}}^{\mathrm{0}}{\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}}}\left(p\right)\cdot \mathrm{WF}\left(p,T\right)\cdot \mathrm{d}p}{{\int }_{p\mathrm{surf}}^{\mathrm{0}}\mathrm{WF}\left(p,T\right)\cdot \mathrm{d}p}},\end{array}$$

where p_surf denotes the target pressure where the laser beam hits the ground, p and T are the pressure and temperature profiles, and ${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}}}\left(p\right)$ is the dry-air volume mixing ratio profile of methane. The weighting function WF(pT) describes the measurement sensitivity of $X_{{\mathrm{CH}}_{\mathrm{4}}}$ along the vertical, and IWF is the integrated weighting function of the column. These quantities are computed from meteorological and spectroscopic data, and the WF is given by the following equation:

$$\begin{array}{}\text{(4)}& \mathrm{WF}\left(p,T\right)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{on}}\left(p,T\right)-{\mathit{\sigma }}_{\mathrm{off}}\left(p,T\right)}{g\left(p\right)\cdot (M_{\mathrm{air}}+M_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\cdot {\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(p,T\right))}}.\end{array}$$

M_X denotes the molecular masses of the chemical species X, ${\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the dry-air volume mixing ratio of water vapour, g(p) stands for the acceleration of gravity (treated as altitude and hence p dependent), and here, σ_on,off are the cross sections for the online or offline wavelengths (not to be confused with the standard deviation notation σ used elsewhere in this article).

As previously mentioned, in order to reach the targeted 1 % relative random error on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ measurements, the signal processing of MERLIN requires a horizontal averaging of the data. However, we will show in next section that the non-linearity of Eq. (2) in combination with the measurement noise and the variability of the observed scene (surface elevation, reflectivity, meteorology) along the averaging window induces biases on the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$.

2.3 MERLIN measurement noise

As will be seen in the following sections, the noise that affects the measurement is one of the factors that induce the averaging bias on the retrieved methane mixing ratio. The noise originates from the detector noise, shot noise and speckle noise. In the case of MERLIN system, the dominant noise is the detector noise which is considered to be normal as it is mainly thermal noise. Then, due to the high number of photons within the signal (approximately 10³ for the dark current and lidar signal), the Poisson statistics approximates a shifted Gaussian distribution very well (central limit theorem). Furthermore, according to Kiemle et al. (2011), the laser speckle is not the dominant source of the statistical fluctuation and is even negligible thanks to the relatively large field of view and surface spot size. The normality of the noise on calibrated signals Q^on,off is also justified by real measurements (out of the scope of this paper). The noise model used to generate the simulated signals is based on MERLIN system parameters and is presented in the Appendix A.

3 Averaging schemes and bias correction: a theoretical approach

3.1 Definitions

In the following sections, we will use triangular bracket notation to denote the arithmetic sample mean $\langle Y\rangle =\frac{\mathrm{1}}{N}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{Y}_i$ of the quantity Y, and $\mathrm{\Delta }{Y}_i=Y_i-\langle Y\rangle$ will represent the deviation of the ith quantity to this arithmetic mean. By extension, when we use a weighted sample mean of the quantity Y, weighted by a quantity Z, we will denote it $\langle Y{\rangle }_{w\left[Z\right]}={\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{w}_i\left[Z\right]\cdot Y_i$, where $w_i\left[Z\right]=Z_i/{\sum }_{k=\mathrm{1}}^{N_{\mathrm{s}}}{Z}_k$ are the normalized weights used. The expected value of a random variable X will be denoted E[X], and the fact that X follows a normal distribution of mean value μ and variance σ² will be denoted $X\sim N\left(\mathit{\mu },{\mathit{\sigma }}^{\mathrm{2}}\right)$.

We are interested in the retrieval of the column-integrated methane concentration on a 50 km horizontal section along the satellite track. This quantity will be hereafter denoted ${\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{T}}$ (where T stands for target). The information that we can compute using the satellite measurements is the shot-by-shot $X_{{\mathrm{CH}}_{\mathrm{4}},i}$ (i is the shot index), which is related to the shot-by-shot volume mixing ratio of methane ${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i}\left(p\right)$ and the shot-by-shot weighting function WF_i(p) by Eq. (3). For the purpose of building the data processing chains, all the quantities must be described on a gridded model (vertical and horizontal discretization) of the atmosphere. This grid is composed of (N_l⋅N_s) cells where N_l is the number of vertical layers of the model and N_s is the number of shots that we want to average along the satellite path. To model the atmosphere, the pressure at the interface of each layer (at each N_l+1 levels) uses a hybrid sigma coordinate system and is denoted P_i,j. Note that the standard notation for indices will be kept consistent throughout this article. The first index (often denoted i) will represent the shot index and the second index (often denoted j) will represent the layer index (or level index). The term “level” stands for the vertical level in pressure units. The pressure thickness of every layer, denoted ΔP_i,j, is then derived from the pressure at every level.

The discrete form of Eq. (3) is

$$\begin{array}{}\text{(5)}& X_{{\mathrm{CH}}_{\mathrm{4}},i}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^{N_{\mathrm{l}}}{\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j}\cdot {\mathrm{WF}}_{i,j}\cdot \mathrm{\Delta }{\mathrm{P}}_{i,j}}{{\sum }_{j=\mathrm{1}}^{N_{\mathrm{l}}}{\mathrm{WF}}_{i,j}\cdot \mathrm{\Delta }{\mathrm{P}}_{i,j}}}.\end{array}$$

In order to define the average value ${\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{T}}$, we must define average values for the volume mixing ratio of methane and the weighting function. As the two quantities are intensive properties, it is necessary to multiply them by the pressure thickness to get the corresponding additive quantity. The average volume mixing ratio and the average weighting function of the jth layer are thus given by

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{{\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}}}}}_j={\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{\mathit{\pi }}_{i,j}\cdot {\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j},\text{(7)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathrm{WF}}}_j={\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{\mathit{\pi }}_{i,j}\cdot {\mathrm{WF}}_{i,j},\end{array}$$

where the weights are defined as

$$\begin{array}{}\text{(8)}& {\mathit{\pi }}_{i,j}={\displaystyle \frac{\mathrm{\Delta }{P}_{i,j}}{{\sum }_{k=\mathrm{1}}^{N_{\mathrm{s}}}\mathrm{\Delta }{P}_{k,j}}}.\end{array}$$

The pressure thickness, as an extensive property, is averaged arithmetically, and the average value is denoted ${\stackrel{\mathrm{‾}}{\mathrm{\Delta }P}}_j$. Then, we can define the average column-integrated methane concentration as

$$\begin{array}{}\text{(9)}& {\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{T}}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^{N_{\mathrm{l}}}{\stackrel{\mathrm{‾}}{\mathrm{vmr}}}_{{\mathrm{CH}}_{\mathrm{4}},j}\cdot {\stackrel{\mathrm{‾}}{\mathrm{WF}}}_j\cdot {\stackrel{\mathrm{‾}}{\mathrm{\Delta }P}}_j}{{\sum }_{j=\mathrm{1}}^{N_{\mathrm{l}}}{\stackrel{\mathrm{‾}}{\mathrm{WF}}}_j\cdot {\stackrel{\mathrm{‾}}{\mathrm{\Delta }P}}_j}}.\end{array}$$

3.2 Averaging schemes and types of biases

There are several ways to average the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ provided the shot-by-shot calibrated signals $Q_i^{\mathrm{on},\mathrm{off}}$. Table 1 presents four different averaging schemes: averaging of columns of $X_{{\mathrm{CH}}_{\mathrm{4}}}$ (AVX – first line of Table 1), averaging of columns of DAOD and IWF (AVD – second line of Table 1), averaging of signals (AVS – third line of Table 1) and averaging of quotients (AVQ – fourth line of Table 1). Since these four averaging schemes do not average the same physical quantity, they are differently biased.

Table 1 Averaging schemes and characteristics of their biases.

There are two main causes of bias on the retrieved $X_{{\mathrm{CH}}_{\mathrm{4}}}$: the statistical bias and geophysical biases. The statistical bias, which affects every shot individually, is not produced by the averaging process and must be taken into account for shot-by-shot measurement. It is induced by the random nature of the measurement of online and offline signals into non-linear equations. Figure 3 illustrates the statistical bias on the DAOD, when online and offline signals follow normal distributions. It highlights that, in this case, the DAOD derived from these signals is no longer normally distributed, and it indicates a bias and a skewness. The second main sources of bias are called geophysical biases. These biases are induced by the process of averaging. The successive averaged shots do not sound the same portion of atmosphere (surface pressure and gas concentrations vary), they are not reflected on the same surface (reflectivity varies) and the elevation of the scattering surface is not constant in general (altitude and hard-target surface pressure vary). All these variations of geophysical quantities induce several biases on the average values.

Figure 3 Effect of the non-linearity on the DAOD distribution for a low reflectivity (0.016 ice and snow cover). Panel (a) shows that the online and offline calibrated signals are normally distributed. A significant part of online calibrated signals (orange) is negative, which makes the corresponding double shots unusable (undefined logarithm). Panel (b) shows that the DAODs corresponding to the usable calibrated signals are not normally distributed, and they present a bias and are skewed. The true DAOD is 0.53 whereas the mean of the distribution is about 0.54, which leads to a bias on the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ of approximately +34 ppb here.

The first scheme, AVX, directly averages the column mixing ratios of methane. Every shot is impacted by the statistical bias developed in Sect. 3.3.1. Furthermore, since a column with a high total molecular content and another with fewer molecules would count the same in the averaged mixing ratio, the uniform weighting of methane concentrations leads to the creation of a bias that is called geophysical bias of type 1, described in Sect. 3.4.1.

The second scheme, AVD, computes the ratio of the mean DAOD and the mean IWF. It is also impacted by the statistical bias (cf. Sect. 3.3.1). However, this scheme takes into account the fact that every column does not present the same molecular content as DAOD and IWF are averaged separately. Thus, it is not impacted by geophysical bias of type 1.

The third scheme, AVS, averages signals before computing relative transmissions, DAOD and $X_{{\mathrm{CH}}_{\mathrm{4}}}$. The statistical bias is only applicable to the resulting averaged signals such that this bias is highly reduced compared to AVX and AVD (cf. Sect. 3.3.2). However, geophysical biases are increased. First, when the DAOD varies from shot to shot – due to altitude (or surface pressure) variations or methane concentration variations for instance – the DAOD computed from averaged signals is not representative of the true mean DAOD. This is called geophysical bias of type 2, presented in Sect. 3.4.2. Secondly, for the AVS scheme, the average DAOD is weighted by the offline signal strength. Consequently, the variance of the average quantities is reduced. However, a correlation between methane concentration and reflectivity adds a bias to the retrieved quantities. This bias is called geophysical bias of type 3 and will also be discussed in Sect. 3.4.2.

The fourth scheme, AVQ, averages transmissions before computing average DAOD (and average $X_{{\mathrm{CH}}_{\mathrm{4}}})$. The transmissions for every column are averaged with a uniform weighting. Note that the major drawback of this scheme is that it mixes several bias sources that cannot be easily corrected. Indeed with the averaging being made inside the logarithm, it is not possible to separate the bias into two terms due to the measurement noise and the variation of geophysical parameters of the scene (cf. Appendix B). This scheme will not be considered in the next sections.

In the following sections, for each averaging scheme of Table 1 (except averaging of quotients), we will quantify separately the statistical bias and the geophysical biases and will in the end combine them in order to determine the total bias for various scenarios.

3.3 Statistical bias

3.3.1 Statistical bias on AVX and AVD

The averaging of columns (either $X_{{\mathrm{CH}}_{\mathrm{4}}}$ or DAOD) needs DAODs to be computed for every couple of calibrated signals (Q^off,Q^on). However, as the measurements are affected by random noise and the IPDA lidar equation (Eq. 2) is not linear, a bias appears when computing the DAOD (Fig. 3). Let us define the estimator of the DAOD $\widehat{\mathit{\delta }}$ as follows:

$$\begin{array}{}\text{(10)}& \widehat{\mathit{\delta }}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{Q^{\mathrm{off}}}{Q^{\mathrm{on}}}}\right).\end{array}$$

The total noise contributions affecting offline and online signals are statistically independent. Thus, for each single shot, the calibrated signals Q^on and Q^off can be considered as independent random variables. Furthermore, due to the relatively high number of photons in a single pulse, we can assume that these random variables are normally distributed around a mean value μ^on,off and with a standard deviation σ^on,off.

Under the normality assumption, Eq. (10) can be decomposed into three terms:

$$\begin{array}{ll}{\displaystyle }\widehat{\mathit{\delta }}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot & {\displaystyle }\ln \left({\displaystyle \frac{{\mathit{\mu }}^{\mathrm{off}}+{\mathit{\sigma }}^{\mathrm{off}}\cdot X^{\mathrm{off}}}{{\mathit{\mu }}^{\mathrm{on}}+{\mathit{\sigma }}^{\mathrm{on}}\cdot X^{\mathrm{on}}}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{{\mathit{\mu }}^{\mathrm{off}}}{{\mathit{\mu }}^{\mathrm{on}}}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }\cdot \ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{off}}}{{\mathrm{SNR}}^{\mathrm{off}}}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{on}}}{{\mathrm{SNR}}^{\mathrm{on}}}}\right),\end{array}$$

where X^on,off follows standard normal distributions. And the signal-to-noise ratios are defined as

$$\begin{array}{}\text{(12)}& {\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}={\displaystyle \frac{{\mathit{\mu }}^{\mathrm{on},\mathrm{off}}}{{\mathit{\sigma }}^{\mathrm{on},\mathrm{off}}}}.\end{array}$$

The first term of Eq. (11) is the parameter that needs to be estimated (i.e. the unbiased DAOD), and the two last terms are error terms that correspond to the bias of $\widehat{\mathit{\delta }}$ due to the non-linearity of the function:

$$\begin{array}{ll}{\displaystyle }{\mathrm{Bias}}_{\mathrm{stat}}& {\displaystyle }\left(\widehat{\mathit{\delta }}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot E\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{off}}}{{\mathrm{SNR}}^{\mathrm{off}}}}\right)\right]-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\\ \text{(13)}& {\displaystyle }& {\displaystyle }\cdot E\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{on}}}{{\mathrm{SNR}}^{\mathrm{on}}}}\right)\right].\end{array}$$

The task is now to evaluate this bias to remove or at least reduce it. Analytically, under the normal distribution hypothesis, the expected values are defined by:

$$\begin{array}{ll}{\displaystyle }E& {\displaystyle }\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{on},\mathrm{off}}}{{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}}\right)\right]={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}}}\underset{-{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}{\overset{+\mathrm{\infty }}{\int }}\\ \text{(14)}& {\displaystyle }& {\displaystyle }\ln \left(\mathrm{1}+{\displaystyle \frac{x}{{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}}\right)\cdot e^{-\frac{x^{\mathrm{2}}}{\mathrm{2}}}\mathrm{d}x.\end{array}$$

Providing that SNR^on,off is high enough, we can use a Taylor expansion of the logarithm around zero so that the bias can be approximated by the following formula (Bösenberg, 1998):

$$\begin{array}{}\text{(15)}& {\mathrm{Bias}}_{\mathrm{stat}}\left(\widehat{\mathit{\delta }}\right)\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\displaystyle \frac{\mathrm{1}}{{\left({\mathrm{SNR}}^{\mathrm{on}}\right)}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{{\left({\mathrm{SNR}}^{\mathrm{off}}\right)}^{\mathrm{2}}}}\right].\end{array}$$

The assumption that the calibrated signals follow a normal distribution does not rigorously hold when the DAOD is computed. Indeed, over dark surfaces (low reflectivity), the SNR may happen to be so low that either one or both calibrated signals Q^off and Q^on takes negative values; hence, DAOD is undefined. This can actually happen as the calibrated signal Q^on,off is computed from the raw signal that corresponds to a photon count (positive quantity) from which the estimated background energy has been subtracted. Whenever one of the calibrated signals is negative, the corresponding couple (Q^offQ^on) must be discarded. And as the lowest calibrated signals are systematically discarded from the averaging set, the average measurement is biased. This bias can be corrected by considering Eqs. (11) and (13) with X^on,off as a left-truncated normal distribution with a mean value of zero, a variance of one and a left-truncation at $-{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}$ (Johnson et al., 1994). When done so, it becomes

$$\begin{array}{ll}{\displaystyle }E& {\displaystyle }\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{on},\mathrm{off}}}{{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}}\right)\right]={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\left[\mathrm{1}-\mathrm{\Phi }\left(-{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}\right)\right]}}\\ \text{(16)}& {\displaystyle }& {\displaystyle }\underset{-{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}{\overset{+\mathrm{\infty }}{\int }}\ln \left(\mathrm{1}+{\displaystyle \frac{x}{{\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}}}\right)\cdot e^{-\frac{x^{\mathrm{2}}}{\mathrm{2}}}\mathrm{d}x,\end{array}$$

where Φ is the standard normal cumulative distribution function.

To correct the bias due to the non-linearity of the IPDA lidar equation, the SNR must be estimated. Once done, the bias correction scheme would either need to estimate the bias directly from the approximate Taylor expansion formula of Eq. (15) or estimate the bias using Eq. (13) and a numerical computation of Eq. (16). Typically, for MERLIN observations, the error made by using the Taylor expansion of Eq. (15) instead of Eq. (16) is lower than 1 ppb on the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ for a surface reflectivity value greater than 0.1 (SNR^off≈16 and SNR^on≈7) as shown on Fig. 4. Table 2 shows the error made by using the Taylor expansion instead of computing the truncated normal integral. For values of reflectivity smaller than 0.1, it would be preferable to use the exact formula for the bias presented in Eqs. (13) and (16). Further study (not presented here) shows that for very low reflectivity, the estimation of the noise induced bias is really sensitive to an error on the SNR, and this correction is no longer applicable in practice. The way statistical bias on the DAOD is translated to bias on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ will be treated in Sect. 3.5.

Figure 4 Statistical bias induced by measurement noise. The online and offline SNRs drive the value of the statistical bias. The blue line is derived from the integration of the truncated normal distribution (Eqs. 16 and 13). The orange line is the Taylor development (Eq. 15), only valid when reflectivity is high enough (i.e. high SNR). The expected bias computed from a simple Monte Carlo simulation (yellow dots) shows that the integration approach is the most accurate. For reflectivity values of 0.1 (vegetation cover), integration (blue) and Taylor development (orange), it differs by about 1 ppb (cf. Table 2 for some values).

Table 2 Error on the statistical bias estimation by using the Taylor expansion instead of using a truncated normal distribution (cf. Fig. 4).

3.3.2 Statistical bias on AVS

The third averaging scheme defined on Table 1, AVS, averages online and offline calibrated signals separately. The corresponding estimator of the average DAOD is written

$$\begin{array}{}\text{(17)}& {\widehat{\mathit{\delta }}}^{\mathrm{avs}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{\langle Q^{\mathrm{off}}\rangle }{\langle Q^{\mathrm{on}}\rangle }}\right).\end{array}$$

Consistently with Sect. 3.3.1, we consider the individual calibrated signals to be normal random variables of mean ${\mathit{\mu }}_i^{\mathrm{on},\mathrm{off}}$ and standard deviation ${\mathit{\sigma }}_i^{\mathrm{on},\mathrm{off}}$. The parameters of the distributions depend on shot i since each shot is considered as the realization of a different distribution, depending on the geophysical parameters of the scene (reflectivity, atmospheric transmission, surface pressure). The successive measurements are considered independent and, as the sum of independent normal random variables, is a normal random variable. We introduce S^on,off, the average random variable:

$$\begin{array}{}\text{(18)}& S^{\mathrm{on},\mathrm{off}}=\langle Q^{\mathrm{on},\mathrm{off}}\rangle \sim N\left(m^{\mathrm{on},\mathrm{off}},{\left({\mathit{ϵ}}^{\mathrm{on},\mathrm{off}}\right)}^{\mathrm{2}}\right),\end{array}$$

where the mean and variance of S^on,off are

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{m}^{\mathrm{on},\mathrm{off}}=\langle {\mathit{\mu }}^{\mathrm{on},\mathrm{off}}\rangle ,\text{(20)}& {\displaystyle }& {\displaystyle }{\left({\mathit{ϵ}}^{\mathrm{on},\mathrm{off}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}{\sum }_{i=\mathrm{1}}^N{\left({\mathit{\sigma }}_i^{\mathrm{on},\mathrm{off}}\right)}^{\mathrm{2}}.\end{array}$$

The empirical estimate of the SNR of the equivalent measurement S^on,off on the whole averaging window can be written

$$\begin{array}{ll}{\displaystyle }{\mathrm{SNR}}_{\mathrm{eq}}^{\mathrm{on},\mathrm{off}}& {\displaystyle }={\displaystyle \frac{m^{\mathrm{on},\mathrm{off}}}{{\mathit{ϵ}}^{\mathrm{on},\mathrm{off}}}}=\left({\sum }_{i=\mathrm{1}}^N{\mathit{\mu }}_i^{\mathrm{on},\mathrm{off}}\right)\\ \text{(21)}& {\displaystyle }& {\displaystyle }{\left({\sum }_{i=\mathrm{1}}^N{\left({\displaystyle \frac{{\mathit{\mu }}_i^{\mathrm{on},\mathrm{off}}}{{\mathrm{SNR}}_i^{\mathrm{on},\mathrm{off}}}}\right)}^{\mathrm{2}}\right)}^{-\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Given these definitions, we can write the bias due to shot random variations as in Eq. (13):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{Bias}}_{\mathrm{stat}}\left({\widehat{\mathit{\delta }}}^{\mathrm{avs}}\right)=E\left[{\widehat{\mathit{\delta }}}^{\mathrm{avs}}\right]-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{m^{\mathrm{off}}}{m^{\mathrm{on}}}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot E\\ \text{(22)}& {\displaystyle }& {\displaystyle }\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{off}}}{{\mathrm{SNR}}_{\mathrm{eq}}^{\mathrm{off}}}}\right)\right]-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot E\left[\ln \left(\mathrm{1}+{\displaystyle \frac{X^{\mathrm{on}}}{{\mathrm{SNR}}_{\mathrm{eq}}^{\mathrm{on}}}}\right)\right],\end{array}$$

Provided an estimation of the shot-by-shot SNR, we can estimate the bias of Eq. (22) with the same methods as in Sect. 3.3.1, either considering the simplified Taylor expansion approximation (Eq. 15) or the more accurate integral of truncated normal distribution (Eqs. 13 and 16). Compared to AVX and AVD schemes, the equivalent SNRs, when averaged horizontally on 150 shot pairs, are considerably larger and as a consequence, the bias is considerably smaller. The Taylor expansion approximation holds really well and an error on the estimation of SNR has a negligible impact on the bias estimation.

3.4 Geophysical bias

3.4.1 Geophysical bias of type 1 on AVX

Considering an arithmetic averaging for both AVX and AVD schemes yields different results, since the former scheme averages concentrations and the latter averages quantities that are proportional to number of molecules of methane. Whereas the AVX scheme computes the arithmetic mean of $X_{{\mathrm{CH}}_{\mathrm{4}}}$ (Eq. 23), the AVD scheme computes average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ weighted by the IWF (Eq. 24):

$$\begin{array}{}\text{(23)}& {\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avx}}=\langle {\displaystyle \frac{\mathit{\delta }}{\mathrm{IWF}}}\rangle =\langle X_{{\mathrm{CH}}_{\mathrm{4}}}\rangle \end{array}$$

and

$$\begin{array}{}\text{(24)}& {\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avd}}={\displaystyle \frac{\langle \mathit{\delta }\rangle }{\langle \mathrm{IWF}\rangle }}=\langle X_{{\mathrm{CH}}_{\mathrm{4}}}{\rangle }_{w\left[\mathrm{IWF}\right]},\end{array}$$

where

$$\begin{array}{}\text{(25)}& w_i\left[\mathrm{IWF}\right]={\displaystyle \frac{{\mathrm{IWF}}_{\mathrm{i}}}{{\sum }_{k=\mathrm{1}}^{N_{\mathrm{s}}}{\mathrm{IWF}}_k}}.\end{array}$$

For the AVX scheme, the quantity that is averaged is the column concentration which is an intensive property. If a uniform weighting is considered, there is the same contribution from columns with many molecules as from ones with less molecules. For this scheme, a variation of IWF from shot to shot (i.e. variation of altitude and/or surface pressure) leads to an overestimation of the methane content of columns that contain fewer molecules in the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$. This bias will be called geophysical bias of type 1 and is simply corrected by introducing the weighted average by the IWF. This has to be taken into account when computing the statistical bias for this scheme, as will be introduced in Sect. 3.5.

On the contrary, the AVD scheme averages the extensive properties of DAOD and IWF separately. Thus, when the DAODs are averaged, the molecule amount is preserved such that the AVD scheme is not affected by a type 1 geophysical bias.

3.4.2 Geophysical bias of type 2 and 3 on AVS

Once the bias induced by the random nature of the measurement has been subtracted, the resulting estimator is still biased by the effects of horizontal variations of geophysical quantities. Indeed, using Eq. (22), we are left with

$$\begin{array}{}\text{(26)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}={\widehat{\mathit{\delta }}}^{\mathrm{avs}}-{\mathrm{Bias}}_{\mathrm{stat}}\left({\widehat{\mathit{\delta }}}^{\mathrm{avs}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\displaystyle \frac{m^{\mathrm{off}}}{m^{\mathrm{on}}}}\right),\end{array}$$

where m^off and m^on, as defined by Eq. (19), are the average of signal expected values. Successive shot pairs are not measuring the same column of atmosphere, such that altitude, reflectivity and CH₄ concentration vary horizontally. Unlike for the AVX and AVC schemes where the ratio is computed separately, for the AVS scheme, the changing reflectivity or atmospheric transmission does not cancel out directly when computing the ratio of signals. Although measurement random noise is significantly reduced, a geophysical noise appears. We can rewrite Eq. (26) as

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{\displaystyle \frac{{\mathit{\mu }}_i^{\mathrm{off}}}{{\sum }_{k=\mathrm{1}}^N{\mathit{\mu }}_k^{\mathrm{off}}}}\cdot \exp \left(-\mathrm{2}\cdot {\mathit{\delta }}_i\right)\right)=\\ \text{(27)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{w}_i\left[{\mathit{\mu }}^{\mathrm{off}}\right]\cdot {\mathit{\tau }}_i^{\mathrm{2}}\right).\end{array}$$

Using a Taylor expansion of Eq. (27), it is possible to show that $\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}$ approximately equals the mean of the single-shot DAODs weighted by the w_i[μ^off]:

$$\begin{array}{}\text{(28)}& \stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}={\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{w}_i\left[{\mathit{\mu }}^{\mathrm{off}}\right]\cdot {\mathit{\delta }}_i+R_{\mathrm{res}},\end{array}$$

where R_res is the residual error of the linear approximations when averaging DAODs instead of transmissions. This term will be called type 2 geophysical bias. Equations (27) and (28) lead to

$$\begin{array}{ll}{\displaystyle }{R}_{\mathrm{res}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{w}_i\left[{\mathit{\mu }}^{\mathrm{off}}\right]\cdot \exp \left(-\mathrm{2}\cdot {\mathit{\delta }}_i\right)\right)\\ \text{(29)}& {\displaystyle }& {\displaystyle }-{\sum }_{i=\mathrm{1}}^N{w}_i\left[{\mathit{\mu }}^{\mathrm{off}}\right]\cdot {\mathit{\delta }}_i.\end{array}$$

Note that when the DAOD is constant all along the averaging scene, R_res is exactly zero. Furthermore, when μ^off is horizontally constant, R_res is approximately the variance of the DAOD. In fact, the term R_res is twofold: on the one hand, it is linked to DAOD fluctuations and, on the other hand, to the correlation between DAOD and reflectivity fluctuations. These correlations might occur, for instance, if there are covariations between topography (and thus DAOD) and surface type (e.g. snow on mountain tops and vegetation in valleys). In the general case, R_res is not zero and can be estimated using $\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}$ from Eq. (26), corrected for the statistical bias only, to compute a first-order estimate for $X_{{\mathrm{CH}}_{\mathrm{4}}}$:

$$\begin{array}{}\text{(30)}& X_{{\mathrm{CH}}_{\mathrm{4}}}^{\left(\mathrm{1}\right)}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}}{\langle \mathrm{IWF}{\rangle }_{w\left[Q^{\mathrm{off}}\right]}}}.\end{array}$$

Using Q^off instead of μ^off which is unknown, we can estimate R_res as

$$\begin{array}{ll}{\displaystyle }{R}_{\mathrm{res}}^{\left(\mathrm{1}\right)}& {\displaystyle }\approx -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \ln \left({\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{w}_i\left[Q^{\mathrm{off}}\right]\cdot \exp \left(-\mathrm{2}\cdot X_{{\mathrm{CH}}_{\mathrm{4}}}^{\left(\mathrm{1}\right)}\right.\right.\\ \text{(31)}& {\displaystyle }& {\displaystyle }\left.\left.\cdot {\mathrm{IWF}}_{\mathrm{i}}\right)\right)-\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{avs}}}.\end{array}$$

This process could be turned into an iterative correction. However, the first-order estimate is sufficiently accurate in all cases (not shown).

According to Eq. (28), we notice that the AVS scheme, corrected for type 2 geophysical bias, computes an average DAOD weighted by the off-signal strength. Since the main cause of variation of the offline received power is the variation of surface or hard-target reflectivity, the transmissions associated to brighter scenes count more in the average than the transmissions of darker scenes. The AVS scheme averages the measurements in such a way that a greater weight is given to high SNR signals. Consequently, this DAOD estimate is more precise (lower standard deviation) but also biased. This bias is called type 3 geophysical bias and will be defined in Sect. 3.5.

3.5 From biases on DAOD to biases on $X_{{\mathrm{CH}}_{\mathrm{4}}}$

In Sect. 3.3 and 3.4, the statistical and geophysical biases on DAOD have been derived. Here we are interested in translating biases on DAOD to biases on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ that we want to estimate. As shown by Eq. (3), $X_{{\mathrm{CH}}_{\mathrm{4}}}$ is obtained by dividing the DAOD by the IWF. This needs the IWF to be averaged horizontally, consistent with the DAOD averaging scheme. Not only is the computation of the average IWF with consistent weights important to compute $X_{{\mathrm{CH}}_{\mathrm{4}}}$, but it is also needed by the data users for the assimilation to transport models.

For the AVD scheme, the DAODs are arithmetically averaged with a uniform weight. Hence, the IWF must be averaged in the same fashion. A shot-by-shot DAOD bias according to Eq. (13) translates into a statistical bias on $X_{{\mathrm{CH}}_{\mathrm{4}}}$ as follows:

$$\begin{array}{}\text{(32)}& {\mathrm{Bias}}_{\mathrm{stat}}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avd}}\right)={\displaystyle \frac{\langle {\mathrm{Bias}}_{\mathrm{stat}}\left(\widehat{\mathit{\delta }}\right)\rangle }{\langle \mathrm{IWF}\rangle }}.\end{array}$$

For the AVX scheme, $X_{{\mathrm{CH}}_{\mathrm{4}}}$ is computed for every shot. The statistical bias on every shot is the quotient of the bias on the shot DAOD over the shot IWF. However, when horizontally averaging the statistical bias on $X_{{\mathrm{CH}}_{\mathrm{4}}}$, the type 1 geophysical bias must be taken into account (Sect. 3.4.1). The average bias should be weighted by the shot-by-shot IWF as in Eq. (24):

$$\begin{array}{}\text{(33)}& {\mathrm{Bias}}_{\mathrm{stat}}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avx}}\right)=\langle {\displaystyle \frac{{\mathrm{Bias}}_{\mathrm{stat}}\left(\widehat{\mathit{\delta }}\right)}{\mathrm{IWF}}}{\rangle }_{w\left[\mathrm{IWF}\right]}.\end{array}$$

For the AVS scheme, the IWF must be weighted consistent with the averaging scheme. Equation (28) shows that the average DAOD is weighted by the offline signal strength. As presented in the third line of Table 1, in order to keep the mixing ratio of methane consistent, the averaging of the IWF must also be weighted by w_i[Q^off]. Consistently, the translation of bias on the DAOD to bias on the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ considers the same weighting for IWF. The statistical bias translates from Eq. (22) to

$$\begin{array}{}\text{(34)}& {\mathrm{Bias}}_{\mathrm{stat}}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avs}}\right)={\displaystyle \frac{{\mathrm{Bias}}_{\mathrm{stat}}\left({\widehat{\mathit{\delta }}}^{\mathrm{avs}}\right)}{\langle \mathrm{IWF}{\rangle }_{w\left[Q^{\mathrm{off}}\right]}}}.\end{array}$$

Concerning geophysical biases, a type 2 geophysical bias (due to the linearization of DAOD variations and the correlation of signal and transmission fluctuations) described by Eq. (31) becomes

$$\begin{array}{}\text{(35)}& {\mathrm{Bias}}_{\mathrm{geo}\mathrm{2}}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avs}}\right)={\displaystyle \frac{R_{\mathrm{res}}^{\left(\mathrm{1}\right)}}{\langle \mathrm{IWF}{\rangle }_{w\left[Q^{\mathrm{off}}\right]}}}.\end{array}$$

The geophysical bias of type 3, caused by the higher sensitivity to the spots with higher reflectivity, could be written as

$$\begin{array}{}\text{(36)}& {\mathrm{Bias}}_{\mathrm{geo}\mathrm{3}}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avs}}\right)={\displaystyle \frac{\langle {\mathit{\delta }}^{\mathrm{true}}{\rangle }_{w\left[Q^{\mathrm{off}}\right]}}{\langle \mathrm{IWF}{\rangle }_{w\left[Q^{\mathrm{off}}\right]}}}-{\displaystyle \frac{\langle {\mathit{\delta }}^{\mathrm{true}}\rangle }{\langle \mathrm{IWF}\rangle }}.\end{array}$$

Indeed, the AVS scheme does not measure the true concentration of CH₄ on the 50 km window. The weighting by w_i[Q^off] implies that greater weight is granted to shots measuring brighter targets. This could be detrimental to the measurement if there were a strong correlation between reflectivity and CH₄ concentration on a global scale, which should not be the case. For assimilation or inverse modelling to models with a higher resolution than 50 km, the weighting could also be taken into account in the forward model for the $X_{{\mathrm{CH}}_{\mathrm{4}}}$.

4 Methodology to test averaging algorithms and their bias corrections

4.1 Data sets (latitude, longitude, altitude, surface pressure, and relative reflectivity)

The three averaging schemes and their associated biases will be tested on scenes modelled from real satellite data in terms of geophysical properties. For this purpose, we are interested in simulating the calibrated signals $Q_i^{\mathrm{on},\mathrm{off}}$ and the integrated weighting function IWF_i, both on a 50 km scale. To be computed, the signals require the weighting functions for every shot (WF_i,j), the volume mixing ratio of methane (${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j})$, both defined on the pressure grid (P_i,j), and the target reflectivity (ρ_i) for every shot. The integrated weighting function is computed from WF_i,j (and P_i,j). The data sets are built from satellite data provided by the SPOT-5 satellite for latitude, longitude and relative reflectivity; the Shuttle Radar Topography Mission (SRTM) digital elevation map data for topography; and the European Centre for Medium-Range Weather Forecasts (ECMWF) analyses for surface pressure, from which we deduce the pressure grid on 150 shots and 19 levels.

SPOT-5 was a CNES satellite launched in 2002 and operated until 2015 (Gleyzes et al., 2003). Amongst the five spectral bands of the High Resolution Geometric (HRG) instrument, it has a spectral band in the short-wave infrared domain (1.55 to 1.7 µm) with a spatial resolution of 20 m. This band includes the MERLIN laser wavelength and, as we expect spectral variations of surface albedo to have rather low spectral variations, we use the Spot SYSTEM SCENE level 1A product (images using radiometric corrections, equivalent radiance in W m⁻² Sr⁻¹ µm⁻¹) as a proxy of surface reflectivity. Indeed, as we were careful to select images with no clouds, we neglect the effect of atmospheric extinction on the SPOT-5 measurements. Note that we are interested here in a description of the reflectivity variations in the 50 km averaging window, not by the absolute value of reflectivity. This is why we consider this SPOT-5 product suitable, and we will scale it to any prescribed value of surface reflectivity in the simulations described hereafter anyway. The topography is taken from the SRTM digital elevation model (Jarvis et al., 2008), which has a spatial resolution of about 90 m. Surface pressure is taken from ECMWF 4D variational analyses from the long window daily archive and interpolated at SRTM grid points. A correction from difference between ECMWF Integrated Forecasting System (IFS) model topography and SRTM altitude is applied. In order to make both SRTM and SPOT-5 data consistent, the three selected SPOT-5 images are first processed by a low pass convolution to obtain a 90 m spatial resolution and then projected into the SRTM geometry. Note that the spatial resolution thus obtained is also close to MERLIN single shot footprint. Table 3 summarizes the data sets' content and resolutions.

Table 3 Data sets' resolution characteristics.

Three sites have been selected to be representative of topographic variability; they are located in the neighbourhood of three French cities: Toulouse, Millau and Chamonix. The different characteristics of the three samples are described in Table 4. Figures 5 and 6 show the variation of surface pressure and relative variations of reflectivity along the averaging scheme. Toulouse presents a medium variation of geophysical parameters (altitude and thus surface pressure), Millau presents a high variation and Chamonix a very high variation. Figure 7 shows the global cumulative distribution of standard deviations of altitude of SRTM database worldwide. We notice that 67 % and 97 % of the scenes present a lower altitude standard deviation than the one considered on the Millau and Chamonix data, respectively.

Figure 5 Surface pressure of the three scenes from the data sets. Toulouse, Millau and Chamonix present medium, high and very high variability, respectively (cf. Table 4).

Figure 6 Relative variations of reflectivity of the three scenes from the data set. Toulouse, Millau and Chamonix present medium, high and high) variability, respectively (cf. Table 4).

Table 4 Characteristics of the data used for the simulation.

For sensitivity study purposes, the reflectivity relative variations from the SPOT-5 data set are multiplied by a reference mean reflectivity that can be chosen to obtain the usable scene reflectivity. Four mean reflectivity values will be considered: 0.1 (vegetation), 0.05 (mixed water and vegetation), 0.025 (sea and ocean) and 0.016 (ice and snow).

Figure 7 Global cumulative distribution of the standard deviation of altitude obtained on SRTM. A total of 46 %, 67 % and 97 % of SRTM boxes present a lower standard deviation than the Toulouse scene, Millau scene and Chamonix scene, respectively. The three scenes are representative of medium, high and very high variations of altitude.

The pressure grid P_i,j and the pressure thickness grid ΔP_i,j are obtained from surface pressure $P_i^{\mathrm{surf}}$ from the ECMWF analyses data set using a hybrid sigma coordinate system.

The methane volume mixing ratio, ${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j}$, is arbitrarily set to values assumed to be realistic. For every shot i and layer j

$$\begin{array}{}\text{(37)}& \left\{\begin{array}{cc}{\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j}=\mathrm{1780}\mathrm{ppb}& \mathrm{if}{P}_{i,j}<(max\left(p_i^{\mathrm{surf}}\right)-min\left(p_i^{\mathrm{surf}}\right))/\mathrm{2}\\ {\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j}=\mathrm{1880}\mathrm{ppb}& \mathrm{otherwise}\end{array}\right..\end{array}$$

This replicates the possible correlation between methane concentration and altitude (more methane in valleys and less over mountain tops).

Finally, the weighting functions are calculated, as described in Eq. (4), from methane absorption cross sections and meteorological data (ΔP_i,j, temperature, humidity). They are computed using CH₄ absorption cross sections from the 4A radiative transfer model (Scott and Chédin, 1981; Chéruy et al., 1995) on a reference winter mid-latitude atmosphere from the Thermodynamic Initial Guess Retrieval (TIGR) data set (Chevallier et al., 2000). The sensitivity to the thermodynamic condition of the atmosphere has been tested and is negligible here (not shown).

4.2 Overall test framework

The aim of the simulation is to compare the biases of the estimated $X_{{\mathrm{CH}}_{\mathrm{4}}}$ for several averaging schemes and to evaluate the accuracy of the bias correction. A global description of the simulation is presented on Fig. 8. Each simulation case considers a typical number of N_s=150 double shots per averaging window, approximately corresponding to 50 km along the satellite ground track. It relies on a description of the geophysical scene in terms of surface pressure $P_i^{\mathrm{surf}}$, reflectivity ρ_i, an arbitrary CH₄ concentration field ${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i,j}$ and weighting functions WF_i,j (cf. Sect. 4.1). Then, the online and offline calibrated signals are computed from surface reflectivity and a random noise simulation, and the weighting functions are integrated (cf. Sect. 4.3). Next, we proceed to the computation of the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ on 50 km resolution with the different averaging schemes (AVQ not simulated) and the correction algorithms, presented in Sect. 4.4 and Table 5.

Figure 8 Global description of the simulation. Data sets (blue) are described in Sect. 4.1. Signals and the IWF computation (orange) are described in Sect. 4.3. Averaging strategies performed and their related bias corrections (green) are described in Sect. 4.4 and Table 4. Target $X_{{\mathrm{CH}}_{\mathrm{4}}}$ computation (red) is described in Sect. 3.1. ΔX is the scheme bias which is the difference between scheme and target $X_{{\mathrm{CH}}_{\mathrm{4}}}$, and Bias_res is the residual bias when evaluated biases have been subtracted from the scheme bias.

Table 5 Computational details about averaging schemes and bias evaluation.

In order to estimate the bias, the computation of an average column-integrated methane concentration $X_{{\mathrm{CH}}_{\mathrm{4}}}$ from the shot-by-shot volume mixing ratio profiles, ${\mathrm{vmr}}_{{\mathrm{CH}}_{\mathrm{4}},i}\left(p\right)$, is needed and will be computed as ${\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{T}}$ in Eq. (9).

In order to assess the performance of averaging schemes and bias correction algorithms, the standard deviation and mean of the difference ΔX between the ${\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{scheme}}$ estimated from one of the studied averaging schemes and the target value ${\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{T}}$ must be computed over a set of M simulations. The number of simulation M has to be high enough to compute the residual bias (empirical mean of ΔX) with sufficient accuracy. Let us denote σ the standard deviation of the distribution of the variable ΔX, S_M=〈ΔX〉 the empirical mean over M samples (i.e. the empirical estimate of the bias of the averaging scheme). To get an estimate of the expected value of ΔX with an accuracy of 0.1 ppb with 90 % confidence, it requires approximately M=300 000 samples according to the central limit theorem.

The typical standard deviation can also be evaluated from the sample and is approximately 22 ppb for the typical case (mean reflectivity of 0.1).

4.3 Simulation of online and offline lidar calibrated signals and IWF

Once the scene parameters are defined on the 50 km averaging window and the atmosphere is modelled, the online and offline calibrated signals must be simulated. We first have to compute the deterministic values of the calibrated signals without noise and simulate the random noise that affects them. The values of the signals are determined by the scene reflectivity (for both online and offline signals) and by the atmospheric transmission (online signals only). From the weighting functions, the methane field and the pressure field, we compute the reference DAOD, denoted ${\mathit{\delta }}_i^{\mathrm{true}}$, as the numerator of Eq. (5). Then we compare the transmission for each double shot according to

$$\begin{array}{}\text{(38)}& {\left({\mathit{\tau }}_i^{\mathrm{true}}\right)}^{\mathrm{2}}=\exp \left(-\mathrm{2}\cdot {\mathit{\delta }}_i^{\mathrm{true}}\right).\end{array}$$

From them, considering the reflectivity, we are able to determine the relative value of the online and offline mean calibrated signals:

$$\begin{array}{}\text{(39)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_i^{\mathrm{off}}={\mathit{\rho }}_i,\text{(40)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_i^{\mathrm{on}}={\mathit{\rho }}_i\cdot {\left({\mathit{\tau }}_i^{\mathrm{true}}\right)}^{\mathrm{2}},\end{array}$$

where i is the shot index, τ_i is the transmission and ρ_i is the reflectivity. Note that any constant affecting both online and offline signals can be disregarded here.

Then, Gaussian random noise has to be added to the values of the signals. It is computed from the SNR that depends on the number of photons reaching the detector (i.e. ${\mathit{\mu }}_i^{\mathrm{on},\mathrm{off}})$. Figure 9 shows the theoretical dependence of the SNR to the reflectivity according to instrument characteristics and Appendix A presents the noise model used.

Figure 9 Online and offline SNR computed from reflectivity according to instrument characteristics.

The IWF_i are simply computed by integrating the WF_i,j on all pressure layers as the denominator of Eq. (5).

4.4 Tested averaging algorithms and bias corrections

The simulation tested the three averaging schemes described in Sect. 3.2: AVX (Table 1, line 1), AVD (Table 1, line 2) and AVS (Table 1, line 3). Table 5 details the computational steps used for averaging, statistical bias evaluation and geophysical bias evaluation for the three schemes. For the AVX and AVD schemes, as explained in Sect. 3.3.1, signal couples with at least one negative calibrated signal must be discarded to compute the shot DAOD. However, since signals are averaged first for the AVS scheme, the probability that one of the averaged signals is negative is extremely small. Thus, no negative calibrated signal discarding is needed for the AVS scheme.

Table 6 Resulting bias (in ppb) for the AVD scheme after noise induced bias correction.

Table 7 Resulting bias (in ppb) for the AVS scheme after noise induced bias correction and geophysical induced bias correction.

Concerning statistical bias evaluation, an SNR estimation is needed. It is directly estimated from instrument parameters and online and offline calibrated signal strength. Once the SNR is estimated, as described in Sect. 3.3, there are two options to evaluate the statistical bias either using the Taylor expansion approximation or the numerical integral of a truncated normal distribution. Contrary to AVS, where Taylor expansion and the numerical integral make a negligible difference, for AVX and AVD, it is better to use the numerical integral as it is more accurate, and this is what is done here.

Type 1 geophysical bias, that affects the AVX scheme, is already compensated by weighting the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ and the average statistical bias by the IWF. The AVD scheme is not affected by geophysical biases. However, type 2 and type 3 geophysical biases affect the AVS scheme. Type 2 geophysical bias is evaluated by Eq. (35) using the first-order $X_{{\mathrm{CH}}_{\mathrm{4}}}$ of Eq. (30). The type 3 geophysical bias is not evaluated and will not be corrected because a correlation between reflectivity and CH₄ concentration is unlikely to occur. Indeed, on a smaller than 50 km scale the typical atmospheric transport should smear out the CH₄ concentration very effectively over areas larger than the small-scale reflectivity jumps even if this is not true for narrow valleys.

5 Results

5.1 Comparison of averaging schemes

5.1.1 Bias of averaging schemes without bias correction

The first results presented here are the respective biases of each averaging scheme without any bias correction. Figure 10 shows the bias on the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ on the three scenes (Toulouse, Millau and Chamonix) for the three averaging schemes that have been studied (AVS, AVD and AVX) without any correction. The bias due to measurement noise and due to geophysical variation appears on the results, as it is not yet subtracted. For the AVS scheme we compare the results with the uniformly weighted average IWF and the average of the IWF weighted by the offline calibrated signal strength: w_i[Q^off]. For the AVD scheme, a uniform weight is considered. And, for the AVX scheme, we compare the uniformly weighted average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ (Eq. 23) and the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ weighted by the IWF: w_i[IWF] (Eq. 24). The mean reflectivity is set to the typical value of 0.1.

Figure 10 Bias before correction for the three studied averaging schemes (red dotted lines: targeted bias ±1 ppb).

For the AVS scheme on Toulouse and Millau scenes, where there are medium to high variations of geophysical quantities, the bias is contained in the ±1 ppb range. However, it is higher on the Chamonix scene, where there are very high variations of geophysical parameters. As expected, the bias for the AVS scheme is mainly impacted by the variations of the geophysical parameters over the scene. Note that on the Chamonix scene, the weighting of the average IWF by the offline calibrated signal strength reduces this bias.

On the contrary, the bias of the AVD and AVX schemes is not affected by the geophysical variations but is mainly driven by the measurement noise, which essentially depends on the scene's mean reflectivity. As shown in Sect. 3.4.1, the AVD scheme with uniform weighting and the AVX scheme weighted by the integrated weighting function (w_i[IWF] weights) show the same bias. Although the comparison between the AVX scheme weighted uniformly (light red on Fig. 10) and the AVD scheme (green on Fig. 10) shows that their biases are close when variations of surface pressure (main cause of variations of IWF) are low (Toulouse, Millau), they become significant when variations are higher (Chamonix).

Without any correction and for the typical reflectivity, the AVS scheme is less biased than the AVD and AVX schemes. However, as we have seen in previous sections, there are ways to estimate the biases and to correct them. The following section will show the results after estimation and correction of the bias induced by the measurement noise.

5.1.2 After correction of statistical bias

As explained in Sect. 3.4, the random nature of the measurement associated with the non-linearity of the measurement equation implies that the estimation of the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ is biased. The statistical bias corrections for AVS, AVD and AVX are based on an estimation of the online and offline SNRs for the measured calibrated signals (cf. Sect. 4.4 and Table 5). Figure 11 shows the residual biases, after subtraction of estimated statistical bias, for the three averaging schemes, with and without relevant weightings, on the three studied scenes and for the typical reflectivity of 0.1. The chosen estimation of the bias is done by numerically computing the integral of the truncated Gaussian distribution (Sect. 3.3).

Figure 11 Residual bias after statistical bias correction for the three studied averaging schemes (red dotted lines: targeted bias ±1 ppb).

We see that the biases of the AVD and AVX schemes are significantly reduced (absolute value decrease by 85 % to 90 %) on every scene. The residual bias is caused by the fact that the SNRs are estimated from the noisy calibrated signals so that the estimation of the bias is not perfectly accurate. This implies that the calibrated signal outcomes from the lower part of the distribution lead to a high error on the estimated bias. This effect could be slightly compensated if, instead of discarding all the negative or null calibrated signals (extremely rare for a reflectivity value of 0.1 over 150), we discarded calibrated signals higher than a strictly positive threshold (e.g. 0.01, not shown). This would lead to a better correction and thus a lower bias, but at the cost of discarding more single-shot observations.

Table 8 Bias (in ppb) of several averaging biases before any bias correction schemes are applied on the three scenes for four reflectivity values.

For the AVS scheme, as the signals are averaged first, the equivalent SNR is very high (${\mathrm{SNR}}_{\mathrm{eq}}^{\mathrm{off}}\approx \mathrm{190}$ and ${\mathrm{SNR}}_{\mathrm{eq}}^{\mathrm{on}}\approx \mathrm{90})$ on the scenes with a mean reflectivity of 0.1. Consequently, the bias due to the equivalent measurement noise is really low (about 0.1 ppb), and this bias correction has only a small effect on the residual bias.

Taking into account the correction of the bias induced by the measurement noise, the AVS scheme still presents a lower bias on Toulouse and Millau scenes than the bias of the AVX and AVD schemes. However, on the Chamonix scene, where the geophysical variations are very high, the AVX and AVD schemes are less biased than AVS.

5.1.3 After correction of geophysical biases

The biases induced by the variation of the geophysical parameters (cf. Sect. 3.4) does not affect the AVD scheme, as the additive properties of DAOD and IWF are averaged separately. The variation of the IWF affects the bias of the AVX scheme and has already been corrected by introducing the w_i[IWF] weights when directly averaging mixing ratios. The AVS scheme is the one most affected by the variations of the geophysical variations, as seen in Sect. 5.1.1.

Figure 12 Residual bias after noise induced bias and geophysical variation induced bias corrections for the three studied averaging schemes (red dotted lines: targeted bias ±1 ppb).

Figure 12 shows the residual bias after the corrections of the statistical bias induced by the measurement noise and the variations of geophysical parameters (cf. Sect. 4.4 and Table 5). We notice that the residual bias for the AVS scheme is considerably reduced when the average weighting function is weighted by the offline calibrated signal strength. Furthermore, the iterative estimation of the bias converges at the first iteration of Eqs. (26) to (31).

Once geophysical biases are subtracted, the three scenes present a low bias. The mean residual bias on the three scenes for the AVD and AVX schemes is approximately $-\mathrm{2.1}\pm \mathrm{0.1}$ ppb, whereas for the AVS scheme, it is approximately $-\mathrm{0.09}\pm \mathrm{0.09}$ ppb. After all corrections, even on highly structured scenes, AVS is the least biased scheme of the three studied schemes. When the average IWF is consistently weighted with the w_i[Q^off] weights, the geophysical induced bias is almost completely removed.

5.2 Impact of the mean reflectivity on the residual bias

All results presented above are computed for scenes with a mean reflectivity of 0.1, which roughly corresponds to vegetation cover. For the purpose of choosing the least biased algorithm to compute average $X_{{\mathrm{CH}}_{\mathrm{4}}}$, it is interesting to test the robustness to reflectivity. Indeed, reflectivity is the main driver for the expected value of SNR; low reflectivity scenes lead to lower SNR and consequently higher bias. Tables 6 and 7 show the residual bias comparing four different reflectivity values: 0.1 (vegetation), 0.05 (mixed sea and vegetation), 0.025 (sea and ocean) and 0.016 (ice and snow). Table 6 gives the residual bias for the AVD scheme (the AVX scheme presents similar results), and Table 7 shows the residual bias for the AVS scheme, where the average IWF is weighted by the offline calibrated signal strength (w_i[Q^off] weights) and corrected of the bias due to geophysical variations from shot to shot. For both tables, the Taylor expansion bias correction (Eq. 15) and numerical computation of the expectation (so-called integral truncated normal distribution Eqs. 16 and 13) are compared.

First, as seen in Table 6 (AVD scheme), the Taylor bias correction does not succeed in quantifying the bias on any of the four mean reflectivity values. The uncertainties are too high and prevent quantitative analysis of the results. This is due to the fact that there are some calibrated signals that are really close to zero and for which the SNR is underestimated; thus the bias (and standard deviation) is overestimated. This could be mitigated by the choice of a higher threshold of the usable calibrated signal before the computation of the DAOD (not shown). The results when using the integral bias correction on AVD are more physical. However, they also show an over estimation of the bias, especially for low reflectivity values. In every case for the AVD scheme, the bias threshold is exceeded.

Table 7 gives the results of the robustness of the AVS scheme to decreasing reflectivity. Unlike the AVD scheme, the AVS scheme, when all corrections are made, presents satisfying results for all reflectivity values, and in every scene the biases remain contained into the threshold interval of ±1 ppb. The effect of the decreasing reflectivity has a very small impact on the residual bias.

To summarize, the best algorithm to limit the bias for MERLIN processing algorithms is clearly the AVS scheme, with an average IWF weighted by the offline calibrated signal strength and both corrections of the geophysical bias and the bias induced by the measurement noise (either Taylor or integral bias correction). On every scene and for all expected reflectivity values, this algorithm is compliant with the averaging bias specifications of the MERLIN mission. Note that this conclusion holds in the case where all the 150 shots are considered; in the case of a partially cloudy window where only a subsample of clear sky shots are averaged, the AVS will still be the best averaging scheme, but the performance will be decreased.

6 Conclusions

The French–German space-borne IPDA lidar mission MERLIN will measure the average integrated column dry-air mixing ratio of methane ($X_{{\mathrm{CH}}_{\mathrm{4}}})$ on a 50 km scale. The processing algorithms must limit both the relative random error (RRE) and the relative systematic error (RSE) on the $X_{{\mathrm{CH}}_{\mathrm{4}}}$. As the IPDA technique relates the signal measurements to the $X_{{\mathrm{CH}}_{\mathrm{4}}}$ by a non-linear equation, a simple and naive averaging can lead to high biases.

Three averaging schemes have been studied: averaging of $X_{{\mathrm{CH}}_{\mathrm{4}}}$ (AVX), averaging of DAOD (AVD) and averaging of signals (AVS). For these averaging schemes, possible sources of bias can either be due to the measurement noise, the variation of the geophysical parameters on the averaging scene or both.

The three schemes are sensitive to the bias induced by the measurement noise even if AVS is far less impacted for the typical reflectivity. This bias can be corrected by a formula introducing the estimated SNR on the measured signals if the SNR is high enough. The bias due to the variation of geophysical parameters does not affect the AVD scheme because it directly averages the desired additive quantities. On the contrary, the AVX scheme must average the concentration weighted by the integrated weighting function (IWF) in order to average a molecule number instead of averaging concentrations. The third scheme AVS measures the average $X_{{\mathrm{CH}}_{\mathrm{4}}}$ weighted by the offline signal strength, which means that more weight to the measurements with a high SNR is given when averaging. The bias of this scheme is sensitive to the variation of geophysical parameters (surface pressure and surface reflectivity). This bias is corrected using an iterative process with the uncorrected $X_{{\mathrm{CH}}_{\mathrm{4}}}$ as first guess.

These averaging schemes and their bias corrections have been tested on scenes modelled from real satellite data in terms of altitude, surface pressure, weighting functions and relative variations of reflectivity. The three scenes present interesting characteristics, as they show different geophysical variations that could impact averaging biases. Besides, the signals and random noise are simulated from geophysical parameters and instrument parameters.

The simulation shows that the lowest biases are obtained for the AVS scheme using appropriate bias corrections and averaging weights. Furthermore, this scheme is robust to low reflectivity values unlike the AVX and AVD schemes, which are highly sensitive to the accuracy of the SNR estimation. The best scheme, AVS, is compliant with the allocated averaging bias requirements (RSE) of 0.06 % (1 ppb for a $X_{{\mathrm{CH}}_{\mathrm{4}}}$ of 1780 ppb) for the whole range of expected reflectivity values (from 0.1 down to 0.016).

A continuation of this study could evaluate the sensitivity of a poor (unprecise or biased) estimation of the SNR on the estimation of the bias due to measurement noise for low reflectivity values. Furthermore, the use of the lidar simulator and processor suites, currently in development at the LMD, could be beneficial to the evaluation of the biases, and more specifically of the averaging biases, on a wider scale (many scenes, atmosphere types, etc.).

Appendix A: Signal generation and noise model

The simulation of calibrated signals requires a noise model. The signal distribution is considered to be Gaussian first, because the number of photons that reaches the photodetector is high enough for the Poisson distribution to be considered as Gaussian. Secondly, the system is limited by the detector noise that is mainly thermal noise, which is normally distributed.

The calibrated signals are produced using a pseudorandom number generator. The expected values of the calibrated signal distributions, μ^on,off, depend on the atmospheric transmission and the surface reflectivity as presented in Eqs. (39) and (40). Then the standard deviation σ^on,off is deduced from the SNR, which is defined as

$$\begin{array}{}\text{(A1)}& {\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}={\displaystyle \frac{{\mathit{\mu }}^{\mathrm{on},\mathrm{off}}}{{\mathit{\sigma }}^{\mathrm{on},\mathrm{off}}}}.\end{array}$$

And the SNR model is described by

$$\begin{array}{}\text{(A2)}& {\mathrm{SNR}}^{\mathrm{on},\mathrm{off}}={\left({\displaystyle \frac{N^{\mathrm{2}}}{a+b\cdot N+c\cdot N^{\mathrm{2}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where N is the number of photoelectrons, and a, b and c are parameters computed from the MERLIN system parameters. This is illustrated on Fig. 9.

The number of photons is computed from the reflectivity and the atmospheric transmission. In the standard case (reflectivity of 0.1 sr⁻¹), its values are approximately 3000 for the offline pulse and 1000 for the online pulse.

The first term of the denominator corresponds to the detector noise, the second to the shot noise and the third to the speckle. Note that the speckle term has been neglected in this article, whereas both detector noise and shot noise have been considered, as they are dominant compared to speckle noise.

Appendix B: Averaging of quotients

The averaging of quotients estimates the average of the shot-by-shot two-way transmissions ${\mathit{\tau }}_i^{\mathrm{2}}$. Due to the measurement noise, we will suppose that, for the shot i, the online and offline measured calibrated signals, $Q_i^{\mathrm{on}}$ and $Q_i^{\mathrm{off}}$, respectively, are outcomes of normal distributions with mean values denoted ${\mathit{\mu }}_i^{\mathrm{on}}$ and ${\mathit{\mu }}_i^{\mathrm{off}}$, respectively, and standard deviations denoted ${\mathit{\sigma }}_i^{\mathrm{on}}$ and ${\mathit{\sigma }}_i^{\mathrm{off}}$, respectively. The $X_{{\mathrm{CH}}_{\mathrm{4}}}$ computed from averaging quotients can be defined as

$$\begin{array}{}\text{(B1)}& {\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avq}}={\displaystyle \frac{-\frac{\mathrm{1}}{\mathrm{2}}\cdot \ln \langle {\mathit{\tau }}^{\mathrm{2}}\rangle }{\langle \mathrm{IWF}\rangle }},\end{array}$$

with the transmission defined as

$$\begin{array}{}\text{(B2)}& \langle {\mathit{\tau }}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\displaystyle \frac{Q_i^{\mathrm{on}}}{Q_i^{\mathrm{off}}}}.\end{array}$$

If we define the standardized signals corresponding to $Q_i^{\mathrm{on}}$ and $Q_i^{\mathrm{off}}$ as $X_i^{\mathrm{on}}$ and $X_i^{\mathrm{off}}$, the average transmission can be written as

$$\begin{array}{}\text{(B3)}& \langle {\mathit{\tau }}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\displaystyle \frac{{\mathit{\mu }}_i^{\mathrm{on}}+{\mathit{\sigma }}_i^{\mathrm{on}}\cdot X_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{off}}+{\mathit{\sigma }}_i^{\mathrm{off}}\cdot X_i^{\mathrm{off}}}}.\end{array}$$

Then we can further separate the random part due to measurement noise and the deterministic part due to varying geophysical parameters as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\langle {\mathit{\tau }}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left[{\displaystyle \frac{{\mathit{\mu }}_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{off}}}}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{on}}}}{X}_i^{\mathrm{on}}\right)\right.\\ \text{(B4)}& {\displaystyle }& {\displaystyle }\left.{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{off}}}{{\mathit{\mu }}_i^{\mathrm{off}}}}{X}_i^{\mathrm{off}}\right)}^{-\mathrm{1}}\right],{\displaystyle }& {\displaystyle }\langle {\mathit{\tau }}^{\mathrm{2}}\rangle =e^{-\mathrm{2}\langle {\mathit{\delta }}^{\mathrm{true}}\rangle }\cdot {\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left[e^{-\mathrm{2}\mathrm{\Delta }{\mathit{\delta }}_i}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{on}}}}{X}_i^{\mathrm{on}}\right)\right.\\ \text{(B5)}& {\displaystyle }& {\displaystyle }\left.{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{off}}}{{\mathit{\mu }}_i^{\mathrm{off}}}}{X}_i^{\mathrm{off}}\right)}^{-\mathrm{1}}\right],\end{array}$$

where 〈δ^true〉 is the average DAOD computed from noiseless mean signals, and Δδ_i the difference compared to the shot-by-shot DAOD. Then we can deduce the error of AVQ scheme as follows:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avq}}& {\displaystyle }=X_{{\mathrm{CH}}_{\mathrm{4}}}^{\mathrm{true}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\langle \mathrm{IWF}\rangle }}\cdot \ln \left({\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left[e^{-\mathrm{2}\mathrm{\Delta }{\mathit{\delta }}_i}\right.\right.\\ \text{(B6)}& {\displaystyle }& {\displaystyle }\left.\left.\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{on}}}}{X}_i^{\mathrm{on}}\right){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{off}}}{{\mathit{\mu }}_i^{\mathrm{off}}}}{X}_i^{\mathrm{off}}\right)}^{-\mathrm{1}}\right]\right).\end{array}$$

The corresponding bias is the expected value of the error term:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{Bias}\left({\stackrel{\mathrm{‾}}{X_{{\mathrm{CH}}_{\mathrm{4}}}}}^{\mathrm{avq}}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\langle \mathrm{IWF}\rangle }}\cdot E\left[\ln \left({\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\right.\right.\\ \text{(B7)}& {\displaystyle }& {\displaystyle }\left.\left.\left[e^{-\mathrm{2}\mathrm{\Delta }{\mathit{\delta }}_i}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{on}}}{{\mathit{\mu }}_i^{\mathrm{on}}}}{X}_i^{\mathrm{on}}\right){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_i^{\mathrm{off}}}{{\mathit{\mu }}_i^{\mathrm{off}}}}{X}_i^{\mathrm{off}}\right)}^{-\mathrm{1}}\right]\right)\right].\end{array}$$

In Eq. (B5), the empirical average transmission is decomposed into two factors. The first is the transmission corresponding to the average DAOD from noiseless signals. The second factor is the average of multiplicative errors that is the deterministic error from geophysical variations on the averaging scene and the random factors due to the presence of measurement noise. As shown in Eq. (B7), the error sources are mixed into the non-linear function, which makes them difficult to evaluate. It is possible to derive a suitable bias correction based on Eq. (B7) for AVQ, but in the end it is not expected to be better than the other ones. Table 8 presents the biases for every averaging scheme presented before any bias correction is applied.