Remote sensing of methane plumes: instrument tradeoff analysis for detecting and quantifying local sources at global scale

Methane (CH₄) is the second most important anthropogenic greenhouse gas with a significant impact on radiative forcing, tropospheric air quality, and stratospheric water vapor. Remote sensing observations enable the detection and quantification of local methane emissions across large geographical areas, which is a critical step for understanding local flux distributions and subsequently prioritizing mitigation strategies. Obtaining methane column concentration measurements with low noise and minimal surface interference has direct consequences for accurately determining the location and emission rates of methane sources. The quality of retrieved column enhancements depends on the choices of the instrument and retrieval parameters. Here, we studied the changes in precision error and bias as a result of different spectral resolutions, instrument optical performance, and detector exposure times by using a realistic instrument noise model. In addition, we formally analyzed the impact of spectrally complex surface albedo features on retrievals using the iterative maximum a posteriori differential optical absorption spectroscopy (IMAP-DOAS) algorithm. We built an end-to-end modeling framework that can simulate observed radiances from reflected solar irradiance through a simulated CH₄ plume over several natural and artificial surfaces. Our analysis shows that complex surface features can alias into retrieved methane abundances, explaining the existence of retrieval biases in current airborne methane observations. The impact can be mitigated with higher spectral resolution and a larger polynomial degree to approximate surface albedo variations. Using a spectral resolution of 1.5 nm, an exposure time of 20 ms, and a polynomial degree of 25, a retrieval precision error below 0.007 mole m⁻² or 1.0 % of total atmospheric CH₄ column can be achieved for high albedo cases, while minimizing the bias due to surface interference such that the noise is uncorrelated among various surfaces. At coarser spectral resolutions, it becomes increasingly harder to separate complex surface albedo features from atmospheric absorption features. Our modeling framework provides the basis for assessing tradeoffs for future remote sensing instruments and algorithmic designs. For instance, we find that improving the spectral resolution beyond 0.2 nm would actually decrease the retrieval precision, as detector readout noise will play an increasing role. Our work contributes towards building an enhanced monitoring system that can measure CH₄ concentration fields to determine methane sources accurately and efficiently at scale.

1 Introduction

Anthropogenic greenhouse gas emissions have been rising continuously, affecting the global climate and the environment (Stocker et al., 2013). Among the most important anthropogenic emissions are carbon dioxide (CO₂) and methane (CH₄). Due to a much shorter lifetime of CH₄ (≈9 years) compared to CO₂ (≈500 years), CH₄ has gained attention as a target for mitigation efforts to achieve short- and medium-term reductions in global warming (Montzka et al., 2011; Prather et al., 2012; Shindell et al., 2012). In general, anthropogenic methane emissions are also much more uncertain than those of carbon dioxide, which can often be characterized to within approximately 10 % just from budget assumptions (Gurney et al., 2019). For instance, just the question of whether or not the leak rate in the natural gas extraction system is 1 % or 2 % is equivalent to a 100 % uncertainty in methane emissions from these leaks. At the same time, leak rate outliers (Frankenberg et al., 2016; Duren et al., 2019; Cusworth et al., 2021) are often local in nature and easily fixable, representing a win–win scenario if faulty equipment or practices can be readily detected and then efficiently mitigated. As CH₄ reduction plays a significant role in climate mitigation efforts, one key step in emission reduction is determining where these emissions are coming from. This is underpinned in the 2018 NASA Decadal Survey, which names the identification and understanding of CH₄ emissions as one of the top priorities in the efforts to improve future climate projections and help lead the way in emission reduction (National Academies of Sciences and Medicine, 2018).

Remote sensing instruments using absorption spectroscopy have been available as one effective solution for measuring atmospheric CH₄ concentration over large geographical areas. Space-based CH₄ retrieval techniques from satellite observations such as the SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY (SCIAMACHY; Frankenberg et al., 2005a, 2011) and the Greenhouse gases Observing SATellite (GOSAT; Parker et al., 2011, 2015; Turner et al., 2015) were dedicated missions with CH₄ as a key target. They used CH₄ absorption features in the 1.6 or 2.3 µm bands to retrieve column CH₄ concentration across the globe. The TROPOspheric Monitoring Instrument (TROPOMI), with a spatial resolution of a few kilometers, has also been shown to be capable of identifying regions of high emissions (de Gouw et al., 2020; Hu et al., 2018). These satellites, which have been designed by the atmospheric community, have particular sets of goals and instrument specifications that are mostly targeted towards obtaining observations of the regional-scale methane distributions with high accuracy and precision. Most of these satellites were designed to measure gradients of methane concentration across hundreds to thousands of kilometers of scale as this enables surface flux inversions at the global scale. Typically, all of these instruments have one feature in common – they have very high spectral resolution (0.05–0.25 nm) to distinguish individual methane absorption lines from spectrally smooth surface albedo variations. However, due to their coarse spatial resolutions, the measurements are not yet at a level where local sources can be identified, attributed to a specific source type (e.g., compressor station or well pad), and mitigated directly.

One potential solution to fill this scale gap is to use an airborne instrument that has a much higher spatial resolution, such as the Methane Airborne MAPper (MAMAP; Gerilowski et al., 2011) or the next-generation Airborne Visible InfraRed Imaging Spectrometer (AVIRIS-NG; Thorpe et al., 2017). The latter is based on the insight that methane column enhancements at high spatial resolution (a few meters) can be so high that the retrieval of the absorbing feature can be done, even with a moderate spectral resolution (5–10 nm). If the methane column is expressed similar to the Dobson unit, i.e., as the thickness of a layer of pure gas which would be formed by the total column amount at standard conditions, then the layer thickness at current background methane conditions would only be about 1.6 cm. Thus, a pure methane layer of only 1.6 mm would enhance the total column by 10 %, which is certainly realistic for measurements of methane point sources at fine spatial resolution.

Bradley et al. (2011) and Thorpe et al. (2014) were among the first to show that moderate resolution instruments can detect methane plumes, even when the strong 2.3 µm methane band is convolved with the AVIRIS (resolution of 10 nm) or AVIRIS-NG (5 nm) instrument line shape functions. While individual lines are hard to resolve, the strong methane band in this range causes enough fine structure in terms of bulk absorptions by a multitude of methane lines within the instrument resolution. Previous studies by Frankenberg et al. (2016), Duren et al. (2019), and Cusworth et al. (2021) have utilized AVIRIS-NG to conduct field campaigns in California and the Four Corners region, where they could create a map of methane enhancements in the area and detect several hundreds of individual methane sources, which followed a heavy tail flux distribution. The concept of this airborne spectrometer provides a promising opportunity for local source detection and quantification. However, the instrument was not originally designed for methane detection, and it does not meet the same precision and accuracy requirements as those satellites from the atmospheric community for methane retrieval at a global scale, which require precision better than 1 %, which is equivalent to enhancements of about 19 ppb (parts per billion) in XCH₄ (or 4×10¹⁷ $\mathrm{molec}..{\mathrm{cm}}^{-\mathrm{2}}$, 0.007 mole m⁻² or 152 ppm m). One significant drawback of coarse spectral resolution is the occurrence of retrieval artifacts that often correlate with specific surface features (see Fig. 1). This can confound the detection and quantification of methane point sources in the analysis and obviate the robust detection of subtle gradients at larger spatial scales (Jongaramrungruang et al., 2021). Even though most strong plumes can be observed, the uncertainties in the overall detection and quantification at the regional level can present persistent problems and often involve human judgment to isolate plumes from artifacts. In fact, during the California survey, the Four Corners study, and the Permian survey (Duren et al., 2019; Frankenberg et al., 2016; Cusworth et al., 2021), human analysts were involved in a manual process to look through each flight line to classify true emission sources from false positives. Similarly, in previous space-based studies to locate and approximate a large emission source, such as a blowout event, prior information about the location of the source is usually already known, making it much easier to find a true methane source from space-based measurements over the area after the fact. There are ongoing efforts to develop automated plume detection for existing instruments. Future real-time monitoring systems would greatly benefit from next-generation instruments that would reduce retrieval artifacts and provide retrievals with improved accuracy, such that remote sensing measurements can be analyzed to locate and quantify plumes automatically and at scale. This motivates the origination of this study.

Figure 1 Example of the systematic outliers from a retrieved AVIRIS-NG scene (b) compared with an RGB image (a).

If we had the opportunity to design a new instrument that is optimized for methane retrievals at fine spatial resolution (sub 50 m), then what would the specifications of this instrument look like? Thorpe et al. (2016) proposed a 1 nm instrument to mitigate the drawbacks of AVIRIS-NG. To fully evaluate optimal performance metrics, we have to consider the tradeoff between spectral and spatial resolutions and concomitant changes in detector noise characteristics. On the one hand, the instrument needs to meet the requirements of the atmospheric community so that it can unambiguously differentiate methane from other confounding factors. On the other hand, the instrument should have an adequate integration time to achieve a high spatial resolution with sufficient signal-to-noise levels. Here, we investigate this tradeoff and evaluate the risks and benefits for methane retrievals at fine resolutions with the purpose of successfully detecting and quantifying local sources in mind. We built an end-to-end modeling framework that can generate reflected solar radiance through a methane plume of known concentration over realistic surfaces and perform the retrieval from the corresponding observed radiance under a given instrument to output the predicted methane concentration in each column. Our model calculates the noise-equivalent spectral radiance (NESR) as a function of incoming radiance and instrument parameters, such as integration time, detector size, quantum efficiency, readout noise, and spectral resolution, rather than prescribing the signal-to-noise ratio (SNR) as an independent variable. By varying the instrument and retrieval parameters, we can derive the associated precision error and bias from the retrieval. We also compare the tradeoff between the two most frequently used fitting windows in the 1.6 and 2.3 µm ranges.

Section 2 outlines the background on radiative transfer, followed by data and methodology on the forward model with realistic surface reflectances, instrument operators, and retrieval setups. The results and discussion are provided in Sect. 3. The final section contains the concluding remarks and future steps.

2 Data and methodology

For the sake of simplicity, we ignore the impact of atmospheric scattering, as Rayleigh scattering is negligible in the near-infrared, and the impact of aerosols is rather small compared to methane enhancements in the near-field of local sources. While aerosols can cause small systematic biases in the retrieved methane amount, their impact on measuring anomalies caused by methane plumes should be rather small. In addition, the precision error is not strongly affected by neglecting atmospheric scattering, and experience with previous moderate resolution methane mapping has shown that surface interferences are more crucial. In the absence of atmospheric scattering, and assuming a Lambertian surface, the reflected radiance as measured by an instrument at the top of the atmosphere in the nadir direction can be modeled as follows:

$$\begin{array}{}\text{(1)}& L_{\mathit{\lambda }}=I_{\mathrm{0},\mathit{\lambda }}\cdot r_{\mathit{\lambda }}\cdot T_{\mathit{\lambda }\uparrow }\cdot T_{\mathit{\lambda }\downarrow }\cdot {\displaystyle \frac{\cos \left(\mathrm{SZA}\right)}{\mathit{\pi }}},\end{array}$$

where I₀ stands for the incoming solar irradiance spectrum, T_λ↓ is the atmospheric transmission along the photon light path downwards to the surface, r_λ is the surface albedo, T_λ↑ is the transmission along the light path on the way up from the surface to the instrument, and SZA is the solar zenith angle. Figure 2 illustrates a schematic for Eq. (1). The subscript λ denotes the wavelength dependence of these variables. The multiplication in Eq. (1) is element-wise for each wavelength in the spectral range of interest.

Figure 2 A schematic for reflected sunlight from the Sun through the atmosphere to a spectrometer in space.

2.1 Incoming solar irradiance

We constructed I_0,λ by multiplying a continuum level spectrum with a high-resolution solar transmission spectrum that includes absorption features in the Sun's photosphere, which are the so-called Fraunhofer lines. These absorption features are caused by trace elements in the solar photosphere. The continuum spectrum is obtained from Meftah et al. (2018) with a 0.2 nm resolution. We fitted a third-order polynomial to this measured spectra in a 1.4–2.5 µm range to obtain a smooth continuum spectrum. A disk-integrated solar transmission spectrum is obtained from a tabulated line list compiled by Toon (2015). We interpolated the baseline and transmission spectra to a common 0.01 nm resolution grid and multiplied them to obtain a high-resolution solar irradiance I_0,λ.

2.2 Atmospheric transmission

The atmospheric transmission can be modeled using the Lambert–Beer law by dividing the atmosphere into vertical layers, each with a constant pressure, temperature, and gas number density. We calculate absorption cross sections for each layer using the HITRAN (high-resolution transmission molecular absorption database) spectral database (Gordon et al., 2017) and a Voigt line shape. The transmission then reads as follows:

$$\begin{array}{}\text{(2)}& T_{\mathit{\lambda }}=\exp \left(-\left(\sum _i\sum _k{n}_{i,k}{\mathit{\sigma }}_{i,k\left(\mathit{\lambda }\right)}\times {\text{AMF}}_k\right)\right),\end{array}$$

where i denotes the ith gas species, k is the kth layer, n is the vertical column density (molecule m⁻²), and σ is the gas absorption cross section (which is a function of pressure, P, and temperature, T). The air mass factor (AMF) per layer denotes the ratio of the integrated number concentration along the actual photon light path and the geometric vertical integration. In the absence of scattering, it is 1 cos (SZA), for the incoming light at SZA, and equal to 1, for the outgoing light as seen in nadir.

The transmission of the atmosphere is calculated from the background gas concentrations from the top of the atmosphere (TOA) to the surface. In addition, we consider the enhancements due to local gas emissions which we primarily considered to reside between the atmospheric boundary layer (BL) and the surface. The instrument is assumed to be located at the TOA.

For the background transmission, we divided the atmosphere into 72 layers and used an atmospheric profile for p and T from the Four Corners area (lat = 36.8${}^{{}^{\circ }}$, long = $-{\mathrm{108}}^{{}^{\circ }}$). We considered H₂O, CO₂, and CH₄ in the background. The concentration of H₂O is obtained from the Modern-Era Retrospective analysis for Research and Applications (MERRA) reanalysis (Rienecker et al., 2011) vertical profile. For simplicity, background CO₂ and CH₄ are set to volume mixing ratios of 400 ppm (parts per million) and 2000 ppb, respectively. For the gas enhancement within the BL, 300 vertical layers are used to model a simulated 3D methane plume enhancement from large eddy simulation (LES) output. The LES enables a realistic simulation of how methane concentrations from a point source evolve in space and time, as it generates the time-resolved three-dimensional CH₄ distribution in the boundary layer. The full description of the LES model setup for the CH₄ plume emanating from a point source can be found in Matheou and Bowman (2016), with the model parameterization and initialization detailed in Jongaramrungruang et al. (2019). We computed the gas absorption cross sections in each layer using an open-source Julia radiative transfer tool that calculates the cross section efficiently using the HITRAN database and graphics processing unit (GPU) capability (Gordon et al., 2017). A Voigt absorption line shape is used in our study. We note that we used 300 vertical layers in the BL in the full forward model to simulate the observed outgoing spectra, but we will use a much smaller number of layers in the retrieval step (more details in Sect. 2.4).

2.2.1 Surface reflectance

To analyze the impact of surface spectral features, we compiled a database of different surface albedos from the ECOSTRESS spectral library (Meerdink et al., 2019), to investigate the impact on our traditional retrieval technique that we use from space, in which the surface is typically characterized by a low-order polynomial in wavelength. In fact, many spectroscopic measurement techniques rely on the fact that atmospheric features exhibit sharp absorption features while the surfaces are spectrally smooth (Platt and Stutz, 2008). From a physical perspective, this is related to more rapid quenching of an excited state in solids or liquids and the suppression of rotational energy levels. At the same time, this separation of high-frequency atmospheric features from low-frequency surface features is at the core of our study, as instruments such as AVIRIS(-NG) have spectral resolutions that can blur the separation between frequencies, allowing the surface features to alias into methane retrievals.

The compiled database contains more than 2000 surfaces over five main categories of rock, soil, mineral, photosynthetic and non-photosynthetic vegetation, and artificial construction materials. Examples of these surface albedos near the 2.3 µm CH₄ absorption range are shown in Fig. 3. The typical spectral resolution of the database is 2 nm, and we resampled all spectra to a common grid, which is subsequently used with spline interpolation in our high-resolution forward model at 0.01 nm.

Figure 3 Albedo spectral variations near 2.3 µm from distinct surface reflectances in our database.

Based on this database, we can create an arbitrarily diverse set of surfaces underlying the simulated 3D methane field. For example, a checkerboard-style tile consisting of 3×10 surfaces can be constructed before we overlay the 3D plume. A retrieval can then be applied for each pixel across the source to visualize the impact on precision error and bias caused by different surfaces. We also used LANDSAT data (Wulder et al., 2019) to represent a natural distribution of surfaces in the Durango area, Colorado. Although it is a 30 m resolution, we can up-sample the surface grids to 5 m each. At each pixel, we matched the measured surface albedo from LANDSAT observations at 0.48, 0.56, 0.65, 0.87, 1.61, and 2.20 µm to the closest possible surface in our database and queried its full albedo spectra for our simulations.

2.3 Instrument operators

2.3.1 Convolution and observed spectra

The actual observed spectrum that is recorded by the instrument is the convolution of the high-resolution incident light L_λ with the instrument line shape, denoted here as the instrument kernel. This convolution is performed in the intensity space given by the following:

$$\begin{array}{}\text{(3)}& \langle L\left(\mathit{\lambda }\right)\rangle =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}L\left({\mathit{\lambda }}^{\prime }\right)\mathit{\varphi }(\mathit{\lambda }-{\mathit{\lambda }}^{\prime })\mathrm{d}{\mathit{\lambda }}^{\prime },\end{array}$$

where L(λ) is the incident spectra on the device, and 〈〉 denotes the convolution with the instrument kernel ϕ. The convolved spectrum can then be interpolated and resampled to the output wavelength grids (λ_out) of the instrument in an observing spectral range of interest. These output spectra are used as the measurement vector in the retrieval process (more details in Sect. 2.4).

The instrument kernel ϕ is modeled using a Gaussian distribution with zero mean and a given full width at half maximum (FWHM). In our experiment, we treated the FWHM as an independent variable that varies between 0.04 and 10.0 nm. The FWHM is a key property of an instrument that determines what spectral variations can be resolved. For instance, if the spectral resolution is coarser than the rotational fine structure of a vibrational–rotational absorption band, the P and R branches of this band will appear as just two separate broadband absorption features. The spectral sampling interval (SSI) varied accordingly with FWHM. Here, we use two cases of SSI equal to FWHM $/$ 2.5 (near-Nyquist sampling in atmospheric sounders) and FWHM/1.0 (critical sampling as in AVIRIS-type imaging spectrometers).

2.3.2 Noise-equivalent spectral radiance

Towards designing the optimal instrument, we have to evaluate the trade space of the spectral resolution, spatial resolution, and detector characteristics. Previous studies have evaluated the trade space between the signal-to-noise ratio (SNR) and spectral resolution, treating the SNR as an independent variable when varying the spectral resolution (Thorpe et al., 2016; Cusworth et al., 2019; Ayasse et al., 2019). However, the SNR deteriorates at a higher spectral resolution, as fewer photons are being counted by each detector pixel. To evaluate this properly, we start working from an instrument model directly where the noise is a function of incoming radiance and parameters, such as integration time, F number, detector pixel, and then the spectral resolution, as described in Strandgren et al. (2020). In this approach, the SNR will be a dependent variable based on our instrument specifications and the actual observed radiance.

The electronic signal measured by each detector pixel can be expressed as follows (Strandgren et al., 2020):

$$\begin{array}{}\text{(4)}& S=\langle L_{\mathit{\lambda }}\rangle {\displaystyle \frac{\mathit{\pi }{A}_{\det }}{\mathrm{4}{f}_{\mathrm{num}}^{\mathrm{2}}}}\cdot \mathit{\eta }\cdot Q_{\mathrm{e}}\cdot \mathrm{\Delta }\mathit{\lambda }\cdot t_{\mathrm{int}},\end{array}$$

where 〈L_λ〉 is our simulated radiance, A_det is the detector pixel area, f_num is the instrument's F number, η is the optical efficiency of the spectrometer, Q_e is the quantum efficiency of the detector, Δλ is the SSI, and t_int is the integration time.

For the NESR, we consider two dominant noise terms, namely shot noise (proportional to $\sqrt{S}$) and effective readout noise σ_ro, as follows:

$$\begin{array}{}\text{(5)}& \mathrm{NESR}={\mathit{\sigma }}_{L_{\mathit{\lambda }}}=\sqrt{S+{\mathit{\sigma }}_{\mathrm{ro}}^{\mathrm{2}}}.\end{array}$$

For the analyses in this study, we varied FWHM from 0.04 to 10.0 nm and integration time from 5 to 105 ms with other default parameter setups, as described in Table 1, unless stated specifically otherwise.

Table 1 A table for default parameter settings in our simulations.

2.4 Retrieval setup

2.4.1 Forward model

The retrieval forward model is similar to that described earlier in previous sections. The main difference is that we now treat the methane concentration in the boundary layer as elements of the state vector that we want to retrieve, as the methane profile is unknown during the observation. A much smaller number of layers to represent methane enhancements within the boundary layers is used. This is an important consideration when dealing with moderate spectral resolution, since the information content is not high enough to discriminate between different layers. Polynomial terms are also used to represent the change in surface reflectance with wavelength. The forward model can be written mathematically as follows:

$$\begin{array}{}\text{(6)}& {\mathit{y}}_{\mathit{\lambda }}=F\left(\mathit{x}\right)=\langle L_{{\mathit{\lambda }}_{\mathrm{0}}}\cdot T_{\mathit{\lambda }\uparrow }\cdot T_{\mathit{\lambda }\downarrow }\rangle \underset{\approx r_{\mathit{\lambda }}}{\underbrace{\sum _{d=\mathrm{1}}^D{a}_d{P}^d\left(\mathit{\lambda }\right)}},\end{array}$$

where y is the measurement vector, P^d(λ) is a Legendre polynomial term at degree d, a_d is a coefficient for each P^d(λ), and D is the number of polynomial degrees used in the retrieval. To evaluate the polynomial degree within the fitting range, we converted the wavelength range within the fitting window to span −1 through 1. Here, the state vector x consists of the vertical column density (molecule m⁻²) of the respective gases in different layers and the polynomial coefficients accounting for low-frequency surface features.

In our experiment, we set the number of water vapor and methane concentration layers to be 2 and 10, respectively. We vary the number of the polynomial degree, D, to investigate its ability to disentangle the surface spectral variations from the atmospheric features while simultaneously considering different instrument spectral resolutions. Spectrally smooth surface albedos might only require a few polynomial coefficients, while more complex surface features within the 260 nm fitting window can require more than 25 coefficients, as will be shown later.

2.4.2 Optimal estimation

To solve for the optimal state vector from this nonlinear system, we used the iterative maximum a posteriori differential optical absorption spectroscopy (IMAP-DOAS) approach (Frankenberg et al., 2005a). Based on maximizing the a posteriori probability density function as introduced to the atmospheric community by Rodgers (2000), the iterative solution can be written as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{x}}_{i+\mathrm{1}}& ={\mathit{x}}_a+({\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i+{{\mathbf{S}}_{\mathbf{a}}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\\ & \cdot [\mathit{y}-F({\mathit{x}}_i)+{\mathbf{K}}_i({\mathit{x}}_i-{\mathit{x}}_a\left)\right],\end{array}\end{array}$$

where x_i is the state vector at the ith iteration, x_a is a priori state vector, S_ϵ is the measurement error covariance matrix, S_a is the a priori covariance matrix, and K_i is the Jacobian of the forward model evaluated at x_i. F(x_i) stands for a forward model at each x_i. When consecutive changes in the reduced χ² of the fit drop below a tolerance level of 10⁻³, we stop the iterations. The Jacobian matrix K_i is computed analytically in each ith iteration, using automatic differentiation techniques. The a priori covariance matrix S_a helps constrain the fit based on the possible range of concentration (thus, S_a is a square matrix with a size equal to the length of x_i). Here, we use loose prior constraints, thus having no significant impact on the retrieved total columns or posterior errors. The measurement error covariance matrix S_ϵ is a matrix in which the diagonal elements are the estimated variances of instrument noise at the observed wavelength grid λ_out (thus, S_ϵ is a square matrix with a size equal to the length of λ_out). These variances are computed from the instrument noise model outlined in Sect. 2.3.2.

2.4.3 Error estimations

The a posteriori error estimate $\widehat{\mathbf{S}}=({\mathbf{K}}^T{\mathbf{S}}_{\mathit{ϵ}}\mathbf{K}+{\mathbf{S}}_a^{-\mathrm{1}}{)}^{-\mathrm{1}}$ provides the full error covariance matrix of the retrieved state vector $\widehat{\mathit{x}}$. The quantity of interest in our application for CH₄ detection and quantification is the total column concentration of CH₄. To obtain this quantity from our retrieval, we can find the summation of the state vector elements over indices corresponding to CH₄. We can, thus, define a summation operator h, of the same size as the state vector, filled with ones where the CH₄ state vector elements are located and the zeros are elsewhere. The summation of the total column is readily derived as ${\mathit{h}}^T\widehat{\mathit{x}}$, and the variance of the total column is computed as ${\mathit{h}}^T\widehat{\mathbf{S}}\mathit{h}$ (Rodgers, 2000). The bias error is obtained as the difference between the best estimated value in the absence of instrument noise and the true CH₄ vertical column enhancement.

2.5 LES CH₄ plumes

The LES is used to generate the time-resolved three-dimensional CH₄ distribution in the boundary layer. This provides a realistic distribution of how methane concentrations from a point source on the ground evolve across the area. The flow in the boundary layer is driven by a constant geostrophic wind in the x direction, influencing the shape of the plume. We conducted LES experiments with a geostrophic wind speed of 4 m s⁻¹, with a range of emission rates from 50 to 5000 kg h⁻¹.

To simulate the reflected sunlight over the area of a CH₄ emission source, we use the output from LES that has a realistic 3D concentration field of CH₄ at a prescribed emission rate and overlay this on top of surface tiles of different albedos. Within each tile, the albedo is chosen from a distinct surface in our database. Figure 4 illustrates this conceptual setup.

Figure 4 A schematic showing the setup for a methane plume over different surface tiles, each of which contains one surface albedo from our database.

At each spatial pixel, the high-resolution outgoing reflected radiance can be calculated and subsequently be convolved with the instrument kernel. This yields the simulated observed radiance at each pixel that we implemented as a measurement vector. Accordingly, the pixel radiances are then converted to signal strength in electrons and respective noise levels using our instrument model. The IMAP-DOAS algorithm is applied to retrieve the column CH₄ and provide error estimates. We can vary the instrument parameters and the number of polynomial degrees in the retrieval to explore their relationship with the associated errors.

3 Results and discussion

3.1 Simulated high-resolution and observed radiance

As a first step towards understanding the effect of the instrument spectral resolution, we simulated high-resolution spectrum 〈L_λ〉 in the 1.6 and 2.3 µm bands and their corresponding Jacobians for instruments, with FWHM of 0.2, 1.5, 5.0 and 10.0 nm. This simulation is based on a CH₄ concentration profile near an origin of the methane plume, with an emission rate of 200 kg h⁻¹, over a construction concrete (albedo index 1 in Fig. 3), with the simulations shown in Fig. 5. As the instrument's FWHM increases, individual absorption lines are increasingly blurred, and the less high-resolution absorption features are recorded. The Jacobian represents the change in radiance with respect to the changes in gas concentration. Here we show the Jacobian for CH₄ and H₂O close to the ground. For the 1.6 µm band, the total radiance is in the range of 20–30 mW m⁻² nm⁻¹ sr⁻¹, and we observed the strongest absorption feature between 1.66 and 1.67 µm, covering the Q branch in the 2ν³ band. On the other hand, for the 2.3 µm band, the radiance varies from around 10 to 3 mW m⁻² nm⁻¹ sr⁻¹ across the band. The CH₄ absorption features are much stronger and more prominent over a wider range from 2.2 to 2.4 µm. To compare the two fitting windows in terms of their effectiveness in CH₄ measurements, we explore errors associated with the CH₄ retrieval using each of these bands in the next section.

Figure 5 Simulated reflected solar radiance, through a plume of 200 kg h⁻¹, over a concrete surface, as observed by instruments of different FWHM in the two fitting windows of 1.6 µm (a–c) and 2.3 µm (d–f). The gray background is the originally calculated spectra at 0.01 nm.

3.2 Comparisons of two CH₄ fitting windows

We investigate the comparison between using two different CH₄ fitting windows near 1.6 and 2.3 µm. At 1.6 µm, the incoming solar irradiance is higher, which could enhance the signal-to-noise ratio. However, at 2.3 µm the CH₄ absorption features are more numerous and prominent over a broader wavelength range, which should increase sensitivity for the CH₄ retrieval. Thus, there is a potential tradeoff between the advantages and disadvantages of both retrieval windows. It is not immediately obvious which fitting window would result in a lower precision error. This could also depend on the spectral resolution, as the spectral fine structure changes with increasing FWHM differently in both windows. Understanding this tradeoff will help guide the development of future instruments.

Resolving which band is better suited for minimizing CH₄ retrieval errors also depends on typical surface albedos at 1.6 and 2.3 µm. We visualized the relative values of surface albedo at 1.6 and 2.3 µm across all surfaces in our database colored by their main surface categories in Fig. 6. Most surfaces lie in the region where their albedos at 1.6 and 2.3 µm are relatively equal (along the 1:1 line). Only a few surfaces have a much higher albedo at 2.3 µm compared to that at 1.6 µm, while there is a cluster of low albedo surfaces for which the 1.6 µm albedo is about twice as high as at 2.3 µm. To compare the two fitting windows, we take two extreme surface examples, i.e., one with an albedo of 0.87 at 1.6 µm and 0.14 at 2.3 µm and another one with and albedo of 0.21 at 1.6 µm and 0.79 at 2.3 µm. We also consider a representative surface with an equal albedo of 0.51 at both 1.6 and 2.3 µm.

Figure 6 Scatterplot showing the relative values of reflectance at 1.6 vs. 2.3 µm from our database of different surface albedos from the ECOSTRESS spectral library. Each point is a distinct surface, and the color shows the type to which it belongs. Cross marks represent three example surfaces that we used in the comparison analyses.

In Fig. 7, we compared the precision errors based on the two fitting windows over the three surfaces using a polynomial degree of 25 and using fitting windows as displayed in Fig. 5. In this simulation, only the FWHM is varied, and an oversampling of 2.5 is used (i.e. FWHM = 2.5 ⋅ SSI). Other instrument parameters are set according to Table 1. The corresponding NESR of each detector pixel is then computed using the instrument noise model, enabling us to compute the changing noise levels with instrument resolution. We note that the number of detector pixels used varies considerably in this simulation, as they correspond to the window length (275 nm in the 2.3 µm band) divided by SSI. Thus, a FWHM smaller than 0.65 in the 2.3 µm band would require more than 1000 detector pixels, which can be hard to achieve from a detector point of view, especially at fast readout rates as required for high spatial resolution. This can preclude very fine spectral resolution if a large spectral bandwidth is required, which can be beneficial, especially if it allows additional species to be measured (such as CO₂). In addition, at a coarser spectral resolution, the full well capacity of the detector might be reached, which puts an upper limit of maximum SNR values per detector pixels (maximum SNR of 1000 for a full well capacity of 10⁶ electrons).

Figure 7 A plot showing precision errors from the CH₄ column retrieval in the 1.6 µm band (dotted) and the 2.3 µm band (solid) over three distinct surfaces (colored), using instruments with different FWHM (in the unit of mole m⁻²).

Here, we observe that the precision error is actually not monotonically improving with finer spectral resolution as the SNR deteriorates in these cases. For the 2.3 µm band, the optimum is actually around 0.2 nm, which may appear surprising. For darker surfaces, the precision, using a 1 nm resolution, can be equal to an instrument with 0.05 nm resolution. Interestingly, the precision errors using the 1.6 µm band deteriorate much more with increasing FWHM than using the 2.3 µm band. The reason for this is related to the methane band structure in both windows, as can be seen in the Jacobians (see Fig. 5). In the 1.6 µm band, most of the fine structure in the P and R branches is lost once the FWHM is coarser than the separation of the rotational fine structure, leaving only the Q branch as spectrally distinct, even at lower spectral resolution. This could also be seen in SCIAMACHY, for which the FWHM in the 1.6 µm band was relatively low, namely about 1.3 nm (Frankenberg et al., 2005b). With the loss of spectral fine structure, the precision deteriorates more rapidly, as low frequency variations in the Jacobians can be confused with surface albedo features, thus not directly constraining methane abundances. The situation in the 2.3 µm band is different as the absorption structures are very different and cover many more methane features that are irregularly spaced. Thus, some unique bulk methane absorption features persist even at spectral resolutions as coarse as 10 nm. In fact, this is the reason methane could be observed with the imaging spectrometer at similar resolutions (Thorpe et al., 2014), enabling high spatial resolution mapping from airborne instruments.

Overall, our retrievals in the 2.3 µm band consistently yield lower precision error compared to that from the 1.6 µm band. This difference is most highlighted for a surface with a much higher albedo at 2.3 µm compared to 1.6 µm. This is because a decrease in solar irradiance at 2.3 µm is now compensated for by increased albedos, resulting in observed radiances being relatively close to those at 1.6 µm, allowing for the effect of absorption depth and structure to be the sole driving force for a better performance. This effect is still seen for most typical surfaces of equal albedos at both 1.6 and 2.3 µm, where the precision error from the 2.3 µm band is consistently lower than that from the 1.6 µm regions for instruments with any FWHM. The only scenario where the retrieval at 1.6 µm could perform better is under an extreme example, where a surface has a much higher albedo at 1.6 µm compared to at 2.3 µm and when FWHM is lower than 0.2 nm. These results indicate that the stronger and broader absorption of CH₄ in the 2.3 µm fitting window plays a more dominant role in the retrieval performance compared to the stronger solar irradiance (and sometimes higher albedo) in the 1.6 µm case. Since the latter condition for a better performance in the 1.6 µm band is much more unlikely, we focus the following analysis on the 2.3 µm fitting window.

3.3 Spectral fit and error analysis with various instruments

Here, we evaluate the instrument performance of hypothetical yet realistic spectrometers covering the 2.3 µm range. If we restrict ourselves to the number of spectral pixels in a fast detector, as used for AVIRIS-NG (480 spectral pixels), we can achieve a FWHM of around 1.5 nm, minimizing the precision error under the hard constraint of a limited detector size. It would also still allow for joint retrievals of CO₂ at 2 µm, as envisioned in Strandgren et al. (2020). This resolution allows us to still resolve significantly more spectral fine structure than current measurements at 5 or 10 nm resolution, which often exhibit retrieval interferences with surface features, as seen in Fig. 1. Our primary focus here is to quantify the impact of the spectral resolution on the ability to unequivocally separate surface spectral features from those in the atmosphere.

To demonstrate the impact of surface reflectance on observed spectra, Fig. 8a shows two reflected spectra through exactly the same CH₄ concentration profile as observed by the same instrument with FWHM of 1.5 nm but over construction concrete as compared to a vegetation surface. Evidently, not only do the observed spectra change in absolute value, but their spectral variations are also different within the fitting window and are more complex for vegetation than for concrete. This exemplifies an important role that surface albedo plays in the retrieval, potentially interfering with the methane absorption lines. To further validate this point, Fig. 8b shows examples of residuals from the best spectral fits for each of these two spectra, using IMAP-DOAS as described in Sect. 2.4.2, with a polynomial degree of 11. The noise level indicates the theoretical 1σ noise level, which is expected from the instrument. In Fig. 8a, these residuals are hardly visible as they are close to the noise level. Clearly, the fit quality is different between the two surfaces, as evidenced by a higher residual for reflected sunlight from the vegetation surface compared to a sample construction concrete. However, differences are subtle and might not be detectable if noise levels are high or fewer detector pixels are available. In practice, given an observed spectrum, our retrieval needs to differentiate the atmospheric feature from the surface feature in order to obtain the best estimate of methane enhancement with minimal bias and precision error. The performance of the retrieval will be influenced by the observing instrument and the representation of surface reflectance in our forward model.

Figure 8 (a) The observed spectra through exactly the same CH₄ source over two different surfaces, namely construction concrete and vegetation, as observed by an instrument with a FWHM of 1.5 nm and the best spectral fit. (b) The associated spectral residuals and an expected 1σ noise level for the instrument measurement.

In this section, we performed the retrieval error analysis for different instrument parameters, namely the FWHM and exposure (integration) times. Since the surface albedo is another important factor in the retrieval performance, as shown in Fig. 8, we also explored using different degrees in the polynomial terms in our forward model (see Sect. 2.4.1), as a large number of polynomial degrees is required to capture the effect of surface albedo spectral variability (e.g., a polynomial degree of 11 still caused subtle yet systematic residuals for vegetation). For a given choice of instrument specification and polynomial degree, the IMAP-DOAS algorithm is performed to predict the column methane concentration and, thereby, derive error estimates. We kept the true methane concentration fixed to a given vertical profile near the source emission with a flux rate of 200 kg h⁻¹. First, to illustrate potential bias that could arise from a variety of surfaces in real-world scenarios, we randomly selected 100 surfaces from our database and used each of them as a surface underlying the CH₄ column. By retrieving the observed radiance using different instrument FWHMs and different numbers of polynomial degrees, we analyzed the range of resulting biases for each case. The bias distribution across the 100 surfaces is shown in Fig. 9. When a polynomial degree of 5 is used in our retrieval, the range of biases observed is (−0.2, 0.8) mole m⁻² compared to just within (−0.15, 0.15) and (−0.02, 0.02) mole m⁻², with a polynomial degree of 11 and 25, respectively. This result suggests that some surfaces interfere strongly with the retrieval, leading to very high biases (>100 % from typical anomalies in plumes) when a low number of polynomial degrees is used, as this causes a forward model error. This implies that the use of a higher polynomial degree has a significant consequence in terms of minimizing the bias in our retrieval results, especially when we have no prior information about surfaces in the vicinity of emission sources. A polynomial degree of about 25 seems to capture most but not all surface effects. However, an instrument such as AVIRIS(-NG) only has 26 (52) detector pixels covering the entire fitting window, thus not allowing us to use such high polynomial degrees, as it would render the problem underdetermined. This clearly illustrates the problem in separating surface and atmospheric features at coarse spectral resolutions, as the problem becomes increasingly ill-posed with a coarser spectral resolution. This is an integral part in obtaining the reliable detection and quantification of local methane sources at a global scale, as is also shown later in Sect. 3.4.1. In general, an instrument with a smaller FWHM leads to smaller observed biases, as expected, by its ability to capture more high-frequency CH₄ absorption features. Furthermore, with an SSI of FWHM $/$ 2.5, biases using higher polynomial degrees show a smaller range compared to the case with an SSI = FWHM, as surface features are harder to discern if atmospheric features are not oversampled. In the following analyses, we, therefore, primarily show the results from the case of SSI = FWHM $/$ 2.5, unless otherwise stated. This also fulfills the Nyquist sampling requirement for typical atmospheric retrievals, which might sometime involve spectral shifts. The impact of spectral shifts is ignored here but would be another reason for both oversampling and higher spectral resolution.

Figure 9 Plots showing a range of biases that occur over 100 randomly sampled surfaces when different FWHM, SSI, and polynomial degrees are used.

Next, we investigated how the precision error varies with the instrument parameters such as FWHM and the exposure (integration) time. While FWHM governs the shape of the methane Jacobians, the exposure time is an important factor determining potential spatial resolution and SNR. Figure 10 shows this relationship using a fixed exposure time of 20 ms, with varying polynomial degrees on Fig. 10a and with fixed polynomial degree (25) and varying exposure time in Fig. 10b. Other parameters are set according to Table 1, and the underlying surface is a construction concrete (albedo index 1 in Fig. 3). If we vary the polynomial degree, the impact on precision is negligible for FWHM <0.5 nm, very small for FWHM <1.5 nm, but diverging from one another for FWHM >2 nm. The reason for this effect is that the polynomial degree determines which spectral variations can be purely attributed to methane and which might be caused by surface features, which eliminates its use to constrain methane. At a fine spectral resolution, the methane fit is mainly driven by the atmospheric high-frequency structure. Thus, the polynomial degree plays no significant role. At a coarser spectral resolution, most atmospheric features are blurred and can be partially confused with the surface, causing a divergence in the precision error with increasing FWHM. A 1.5 nm FWHM still allows us to sample sufficient atmospheric fine structure to minimize the impact of surface interferences and will achieve a sub-percent precision error for a wide range of exposure times.

Figure 10 (a) Precision error as a function of instrument FWHM for different numbers of polynomial degrees used in the retrieval (with an exposure time of 20 ms). (b) Precision error as a function of instrument FWHM for different exposure times (with a polynomial degree of 25).

If we vary the exposure times, the precision error decreases with the increase in the exposure time as the device can collect more photons resulting in an overall stronger signal. There exists a range of FWHM values that minimize the precision error for each exposure time. Generally, the range lies between 0.1 and 0.3 nm. At low exposure times, readout noise becomes increasingly important and leads to a larger precision change when exposure times are varied. Thus, it is vital that an appropriate value of FWHM is chosen in order to achieve low precision error while we use a high degree of the polynomial, such as 25, in our retrieval to simultaneously reduce the bias from surface interferences. It is interesting to note that the FWHM with the best precision moves towards higher FWHM with decreasing exposure times, which is a consequence of the increasing role of the readout noise at the detector.

For example, using a FWHM of 1.5 nm and an exposure time of 25 ms, a precision error of 0.007 mole m⁻² can be achieved, which is about 1 % of the background total column amount of 0.7 mole m⁻². In low Earth orbit, with a satellite speed of 7 km s⁻¹, 25 ms corresponds to a spatial resolution of 175 m. Spacecraft nodding would allow us to slow down the effective ground speed by about a factor of 10, rendering a 1 % total column precision for <20 m spatial resolution feasible from space using existing fast readout short-wave infrared (SWIR) detectors. This would be equivalent to measuring a pure methane layer of only 0.16 mm thickness.

3.4 2D retrieval over realistic surfaces

In Sects. 3.1, 3.2, and 3.3, we analyzed the impact of surface interferences and instrument specifications on the quality of the CH₄ retrieval. The choice of parameters such as FWHM, exposure time, and the number of polynomial degrees leads to significantly different precision errors and biases. These errors in each retrieval column ultimately affect the detection and quantification of CH₄ source in 2D scenes observed over various geographical areas across the globe. In this section, we illustrate how the retrieved CH₄ plume appears using different instruments and retrieval choices.

3.4.1 Occurrence of false positives and false negatives

To explore the 2D pattern of a retrieved methane plume over a variety of realistic surfaces, we overlaid an LES simulated CH₄ plume on top of a checkerboard-style land consisting of 30 surfaces. We then apply the instrument kernel and IMAP-DOAS algorithm to retrieve column CH₄, pixel by pixel, across the 2D scene. A combination of a FWHM of 1.5 and 5.0 nm, with an SSI of FWHM $/$ 2.5, and polynomial degrees of 5, 11, 25, and 50 are adopted. A case of SSI = FWHM is also used for a FWHM of 5.0 nm, which is equivalent to the AVIRIS-NG instrument. The emission rate of the CH₄ source is equal to 200 kg h⁻¹. The spatial distribution of the true CH₄ distribution is shown in Fig. 11, showing local enhancements than can exceed 10 % of the total background atmospheric column. The corresponding retrievals under different surfaces and instrument scenarios are shown in Fig. 12. The deviation of the predicted CH₄ from the true column value in each pixel is a combination of precision error and bias. In this result, the overall mean enhancement that emerges over each surface type, in contrast to the true plume in Fig. 11, could be interpreted as a retrieval bias, while the presence of a speckle-like texture over each surface can be viewed as the retrieval precision error driven by instrument noise for a given surface albedo (larger for dark surfaces). Meanwhile, the bias is related to a systematic shift in the retrieved methane column enhancement from the true value due to surface interference in resolving the methane absorption features. In general, based on the visualization in these 2D plots, the more the retrieved enhancement scene resembles the true CH₄ concentration map, the better the performance of the instrument and retrieval is. Specifically, when a low polynomial degree of 5 was used in the retrieval, we observed significant retrieval biases (both positive and negative) over various surfaces. Evidently, these biases can act to deceive the true location and enhancement of an actual methane plume, especially if the surfaces have elongated shapes like a plume (not like a checkerboard here). As the number of polynomial degrees in the retrieval increases, the level of biases decreases over the scene, enabling the actual methane plume enhancement to be better identified. This is also manifested in the reduced χ² values which describe how well an observed spectrum was fitted (i.e., the smaller, the better). As shown in Fig. 13, the value of reduced χ² drops from 5 to become increasingly closer to 1 when a polynomial degree is changed from 5 to 11, 25, and 50. We note that the 10 surfaces in the bottom row of this checkerboard-style tile, from left to right, are the 10 surfaces, with their albedo spectral variations shown earlier in Fig. 3. The pink quartzite (sixth on the list of these 10 surfaces) is an extreme case, where we see an unusually strong variation near 2.2 µm, resulting in a persisting bias even at a polynomial degree of 25. Nevertheless, generally, when a polynomial degree of 25 is used, most of the biases across surfaces seem to disappear, but surprisingly complex surfaces, such as quartzite, can occur across various natural landscapes and artificial surfaces in cities.

Figure 11 A nadir view of a simulated CH₄ plume from the LES with a prescribed flux rate of 200 kg h⁻¹. The background wind speed is 4 m/s⁻¹ along the x axis.

Figure 12 Retrieved CH₄ column over 30 different surfaces with varying instrument scenarios. The top row shows the results for an instrument with a FWHM of 1.5 nm, while the second row shows the results for an instrument with a FWHM of 5.0 nm. Both the top and middle rows have SSI = FWHM $/$ 2.5. The bottom row shows the results for an instrument with FWHM of 5.0 nm and SSI = FWHM. The biases drop significantly as the polynomial degrees increase.

Figure 13 A plot showing the values of reduced χ² from the retrieval for different instrument FWHM and polynomial degrees. As the polynomial degree becomes higher, the reduced χ² decreases, implying a better spectral fit.

At the same time, using higher polynomial degrees results in higher precision errors, as can be seen from the speckle-like texture in the retrieved scene. To illustrate how precision deteriorates as we increase the polynomial degree for a given instrument FWHM, Fig. 14 shows the rise in precision error when using polynomial degrees of 50, 25, and 11 relative to using the polynomial degree of 5. We can clearly see that the precision error deteriorates with a higher polynomial degree, particularly at coarser FWHM. This is consistent with what we observed earlier in Sect. 3.3. Based on this analysis, we found that using an instrument with a FWHM of 1.5 nm would allow for higher polynomial degrees, such as 25, to be utilized with a relatively small increase in the precision error. For an instrument with a FWHM of 5 nm, using a polynomial degree of 25 or 50 results in a larger precision error increase, underlining the potential problems that could occur when complex surface albedo features exist. The main physical reason for the deterioration in the precision of low spectral resolution instruments is that a lower degree polynomial in the fitting routine is equivalent to a hard a priori constraint suggesting that only spectrally smooth surfaces exist. The retrieval itself then attributes some of the broadband variations in the methane Jacobian to be only attributable to changes in methane and not surface albedos, thus providing a tighter constraint on the methane abundances. For higher resolution instruments, most of the information content for methane is located within the fine structure of the methane absorption lines and less on the broadband variations, causing a much smaller increase in the precision errors if higher polynomial degrees are used.

Figure 14 The increase in precision errors when the polynomial degrees of 50, 25, and 11 are used, as compared to the polynomial degree of 5 case. The increase is computed for a given FWHM and SSI, according to each row.

These results indicate that a narrow FWHM (such as 1.5 nm) and a high number of polynomial degree (at least 25) are needed to reduce both precision errors and biases due to surface interference. In this way, we can obtain a higher quality retrieved scene in order to effectively identify and quantify emission sources over the majority of the surfaces. If low polynomials are used, high biases are likely to occur over certain surface types. These can cause false positives or negatives in the observational systems, complicating the analysis of the locations and the emission rates of the CH₄ sources. This can cause a significant problem for both human analysts and AI models (Jongramrungruang et al., 2021) that rely on the spatial distribution of observed enhancement to make predictions.

To further show an example over a real-world high emission area, we queried a realistic surface distribution over a well pad in the Durango area, Colorado, from LANDSAT. At each location, this data set provides surface albedos at the wavelengths of 0.48, 0.56, 0.65, 0.87, 1.61, and 2.20 µm. The RGB image of this particular location and its corresponding albedo near 2.2 µm are shown in Fig. 15. Based on the albedos in the seven bands available in this LANDSAT scene at each pixel, we found the best-matching surface in our high-solution surface database from the ECOSTRESS spectral library. We used this surface to simulate a semi-realistic CH₄ scene with an emission emerging from the ground in the vicinity of a well pad. The retrieved CH₄ enhancement is shown in Fig. 16. Once again, when an instrument with a FWHM of 5.0 nm is used in conjunction with a low polynomial degree of 5, high biases occur over most surfaces across the area. In particular, oil shale and sandstone are the two surfaces that exhibit strong surface interference, as observed in Fig. 16f. Evidently, the resulting bias occurs at a level that dwarfs the true plume enhancement, rendering it impossible to identify the location of the emission sources, let alone the ability to obtain an accurate quantification of total emission in the area. However, by using a higher polynomial degree of 25, the biases across surfaces are greatly reduced, albeit with slightly higher noise. Nonetheless, this noise is reduced when FWHM decreases from 5.0 to 1.5 nm. We note that this reduction in precision error when we decrease the instrument with a FWHM from 5.0 to 1.5 nm could be even more apparent over some other surfaces with lower albedos compared to the ones in this scene. Visually, we can already see that the actual CH₄ plume in terms of its location and strength can be much more easily identified and distinguished from the surface artifacts with the FWHM of 1.5 nm and polynomial degree of 25. This finding demonstrates how achieving low bias and precision error in the observation and the retrieval process across diverse surfaces profoundly benefit the detection and quantification of true CH₄ sources. This analysis provides insight to how future instruments can be designed to enable an effective and accurate CH₄ source detection and quantification across the globe. In the next section, we further illustrate the retrieval performance when a local CH₄ plume of various emission rates is observed by an instrument with a higher spatial resolution (such as 30 m).

Figure 15 RGB and 2.2 µm albedo images of a realistic surface distribution in the Durango area, Colorado, from LANDSAT.

Figure 16 Retrieved CH₄ column over a well pad in the Durango, Colorado, area with varying instrument FWHM and polynomial degrees in the retrieval. The realistic surface distribution is based on LANDSAT data. The use of a polynomial degree of 25 and a FWHM of 1.5 nm can reduce biases and precision error over this scene.

3.4.2 Effect of spatial resolutions and flux rates

In the previous section, we have shown a 2D retrieved scene, assuming that an observing instrument has a spatial resolution of 5 m. In this section, we repeated the 2D scene retrieval analysis with a spatial resolution of 30 m by averaging the reflected sunlight through a CH₄ plume simulated at 5 m spatial resolution into 30 m spatial resolution prior to applying an instrument operator. We present this analysis with this design consideration for the 30 m spatial resolution and the exposure time of 40 ms to evaluate the potential of future spectrometers on board satellites in the coming years. The 2D scenes retrieved at a 30 m spatial resolution by instruments of different FWHM and polynomial degrees are illustrated in Fig. 17. The choice of a FWHM of 1.5 nm and a polynomial degree of 25 remains very effective in removing surface biases across the scene, and the location of CH₄ plume can be distinguished. It is important to note, however, that the retrieved CH₄ column concentration near the source pixels becomes more diluted as the spatial resolution decreases. This is because the local CH₄ plume distribution at an emission rate of, e.g., 200 kg h⁻¹ varies greatly on scales of just a few meters. Having demonstrated that an instrument with a FWHM of 1.5 nm and a polynomial degree of 25 can significantly reduce precision error and biases due to surface interference, we use this setup to investigate what the 2D retrieved scenes look like for sources of different emission rates to understand the lower limit of CH₄ emission rates that can potentially still be detected.

Figure 17 A plot similar to Fig. 16 but at 30 m resolution. Retrieved CH₄ column over a well pad in the Durango, Colorado, area with varying instrument FWHM and polynomial degrees in the retrieval. The realistic surface distribution is based on LANDSAT data. The use of a polynomial degree of 25 and a FWHM of 1.5 nm can reduce biases and precision errors across the scene.

Figure 18 Retrieved CH₄ column over a well pad in the Durango, Colorado, area for plumes of various emission rates. The spatial resolution is 30 m, the instrument FWHM is 1.5 nm, and the polynomial degree in the retrieval is 25. The realistic surface distribution is based on LANDSAT data. Column enhancement in the vicinity of the CH₄ plume is increasingly visible as the source emission rate becomes larger.

Figure 19 Ratio of retrieved CH₄ excess column, divided by the posterior precision error, over a simulated scene in the Durango, Colorado, area for plumes of various emission rates. The spatial resolution is 30 m, the instrument FWHM is 1.5 nm, and the polynomial degree in the retrieval is 25. The realistic surface distribution is based on LANDSAT data. Ratios higher than 4 imply a probability, of lower than 1 in 15 000, that the pixel enhancement happens by random noise.

The retrieved scenes for the CH₄ emission rates, from 50 to 2000 kg h⁻¹, are shown in Fig. 18, and the corresponding scenes showing the ratio of retrieved methane concentration and precision error in each pixel are given in Fig. 19. The ratio, n, of the pixel enhancement to the precision error represents a n−σ probability event that this pixel enhancement would have randomly happened, purely due to noise, in the absence of a true CH₄. A ratio value above 4 would imply that there is only a probability of lower than 1 in 15 000 that this would happen by chance due to random noise. Thus, the ratio of 4 can be a simple and useful metric to imply where the actual CH₄ enhancement pixels are. Based on this metric, we found that CH₄ source detection can still be possible for plumes with emission rates as low as approximately 50–100 kg h⁻¹.

4 Conclusions

We built an end-to-end modeling framework that can simulate radiances from reflected sunlight through methane plumes over a variety of surfaces. In this study, we simulated a realistic 3D CH₄ concentration field from a point source, using an LES, and varied the underlying surfaces where the emission occurs using a comprehensive surface albedo database from the ECOSTRESS spectral library consisting of over 2000 surface types. The observed radiances and their NESR for various instrument configurations are modeled directly as a function of incoming radiance and instrument parameters, such as the FWHM of the line shape function and the integration time, without having to prescribe the SNR a priori. Based on the modeled radiance, we applied the IMAP-DOAS algorithm to retrieve methane column enhancements.

We compared the tradeoff between the two most frequently used fitting windows for CH₄ in the 1.6 and 2.3 µm ranges. Our analysis has shown that, despite a higher solar radiance near 1.6 µm, the stronger absorption feature of CH₄ near 2.3 µm leads to a consistently lower precision error for the 2.3 µm fitting band. The rare occasion of a 1.6 µm fitting window outperforming the 2.3 µm retrieval band happens only when both FWHM is lower than 0.2 nm and, simultaneously, a surface in consideration has a much higher albedo near 1.6 compared to 2.3 µm. For the purpose of building an instrument to detect methane emissions accurately at sufficiently fine spatial resolutions across most global surfaces, we believe that the 2.3 µm band can perform better than the 1.6 µm band and should be prioritized in most scenarios. We primarily considered the fitting window in the 2120–2395 nm range to study the impact of instrument parameters and retrieval choices on the retrieval bias and precision error.

To highlight the impact of surface interferences, the number of polynomial degrees is varied in the IMAP-DOAS retrieval experiments. This framework allows us to derive the corresponding precision error and bias when different sets of instrument parameters and the number of polynomial degrees are used in the retrieval of CH₄ column concentration. Our analysis shows that the number of polynomial degrees used to represent surface spectral variations in the retrieval algorithm has a significant impact on the bias of the retrieved methane columns, causing a positive bias as large as 0.8 mole m⁻² for a retrieval with a polynomial degree of 5 compared to 0.2 and only 0.02 mole m⁻² for degrees of 11 and 25, respectively, across the majority of surfaces. Using a higher polynomial degree, however, is found to simultaneously increase the precision error for methane retrieval, as this relaxes the constraints on the possible methane absorption contribution in spectral variations of the observed radiance. This is particularly evident at a FWHM of greater than 2 nm. Thus, using an instrument with a lower FWHM, such as 1.5 nm, will allow for a high number of polynomial terms to be used while inducing a smaller deterioration in precision error. For example, we found that an instrument with a FWHM of 1.5 nm and the exposure time of 20 ms can achieve a precision error of less than 0.007 mole m⁻² (or less than 1.0 % of the total column in the atmosphere) over a typical construction concrete surface with an albedo of 0.35, even when a polynomial degree of 25 is used.

Having low bias in the retrieval is integral to removing correlated surface features in the retrieval enhancement map. These surface features from retrieval errors likely appear as false positives and, subsequently, cause significant impacts on the detection and quantification of true CH₄ sources. To demonstrate the significance of surface interferences, we used a realistic surface distribution over a well pad from the Durango area, Colorado, as a background surface with an LES methane plume of 200 kg h⁻¹ to create a synthetic emission in a real-world environment. Our results illustrated that, when an instrument with a high FWHM (such as 5 nm) is used with a low polynomial degree (such as 5), a large retrieval bias appears broadly across the 2D scene. These interferences occur severely over surfaces such as oil shale and sandstone, resulting in difficulties in clearly distinguishing a true plume from areas of systematic biases. Nevertheless, by using a lower FWHM value, such as 1.5 nm, and the polynomial degree of at least 25, we have illustrated the ability to obtain a low retrieval bias across the entire scene and to effectively differentiate the source location from the background. We also repeated the 2D retrieval analysis for 30 m spatial resolution, by averaging the radiance per unit area from 5 to 30 m², and adjusted the exposure time to 40 ms, which could be achievable for future satellites. Again, our results have shown that using a FWHM of 1.5 nm and a polynomial degree of 25 plays a crucial role in resolving surface features and removing false positives, ultimately giving us the ability to distinguish the true emission location. In the absence of bias, the ratio of the retrieved column enhancement and the retrieval precision error in our retrieved 2D scenes indicates that it might be possible to detect a CH₄ emission source with a flux rate as low as 50–100 kg h⁻¹.

This study highlights the effect of changing the instrument FWHM, exposure times, and the polynomial degree on minimizing retrieval errors. The FWHM and exposure time are intrinsic to how a spectrometer is designed, as opposed to describing a spatial resolution, which depends on external factors such as the viewing geometry and the speed of a remote sensing platform. Further studies will be required to determine how these variables are implemented into an observing system that can achieve specific spatial resolutions of interest. Additional considerations, such as the device saturation constraint, will also influence the ultimate achievable exposure time for a newly designed instrument. Our end-to-end simulator, which includes an instrument model and retrieval, can be generalized to study the performance of future instruments with specific engineering requirements.

The findings in this study can inform future satellite instrument designs and the retrieval algorithm in order to have robust observations capable of separating real plumes from surface interference. Reducing both bias and precision error can have a profound benefit, both for manual analysis by humans and for an automated model plume detection, such as an artificial neural network (ANN) approach. This will enable the analytic chain to have a higher accuracy and level of confidence in detecting and quantifying more subtle methane sources from observed scenes across large geographical areas.

The modeling framework in this work could also be generalized to improve the detection and quantification of CO₂, with minor modifications on a different fitting window and a different magnitude of flux rate in the emission simulations.