The GHGSat-D imaging spectrometer

The demonstration satellite GHGSat-D, or “Claire”, launched on 21 June 2016, is the first in a planned constellation of small satellites designed and operated by GHGSat, Inc. to measure greenhouse gas emissions at the facility scale from space. Its instrument measures methane concentrations by collecting and spectrally decomposing solar backscattered radiation in the shortwave infrared using a compact fixed-cavity Fabry–Pérot imaging spectrometer. The effective spatial resolution of 50×50 m² over targeted 12×12 km² scenes is unprecedented for a space-based gas-sensing spectrometer. Here we report on the instrument design and forward model and retrieval procedure, and we present several examples of retrieved methane emissions observed over industrial facilities. We discuss the sources of error limiting the performance of GHGSat-D and identify improvements for our follow-on satellites. Claire's mission has proven that small satellites can be used to identify and quantify methane emissions from industrial facilities, enabling operators to take prompt corrective action.

1 Introduction

GHGSat is a Canadian company incorporated in 2011 with the goal to provide a precise, scalable, and economical method of measuring greenhouse gas (GHG) emissions from industrial facilities worldwide. The first GHGSat satellite, GHGSat-D or “Claire”, is a demonstration small satellite that measures surface-level methane emission plumes with high spatial resolution for facility-scale attribution. The instrument on GHGSat-D is a wide-angle fixed-cavity Fabry–Pérot (F–P) imaging spectrometer able to resolve methane (CH₄) absorption lines in the shortwave infrared (SWIR; Sloan et al., 2016), where water (H₂O) and carbon dioxide (CO₂) absorption lines are also present. The distinguishing features of GHGSat-D compared with other GHG remote sensing missions are its (1) combination of high spatial resolution (∼50 m) and fine spectral resolution (∼0.1 nm) and (2) compact package (Table 1). To this day GHGSat-D remains the only gas-sensing satellite with such high spectral and spatial resolution. This combination enables unique capabilities for imaging and quantifying emission plumes with unambiguous attribution at the facility scale. Delivering this performance in a compact package enables low-cost launch and operation of the satellite so that a constellation can be used for high-density coverage. To date, GHGSat-D has performed over 5000 observations of commercial facilities in oil–gas, power generation, coal mining, waste management, and agriculture sectors around the world.

Table 1 Satellite parameters for GHGSat-D.

1.1 Methane monitoring for industry

Methane is a potent greenhouse gas whose concentration has increased from 720 to 1800 ppb since pre-industrial times (Hartmann et al., 2013) and is responsible for a radiative forcing of 0.97 W m⁻², second only to CO₂ (Myhre et al., 2013). Because methane has a relatively short atmospheric lifetime of ∼10 years compared to CO₂, actions taken that reduce methane emissions will have a significant effect on the near-term warming rate. Anthropogenic methane emissions originate from a very large number of point sources including oil–gas facilities, coal mines, landfills, wastewater treatment plants, and confined livestock operations. Studies have shown that a relatively small fraction of sources are responsible for the majority of methane point source emissions, with 60 %–90 % of overall emissions coming from emitters with flux Q₀≥100 kg h⁻¹ (Brandt et al., 2016; Duren et al., 2019). A fleet of satellites with a detection threshold at or below Q₀ could lead to efficient emission abatement – giving operators the opportunity to take corrective action, often at no net cost (IEA, 2018).

1.2 Space-based GHG monitoring

Facility-scale greenhouse gas emissions can be monitored in a variety of ways including stationary ground-based systems (Chen et al., 2016; Robinson et al., 2011), mobile ground-based measurements (Yacovitch et al., 2015), and aircraft observations (Conley et al., 2016; Duren et al., 2019; Sherwin et al., 2020; Wolff et al., 2020). Satellites are an attractive complementary observation platform since they have global coverage, employ the same measurement method for any observation site in the world, and can repeatedly revisit any facility in the world. In recent years, building on the pioneering work of the global mapping missions SCIAMACHY, GOSAT, OCO2, and S5P/TROPOMI (Burrows et al., 1995; Hamazaki et al., 2005; Veefkind et al., 2012), satellites have emerged as a candidate technology to measure individual anthropogenic emission sources (Nassar et al., 2017; Pandey et al., 2019). For example, the TROPOMI instrument on board the Sentinel-5P satellite provides daily global coverage with a spatial resolution of several kilometres and can be used to “tip and cue” GHGSat satellites to locate facility-scale emission sources. TROPOMI has also been used to detect and quantify very strong emitters (≳10 t h⁻¹), including an industrial “blow-out” event (Pandey et al., 2019). In another example, the Hyperion imaging spectrometer with 30 m spatial resolution but coarse spectral resolution (10 nm) was used to quantify the Aliso Canyon methane blowout event (Thompson et al., 2016). However, until the launch of GHGSat-D, instruments with fine spectral resolution had spatial resolutions limited to the kilometre scale and above. Current and future hyperspectral imagers PRISMA (Loizzo et al., 2018), EnMAP (Guanter et al., 2015), and EMIT (Green et al., 2018) have spatial resolution similar to GHGSat (30–60 m) but much coarser spectral resolution (7–10 nm). Recently, the discovery of a number of large methane leaks in an oil–gas production area was reported in Varon et al. (2019). GHGSat-D's high spatial resolution enabled the attribution of the leaks to specific pieces of equipment within the industrial site, the locations of which were promptly communicated to the site operator. GHGSat-D has also combined multiple single-pass measurements to quantify time-averaged methane emission rates from coal mine vents (Varon et al., 2020).

Whereas satellite instruments can detect the presence of emissions by measuring the column density enhancement relative to background, emission quantification requires a method to model the local transport using meteorological information (Jongaramrungruang et al., 2019; Nassar et al., 2017; Pandey et al., 2019; Varon et al., 2018). Locally measured meteorological information is usually not available, so wind speed and direction data must be inferred from the plume observations and/or drawn from meteorological databases like the NASA Goddard Earth Observation System Fast Processing (GEOS-FP) reanalysis product (Molod et al., 2012), or similar sources.

1.3 Advantage of high spatial resolution

The primary motivation for a high-spatial-resolution methane measurement is to identify the industrial facility, or even the piece of infrastructure within the facility, responsible for the detected emissions, to provide actionable information to the operator. High-spatial-resolution measurements offer additional advantages: (1) the ability to identify most types of clouds within a scene mitigates the need for complex cloud screening algorithms; (2) the ability to image the shape of the emission plume can help infer wind direction and speed (Jongaramrungruang et al., 2019); and (3) the measured column density enhancement ΔΩ in a square ground cell with side length L generally increases for decreasing L.

Figure 1 Plume geometry scenarios for illustrating the dependence of detected enhancement ΔΩ on pixel size L.

This last point deserves elaboration. The impact of pixel size on the ability to monitor CO₂ emissions was studied in Hill and Nassar (2019). A useful heuristic introduced by Jacob et al. (2016) asserts that $\mathrm{\Delta }\mathrm{\Omega }\propto \mathrm{1}/L$. However, this scaling does not hold for all plume geometries. Three scenarios are illustrated in Fig. 1, where b is the scenario as considered by Jacob et al. (2016). Here, for constant emission rate and wind speed, the plume enhancement is contained within the pixel boundary along one axis as L is reduced but “escapes” the pixel boundary along the other axis. Therefore, the amount of excess gas in a pixel of size L scales linearly with L, giving the $\mathrm{\Delta }\mathrm{\Omega }\propto \mathrm{1}/L$ scaling relationship. In Fig. 1a, on the other hand, the emitted gas remains entirely within the boundary of the smaller pixel L₁, leading to an excess mass invariant with L, and $\mathrm{\Delta }\mathrm{\Omega }\propto \mathrm{1}/L^{\mathrm{2}}$. The scenario in Fig. 1a could occur for very low wind speed and/or a transient event where emissions begin just before the satellite observation. Local wind eddies could also produce point-like enhancements. While geometry a is likely less common than b, it is notable since it gives a scaling more favourable than b to smaller pixels. Finally, scenario c considers a quasi-uniform density field (plume extent greater than L₂) – in this case ΔΩ is approximately invariant with L, and hence there is no advantage to smaller pixel sizes. Geometry c might describe widely dispersed plume enhancements far downwind of the source location or the enhancement from an area emission source. This simplified treatment is based on averaging the excess density within a cell, which neglects the nonlinear character of Beer's law when averaging radiance spectra. When this effect is included, it suppresses ΔΩ relative to the simple scaling relationships above when the peak enhancement within a cell significantly exceeds the mean. This leads to an additional, if modest, advantage to smaller pixel sizes.

Space-based methane sensors with kilometre-scale spatial resolutions like TROPOMI have focused on achieving high-precision column measurements (<0.95 % single measurement) coupled with high absolute accuracy (<0.4 % single measurement; Hu et al., 2018). In contrast, a methane sensor with high spatial resolution designed to detect local enhancements above background can have relaxed absolute accuracy requirements. For example, a 5 % bias in the absolute methane column density – on the upper end of what is estimated from neglecting scattering due to aerosols, for instance (Aben et al., 2007; Butz et al., 2009; Houweling et al., 2005) – would lead to a 5 % bias in the local methane enhancement and corresponding emission estimate. This error level is much smaller than the total error in a typical emission estimate (Varon et al., 2019), which is usually dominated by wind uncertainty. Furthermore, any additive scene-wide bias in the absolute methane column density has negligible impact on the ability of the imaging spectrometer to distinguish the local methane enhancement from background and to quantify its magnitude. This holds true even if the additive bias is not consistent between observations – the key requirement is that the retrieval errors not have sharp gradients within the scene.

2 Instrument overview

The GHGSat-D satellite (Table 1) uses the NEMO platform from the University of Toronto Institute for Aerospace Studies (UTIAS) Space Flight Laboratory. The NEMO platform can host a payload with a mass of up to 6 kg and provide up to 45 W of power. GHGSat-D was launched on Indian Polar Satellite Launch Vehicle C-34 on 21 June 2016 at 23:56 EDT. It was fully commissioned in July 2016 with a 3–5-year expected lifetime. GHGSat-D remains operational at the time of writing.

The patented GHGSat-D instrument operates in the shortwave infrared (SWIR) between 1630–1675 nm where methane, CO₂, and water vapour absorption lines are present. Critical instrument parameters are listed in Table 2. The spectral resolution is determined by the F–P gap spacing, its uniformity, and coating reflectivities. The F–P coating reflectivity is a critical parameter that leads to a trade-off between the per-pixel signal and spectral resolution.

Table 2 Instrument parameters.

Figure 2 (a) The GHGSat-D spacecraft with the imaging spectrometer on board. (b) The mounted Fabry–Pérot interferometer. (c) Schematic of the unfolded optical system with the (i) OSF, (ii) F–P, and (iii) detector identified. The red, blue, and green rays originate from different ground locations.

2.1 Optical system

The GHGSat-D optical system (Fig. 2) is composed of three lens assemblies with focal lengths f_t1,f_t2, and f_im: the first two lenses, in confocal arrangement, constitute the telescope, and the last lens, the imaging assembly, forms a two-dimensional image of the ground on the detector. The F–P is placed in the Fourier plane of the optical system, between f_t2 and f_im. An order sorting filter (OSF), placed between the lenses that make up the f_t2 assembly, defines the spectral bandpass region. The 1630–1675 nm bandpass is chosen to allow transmission of the methane R and Q branches. The choice of focal lengths must balance several considerations simultaneously. First, the angular magnification ratio of the telescope $|f_{\mathrm{t}\mathrm{1}}/f_{\mathrm{t}\mathrm{2}}|$ must be consistent with the mechanical constraints on input aperture and F–P size and the spectroscopic constraints that dictate the desired range of incident ray angles on the F–P. Second, the effective focal length $f=f_{\text{im}}|f_{\mathrm{t}\mathrm{1}}/f_{\mathrm{t}\mathrm{2}}|$ is constrained by the choice of spatial resolution. The ground sampling distance (GSD) of GHGSat-D was chosen to be L≈25 m: small enough to resolve facility features, yet large enough to image a ∼5 km sized facility and the full extent of an emission plume. This constrains the effective focal length f of the optical system given the camera pixel size a and orbiting altitude h: $f=ah/L\approx \mathrm{500}$ mm. Optical aberrations in the imaging system limit the spatial resolution of GHGSat-D to approximately 50 m.

For a given image, a polychromatic light ray originating from a specific ground location enters the optical system through the input aperture at some angle pair (ψ,ϕ), where ψ is the small elevation angle and ϕ is the azimuthal angle. In the paraxial approximation, the light ray emerges from the telescope with angle pair $\left(\left|\frac{f_{\mathrm{t}\mathrm{1}}}{f_{\mathrm{t}\mathrm{2}}}\right|\mathit{\psi },\mathit{\varphi }\right)$ and is incident on the F–P. The imaging assembly then focuses the light ray to detector pixel

$$\begin{array}{}\text{(1)}& \begin{array}{rl}(i,j)=(& i_{\mathrm{0}}+b\left(\left|{\displaystyle \frac{f_{\mathrm{t}\mathrm{1}}}{f_{\mathrm{t}\mathrm{2}}}}\right|\mathit{\psi }\right)\cos \left(\mathit{\varphi }\right),\\ & j_{\mathrm{0}}+b\left(\left|{\displaystyle \frac{f_{\mathrm{t}\mathrm{1}}}{f_{\mathrm{t}\mathrm{2}}}}\right|\mathit{\psi }\right)\sin \left(\mathit{\varphi }\right)),\end{array}\end{array}$$

where b is a proportionality constant relating angle to pixel radius. The optical axis intercepts the 640×512 temperature-controlled InGaAs detector array at pixel (i₀,j₀).

2.2 Fabry–Pérot

The GHGSat-D F–P element consists of two optical flats mounted within a mechanical enclosure. The two optical surfaces are positioned such that the inner surfaces, with reflectivity R and spaced a distance d apart, form an optical cavity. Light with wavelength λ and incident angle $\mathit{\theta }=\left|\frac{f_{\mathrm{t}\mathrm{1}}}{f_{\mathrm{t}\mathrm{2}}}\right|\mathit{\psi }$ with respect to the F–P surface normal is transmitted according to

$$\begin{array}{}\text{(2)}& T_{\text{FP}}\left(\mathit{\theta },\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{\mathrm{2}F}{\mathit{\pi }}\right)}^{\mathrm{2}}{\sin }^{\mathrm{2}}\left(\frac{\mathrm{2}\mathit{\pi }nd\cos \left(\mathit{\theta }\right)}{\mathit{\lambda }}\right)}},\end{array}$$

where n is the index of refraction of the medium within the optical cavity and $F=\mathit{\pi }\sqrt{R}/(\mathrm{1}-R)$ is the reflectivity finesse. For each value of θ, the transmission spectrum of the F–P is a series of peaks that are spaced in wavelength by the free spectral range $\text{FSR}={\mathit{\lambda }}^{\mathrm{2}}/\left(\mathrm{2}d\right)$ with a spectral width characterized by the full width at half maximum $\text{FWHM}=\text{FSR}/F$. Because the F–P accepts a continuum of θ values, it samples a continuum of wavelengths within the passband.

In applications using F–P-based spectrometers, the F–P gap spacing d is often dynamically scanned during the measurement (Reay et al., 1974). In contrast, GHGSat-D uses a fixed gap spacing and exploits the angular dependence of the mth F–P transmission mode's spectral position ${\mathit{\lambda }}_m=\mathrm{2}\mathit{\pi }nd\cos \left(\mathit{\theta }\right)/m$ to measure the spectrum of the incident light. This approach simplifies the mechanical design and makes it much easier to meet stringent stability requirements. The FSR is chosen so that there are enough F–P transmission modes within the spectral bandpass (usually three at a given θ) to sample the entirety of the spectra within the imaging field of view. The F–P is temperature-controlled to keep thermal mechanical drift to a minimum.

3 Measurement concept

The GHGSat-D spectrometer is based on a wide-angle Fabry–Pérot (WAF–P) imaging concept (Sloan et al., 2016). A programmable number of closely overlapping two-dimensional images are taken (typically 200 in nominal operations) in which the atmospheric absorption spectrum is “imprinted” on the images in the form of spectral rings due to the angle-dependent Fabry–Pérot transmission spectrum. During the observation sequence, the ground target traverses the field of view, sampling the full extent of the spectral information contained in the images. The instrument operates in “target” mode in which the satellite attitude is adjusted to the keep the facility of interest in the field of view for much longer than it would if operated in nadir (downward-pointing) mode, thereby increasing the available integration time and hence signal-to-noise ratio (SNR).

Figure 3 The plots in (a, b) illustrate how the instrument signal in (c) is produced. In (a), the F–P and OSF transmission spectra are shown alongside the normalized top-of-atmosphere spectral radiance (TOASR) for the case of light rays incident on the F–P at normal incidence (i.e. $\left(r,\mathit{\theta }\right)=\mathrm{0}$). In (b), the location of the F–P transmission peaks is shown as a function of radius (gold lines) overlaid on the normalized TOASR (grayscale background image). The horizontal dark bands at the top and bottom of (b) illustrate wavelengths where the OSF transmission is reduced to zero. By sampling a continuum of incident angles θ, the transmitted F–P transmission peaks measure a continuum of wavelengths within the passband. The instrument signal (c) results from integrating the multiplied signals in (b) along the vertical (wavelength) axis. The change in instrument signal with respect to a change in the methane vertical column density is plotted in (d). The TOASR, instrument signal, and change in signal are calculated assuming a target elevation at sea level, a solar zenith angle of 40^∘, and vertical column densities of 0.68 (methane), 160 (CO₂), and 830 mol m⁻² (water vapour).

Figure 3 illustrates how the F–P samples the backscattered solar radiance spectra to generate spectral rings in the image. Multiple F–P transmission modes are allowed through the OSF bandpass. Because the F–P transmission function depends only on θ, it is circularly symmetric and so can be expressed as a function of radius $r=\sqrt{{\left(i-i_{\mathrm{0}}\right)}^{\mathrm{2}}+{\left(j-j_{\mathrm{0}}\right)}^{\mathrm{2}}}$. The r=0 F–P transmission spectrum is shown in Fig. 3a alongside the OSF transmission function and the normalized backscattered top-of-atmosphere spectral radiance (TOASR). For larger radii – and thus larger θ – the F–P spectrum shifts to lower wavelengths with a cos (θ) dependence, allowing us to sample different regions of the TOASR. Figure 3b shows the location (in wavelength space) of the F–P transmission peaks as a function of radius overlaid on the normalized TOASR. The instrument signal is shown in Fig. 3c and the sensitivity of this signal to a change in methane vertical column density is shown in Fig. 3d. At each radius, the signal on the detector array is the wavelength integral of the TOASR multiplied by the F–P and OSF transmission spectra, i.e. the result of integrating Fig. 3b along the vertical axis. A mathematical description of the forward model is given in Sect. 4.1.

In order to measure the spectrum of solar radiation backscattered from a specific ground cell, the location of the ground cell in each image must be known. This is done with an image co-registration algorithm. We then construct a spectrum for each ground cell along the image frame axis by recording the measured signal as a function of the ground cell's radial position with respect to the spectral ring center. Figure 4a–d show an observation where the location of an example ground cell has been tracked in each frame. The constructed spectrum is shown in Fig. 4e. The colour of the data point represents the image frame from which the data was obtained. We construct approximately 200 000 of these spectra in order to retrieve the methane column density for each ground cell within the retrieval domain.

Figure 4 The panels (a, b, c, d) show a selection of frames from an observation over the Lom Pangar hydroelectric reservoir in Cameroon taken on 20 April 2017 with an example ground location (denoted by an orange “x”) tracked in each frame. The image axes are in pixels, with each pixel representing a 24×24 m² area on the ground. The plot in (e) shows the signal (circles) from the example ground location as a function of the image frame (circle colour) and radius from spectral ring center (horizontal axis). The forward model (black line) is plotted alongside the signal data, and residuals between model and data are shown in (f).

4 Retrieval method

The goal of the retrieval algorithm is to estimate the instrument and atmospheric state vector x from a measurement vector y. This is done by constructing a combined forward model F(x) of the instrument and atmosphere and making the association

$$\begin{array}{}\text{(3)}& \mathit{y}=\mathit{F}\left(\mathit{x}\right)+{\mathit{ϵ}}_y+{\mathit{ϵ}}_{\mathrm{F}},\end{array}$$

where ϵ_y represents the measurement error and ϵ_F represents error in the forward model. A retrieval of x requires that we have accurate knowledge of both the forward model and the errors in the measurement system. Because F(x) is a nonlinear function of x, we must solve for the state vector iteratively. This requires knowledge of the Jacobian $\mathbf{K}\left(\mathit{x}\right)=\frac{\partial \mathit{F}\left(\mathit{x}\right)}{\partial \mathit{x}}$ to weight the state vector step Δxⁱ taken during the ith iteration. In this section, we describe the instrument and atmospheric forward model and outline the inversion procedure used to estimate x.

4.1 Forward model

The forward model represents our best knowledge of the instrument and atmosphere, with approximations used to evaluate the model more efficiently when performing retrievals. The camera signal F_i,j at detector pixel (i,j) in photocurrent units [${\mathrm{e}}^{-}{\mathrm{s}}^{-\mathrm{1}}$] is given by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{F}_{i,j}\left(\mathit{x}\right)& =\int L\left(\mathit{x},\mathit{\lambda }\right)\cdot C\left(\mathit{\lambda }\right)\cdot \text{QE}\left(\mathit{\lambda }\right)\cdot T_{\text{OSF}}\left(\mathit{\lambda }\right)\\ & \cdot T_{\text{FP}}\left(\mathit{\theta },\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda },\end{array}\end{array}$$

where L(x,λ) is the spectral radiance as a function of the state parameterx and wavelength λ, C(λ) is the radiometric conversion factor that converts spectral radiance to the number of photons on a pixel per unit time, QE(λ) is the quantum efficiency with which the camera converts a photon to electric charge, T_OSF(λ) is the transmission of the order-sorting filter that defines the spectral bandpass region, and T_FP(θ,λ) is the F–P transmission function defined in Eq. (2). The camera signal in Eq. (4) is plotted as a function of radius in Fig. 3c.

4.2 Atmospheric model

The spectral radiance L(x,λ) is calculated from the spectral irradiance I(λ) assuming Lambertian surface reflectance:

$$\begin{array}{}\text{(5)}& L\left(\mathit{x},\mathit{\lambda }\right)={\displaystyle \frac{a\left(\mathit{\lambda }\right)\cdot \cos \left({\mathit{\theta }}_{\text{sza}}\right)}{\mathit{\pi }{R}_{\text{E–S}}^{\mathrm{2}}}}I\left({\mathit{x}}^{\prime },\mathit{\lambda }\right),\end{array}$$

where a(λ) is the spectrally dependent surface albedo, θ_sza is the solar zenith angle, R_E–S is the relative Earth–Sun distance, x^′ is the state parameter vector without the albedo, and the spectral irradiance is the solution to a simplified radiative transfer equation where thermal emission, aerosol, and molecular scattering have been neglected (Chandrasekhar, 1950):

$$\begin{array}{}\text{(6)}& \mathit{\mu }{\displaystyle \frac{\partial I\left({\mathit{x}}^{\prime },\mathit{\lambda }\right)}{\partial z}}=-{\mathit{\alpha }}_{\text{abs}}I\left({\mathit{x}}^{\prime },\mathit{\lambda }\right).\end{array}$$

This equation is integrated along the downwelling and upwelling light path. Here μ=cos (θ); θ is the angle at which the light travels through the atmosphere with respect to the Earth's surface normal; z is the altitude; α_abs is the pressure-, temperature-, wavelength-, and species-dependent absorption coefficient calculated using the HITRAN API (Kochanov et al., 2016); and the solar irradiance is introduced through a boundary condition and generated from the AER solar irradiance model (Clough et al., 2005). We integrate the radiative transfer equation discretely assuming 100 atmospheric layers that are evenly spaced in pressure. The pressure, temperature, and molecular mixing ratio profiles are taken from the US standard data set (U.S. COESA, 1976). The target elevation is determined from the SRTM 30 m product (Farr et al., 2007).

We justify excluding thermal emission and molecular scattering from the atmospheric model because both are small effects for the wavelengths within our spectral bandpass region. Previous studies of simulated carbon dioxide retrievals using only the 1563–1585 nm band have found that neglecting aerosol and molecular scattering can lead to a few percent error, depending on the surface albedo and aerosol optical depth (Aben et al., 2007). This error can be either positive or negative, depending on whether the presence of aerosols – in combination with the surface albedo – leads to an increase or decrease in the average optical path length. An analysis of SCIAMACHY retrievals in this same spectral region found that errors decrease for aerosol vertical distributions that are narrow and closer to the Earth's surface (Houweling et al., 2005). As mentioned in Sect. 1.3, GHGSat retrievals are primarily intended to measure local plume enhancements. Therefore, we are especially concerned with any unmodelled effects with spatial structure on the length scales of emission plumes. This could potentially include aerosol scattering, such as aerosols that might conceivably be co-emitted with methane plumes. However, since the presence of these aerosol plumes would be much closer to the surface and narrower in vertical profile than the aerosol profiles retrieved in Aben et al. (2007) and Houweling et al. (2005), we expect that errors arising from neglecting scattering should be small compared with other sources of measurement error, similar to what is assumed in AVIRIS airborne methane retrievals (Thorpe et al., 2014). We also note a recent study in which the effect of neglecting aerosols in an AVIRIS airborne methane plume retrieval was shown to be <5 % for representative methane plume enhancements and a significant aerosol optical depth (Huang et al., 2020).

4.3 Inversion procedure

For any ground cell (p,q) in a reference frame, we can compare the observation data vector ${\mathit{y}}^{\left(pq\right)}=\left\{y_{i_{\mathrm{1}},j_{\mathrm{1}}}^{\left(pq\right)},y_{i_{\mathrm{2}},j_{\mathrm{2}}}^{\left(pq\right)},\mathrm{\dots },y_{i_k,j_k}^{\left(pq\right)}\right\}$ to the forward model vector $\mathit{F}\left({\mathit{x}}^{\left(pq\right)}\right)=\left\{F_{i_{\mathrm{1}},j_{\mathrm{1}}}\left({\mathit{x}}^{\left(pq\right)}\right),F_{i_{\mathrm{2}},j_{\mathrm{2}}}\left({\mathit{x}}^{\left(pq\right)}\right),\mathrm{\dots },F_{i_k,j_k}\left({\mathit{x}}^{\left(pq\right)}\right)\right\}$ and infer the state vector x^(pq) using a variant of standard inverse methods described below. The subscripts refer to the pixel indices for this ground cell within the respective frames of the image sequence from 1 to k. The reference frame coordinate system can then be georeferenced with the appropriate rotational and scale transformation.

We use optimal estimation (Rodgers, 2000) to infer the posterior distribution of the state vector given the observation data, an error model, and a prior distribution for the state vector. Assuming a Gaussian form for the measurement and prior probability density functions, maximizing the joint probability density function amounts to minimizing the cost function:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\chi }}^{\mathrm{2}}\left(\mathit{x}\right)& ={\left(\mathit{y}-\mathit{F}\left(\mathit{x}\right)\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{o}}^{-\mathrm{1}}\left(\mathit{y}-\mathit{F}\left(\mathit{x}\right)\right)\\ & +{\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right),\end{array}\end{array}$$

where S_o is the observation error covariance matrix, S_a is the prior covariance matrix, and x_a is the prior state vector. The Gauss–Newton procedure for minimizing the cost function requires that we update the state vector at each iteration by a step:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{x}}^{i+\mathrm{1}}& ={\left({\mathbf{K}}^{\mathrm{T}}\left({\mathit{x}}^i\right){\mathbf{S}}_{\mathrm{o}}^{-\mathrm{1}}\mathbf{K}\left({\mathit{x}}^i\right)+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\\ & \cdot \left({\mathbf{K}}^{\mathrm{T}}\left({\mathit{x}}^i\right){\mathbf{S}}_{\mathrm{o}}^{-\mathrm{1}}\left(\mathit{y}-\mathit{F}\left({\mathit{x}}^i\right)\right)+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)\right),\end{array}\end{array}$$

where K(xⁱ) is the Jacobian of the forward model evaluated at xⁱ.

At each iteration of the Gauss–Newton procedure, the forward model and Jacobian must be evaluated. This is computationally expensive for a single cell, and evaluating it for the ∼200 000 ground cells in our field of view is impractical. Instead, we use a two-step procedure: (1) a scene-wide average retrieval using the full forward model to estimate the scene-wide average state vector $\widehat{\mathit{x}}$ and (2) a per-cell retrieval done using a linearized forward model (LFM) evaluated at the linearization point $\widehat{\mathit{x}}$. A full retrieval takes approximately 30 min using 48 cores, with approximately 10–20 iterations required to meet the convergence criteria for the scene-wide average retrieval and approximately 10 iterations required for the per-cell retrieval. The separate retrieval steps are described in the following sections.

Table 3 Scene-wide average atmospheric state vector elements.

4.4 Scene-wide retrieval

The goal of the scene-wide retrieval is two-fold: to retrieve scene-wide averaged surface and atmospheric parameters such as albedo and molecular column density and to retrieve the F–P gap spacing. Even though the F–P is thermally stabilized, residual drift in the F–P gap spacing can occur between observations. Because the F–P gap spacing directly affects the signal on each detector pixel, we retrieve d for each observation. The scene-wide retrieval uses the full instrument model F_i,j(x) from Eq. (4) in the optimal estimation procedure. The data vector y in the scene-wide retrieval is the radial average of the average of all image frames. The prior x_a uses same-scene information from Landsat 8 for the albedo parameter (Roy et al., 2014) and closest-in-time methane, CO₂, and water vapour values from AIRS (Chahine et al., 2006). The result of the scene-wide retrieval is the state parameter estimate $\widehat{\mathit{x}}$ which includes not only retrieved instrument parameters but also the scene-wide averaged methane, CO₂, and water vapour column densities as well (see Table 3). The elevation is assumed to be flat within the scene.

4.5 Spatially resolved column retrieval

A linearized forward model (LFM) is constructed at the linearization point $\widehat{\mathit{x}}$:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{F}}^{\text{LFM}}\left({\mathit{x}}^{\left(pq\right)}\right)& ={\mathit{x}}^{\left(pq\right)}\left(\mathrm{1}\right)({\mathbf{K}}_{\mathrm{1}}\left(\widehat{\mathit{x}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\widehat{\mathit{x}}\left(\mathrm{1}\right)}}\sum _{l=\mathrm{2}}^n{\mathbf{K}}_l\left(\widehat{\mathit{x}}\right)\left({\mathit{x}}^{\left(pq\right)}\left(l\right)-\widehat{\mathit{x}}\left(l\right)\right))\\ & +\sum _{l=n}^{n+m}{\mathbf{K}}_l\left(\widehat{\mathit{x}}\right)\left({\mathit{x}}^{\left(pq\right)}\left(l\right)-\widehat{\mathit{x}}\left(l\right)\right),\end{array}\end{array}$$

where the first element of the state vector is taken to be the surface reflectance. The methane, CO₂, and water vapour column densities are all retrieved for each ground cell (p,q). Here, we set the molecular column density prior variances to be very large such that the retrieved parameters are almost entirely determined by the data. There are two terms in the LFM: one with n terms that includes state parameters whose Jacobians scale with the surface reflectance (for example the molecular column densities) and another with m terms for state parameters whose Jacobians are not scaled by the reflectance. In particular, the molecular components of the Jacobian are calculated using only a column density enhancement in the lowest atmospheric layer that extends from the surface to approximately 100 m altitude. The primary advantage of using the LFM is that we only compute the forward model and Jacobians once at the beginning of the per-cell retrieval. A disadvantage is that for the parameters that are nonlinear in the forward model, a retrieval using the LFM will introduce systematic biases for deviations far from the linearization point. For the particular case of molecular column densities, this leads to an underestimation that is corrected in post-processing using a non-linear correction function determined from a comparison of F^LFM(x) with F(x) at the appropriate linearization point (Varon et al., 2019).

4.6 Data processing

The inversion procedure is performed after data downlinked from the satellite have been processed to the following data levels (Kobler et al., 1995).

• Level 0. Information received from the satellite is removed of all communication-related artefacts. Telemetry data are parsed and stored separately from the image observation data.

• Level 1A. Telemetry data are processed to provide instrument position, orientation, solar zenith and observation angles, etc.

• Level 1B. The recorded digital signal values in the image observation are converted to photocurrent units [${\mathrm{e}}^{-}{\mathrm{s}}^{-\mathrm{1}}$] after correction for pixel offset and dark current. Corrections are also applied to mitigate optical ghosting and detect and flag dead, hot, or otherwise misbehaving pixels.

Because GHGSat-D does not contain an on-board calibration unit, changes in pixel offset signal and dark current are measured by taking frequent observations over a dark ocean scene. The observation data vector y that the optimal estimation procedure is performed on is then the Level 1B data that have been “dark” corrected by the most recent ocean measurements. The optimal estimation procedure generates the spatially resolved state vector elements x^(pq). These are then used to generate the Level 2 product, consisting of georeferenced arrays of methane abundances and their uncertainties.

5 Data and measurement performance

An example comparison between data and retrieved forward model from a single 24×24 m² ground cell is shown in Fig. 4e. The ground cell samples approximately the same radius – and thus the same wavelengths – twice as the ground location traverses the imaging field of view. However, it is evident that the signal at the same radius can have different values. This is primarily due to two effects: (1) the bidirectional scattering distribution function of the piece of terrain within a ground cell is sampled at two different observation angles during the measurement sequence, and (2) per-pixel signal changes due to satellite motion during the image integration time. The forward model accounts for this slow signal variation by replacing the constant albedo with a second-order polynomial that is a function of the image frame index, similar in concept to what is done in differential optical absorption spectroscopy (DOAS) retrievals (Plattand Stutz, 2008) where the polynomial is a function of wavelength. The forward model also accounts for the fact that the ground cell overlap with a given camera pixel changes between different frames during the observation sequence. Knowledge of the ground cell–camera pixel overlap is provided by the image co-registration algorithm and allows us to account for data effects that, if not treated, would result in erroneous high-frequency error. Residuals between data and the retrieved forward model with a standard deviation of 0.5 % are representative.

When the methane enhancement component of the state vector x is retrieved from the ∼200 000 ground cells within a measurement domain, a retrieved methane map can be plotted as in Fig. 5, which shows a selection of GHGSat-D measurements taken over various types of industrial sites around the world. In each retrieved methane enhancement field, we observed a localized methane plume whose point of origin coincides with a facility. For each observation, the plume is also displayed as an overlay on the retrieved SWIR surface reflectance image using thresholding and a spatial correlation criterion that counts a downwind enhancement as real if it is close (within a few pixels) to another enhancement that was previously determined to be real, with the requirement that the furthest upwind enhancement closely overlaps the location of a facility source. The retrieved SWIR reflectance can then be used to geolocate the retrieved methane enhancement field. The units of the excess methane vertical column density (VCD) are (mol m⁻²). For comparison, the nominal VCD background value is approximately 0.67 mol m⁻², corresponding to a column-averaged mixing ratio of XCH₄ ≈1.9 ppm. The peak plume enhancements in the examples presented are 20 %–120 % above background. In each observation, the origin of the methane plume can be clearly attributed to a ground location with an uncertainty of approximately 30 m. This location accuracy is sufficient to provide actionable information to facility operators.

Figure 5 Retrieved methane enhancement fields. In the left plot of each plot pair, the methane enhancement above the local background value is shown in a 7×7 km² region of interest centered on the plume. In the right plot of each plot pair, the extracted plume is overlaid on top of the retrieved surface reflectance. The methane enhancement colour scale is in units of mol m⁻². Description, location, and date of observations: (a) hydroelectric reservoir at the Lom Pangar dam, Cameroon, 20 April 2017; (b) suspected liquid unloading event in the Permian basin, Texas, USA, 17 October 2018; (c) underground coal mine vent near Camden, Australia, 18 October 2018; (d) underground coal mine vent near Farmington, New Mexico, USA, 18 September 2018; (e) natural gas compressor facility near Korpezhe, Turkmenistan, 24 February 2019; (f) same natural gas compressor facility near Korpezhe, Turkmenistan, 9 March 2019. In all retrieval fields, the plume enhancement competes against retrieval artefacts that are usually oriented along the direction of the satellite orbital direction.

In the example observations shown in Fig. 5, methane plumes were observed over a variety of industrial facilities: a hydro-electric reservoir, coal mine vents, and natural gas sites. On 20 April 2017 we observed a methane plume over the dam vanes of the Lom Pangar hydro-electric reservoir (Fig. 5a) in eastern Cameroon that was flooded the previous year. Hydro-electric reservoirs are a known source of methane and carbon dioxide emissions, especially those in tropical climates that have been recently flooded (Barros et al., 2011; Rosa et al., 2004). On 17 October 2018 we observed a methane plume over a natural gas facility in the Permian basin, TX, USA (Fig. 5b). This was a suspected liquid unloading event in which liquid in the well is removed to keep gas flowing to surface facilities, often resulting in a large, but temporary, pulse of methane emissions. On 18 October and 18 September 2018 we observed methane plumes over vents in coal mining operations near Camden, NSW, Australia (Fig. 5c), and Farmington, NM, USA (Fig. 5d), respectively. Methane can be released from coal and surrounding rock strata during mining operations, and large methane emissions have been observed for both of these sites in previous studies (Frankenberg et al., 2016; Ong et al., 2017). In the Farmington, NM, observation, we can see that the magnitude of the plume enhancement is commensurate with those of neighbouring enhancement artefacts in the methane retrieval field. This fact highlights the advantage of incorporating prior information about source locations when determining whether a measured enhancement at or near the detection limit is real. On 24 February and 9 March 2019 we observed large methane plumes over the Korpezhe oil–gas field in western Turkmenistan (Fig. 5e, f). Large methane emissions from this area have been previously reported in Varon et al. (2019).

The error in methane retrieved from a single-pass GHGSat-D observation is typically between 8 %–25 % of the background value, depending on factors such as the time of year and complexity of the albedo field. This error is due to various GHGSat-D instrument imperfections, including out-of-field straylight, in-field optical ghosting, and memory lag effects in the camera response. Minor imperfections and uncertainties in our instrument model can lead to significant systematic errors in the methane retrievals. One of the most obvious errors in the methane retrieval is “streaking” in the direction of the along-track satellite motion. The high-frequency streaking is likely due to camera pixels which are dead, “hot” (always on), blinking, or have mischaracterized offset or gain that have not been flagged by the bad-pixel detector. The low-frequency character of the streaking can be explained by unwanted optical signal (from straylight or ghosting) coupling to the spectral dips, especially those that occur near tangents of the spectral rings. It should also be noted that the measurement precision varies along the cross-track dimension of the retrieval domain due to different sampling of the spectral absorption rings.

It is useful to compare the observed measurement error levels to the limit set by random noise on the camera. The GHGSat-D per-pixel signal-to-noise ratio is $\text{SNR}=I_{\mathrm{0}}/\mathit{\sigma }=\mathrm{200}$, with I₀ the mean per-pixel signal level and σ² the sum of the mean photon shot noise, camera dark noise, and read noise variance. To compare to the methane measurement error, we define the methane SNR, ${\text{SNR}}_{{\mathrm{CH}}_{\mathrm{4}}}$, using the posterior error $\widehat{\mathbf{S}}={\left({\mathbf{K}}^{\mathbf{T}}{\mathbf{S}}_{\mathrm{o}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}$ derived from optimal estimation theory (Rodgers, 2000), which has a diagonal error covariance matrix S₀ populated with σ². Our definition of ${\text{SNR}}_{{\mathrm{CH}}_{\mathrm{4}}}$ is the inverse square root of the molecule state parameter element in ${\widehat{S}}^{-\mathrm{1}}$ with all elements of S_a set to zero except for methane, which is set to infinity. This describes the situation where we have perfect knowledge of all non-methane state vector elements, and our methane knowledge is determined solely from the data:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\text{SNR}}_{{\mathrm{CH}}_{\mathrm{4}}}& ={\mathrm{\Omega }}_{{\mathrm{CH}}_{\mathrm{4}}}\cdot {\left({\mathit{\sigma }}^{\mathrm{2}}{\left({\mathbf{K}}_{{\mathrm{CH}}_{\mathrm{4}}}^{\mathbf{T}}\cdot {\mathbf{K}}_{{\mathrm{CH}}_{\mathrm{4}}}\right)}^{-\mathrm{1}}\right)}^{-\frac{\mathrm{1}}{\mathrm{2}}}\\ & =\text{SNR}\cdot {\mathrm{\Omega }}_{{\mathrm{CH}}_{\mathrm{4}}}\cdot {\left({\mathbf{k}}_{{\mathrm{CH}}_{\mathrm{4}}}^{\mathbf{T}}\cdot {\mathbf{k}}_{{\mathrm{CH}}_{\mathrm{4}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where ${\mathrm{\Omega }}_{{\mathrm{CH}}_{\mathrm{4}}}$ is the total methane VCD, ${\mathbf{K}}_{{\mathrm{CH}}_{\mathrm{4}}}$ is the methane Jacobian, and ${\mathbf{k}}_{{\mathrm{CH}}_{\mathrm{4}}}={\mathbf{K}}_{{\mathrm{CH}}_{\mathrm{4}}}/I_{\mathrm{0}}$ is the Jacobian normalized by signal level. This expression can also be derived from a linear least-squares method (Adler et al., 2010). For a nominal albedo of 0.2, solar zenith angle of 40^∘, and US standard temperature, pressure, and mixing ratio profiles (U.S. COESA, 1976) scaled to present values of the total column densities, ${\text{SNR}}_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{15.4}$. The 8 %–25 % methane measurement error can be compared to the ideal limit ${\text{SNR}}_{{\mathrm{CH}}_{\mathrm{4}}}^{-\mathrm{1}}=\mathrm{6.5}$ %. We see that the lowest-error GHGSat-D observations approach the ideal limit. These tend to be observations over bright scenes with quasi-uniform surface reflectance, such as deserts or plains. Conversely, observations with larger error in the methane retrieval tend to occur over areas with highly non-uniform and/or low surface reflectance, such as urban scenes.

Given the range of column error levels of 8 %–25 %, our experience with source rate retrievals (Varon et al., 2018, 2019) suggests that GHGSat-D is sensitive to point emitters with Q>1000–3000 kg h⁻¹. Using a typical value of wind speed (3 m s⁻¹), this is consistent with a simple model for point source detection threshold based on excess methane in the source pixel column (Jacob et al., 2016).

6 Summary and future plans

GHGSat has developed, built, and successfully launched a demonstration satellite, GHGSat-D, continuously operational since 2016. GHGSat-D uses a wide-angle Fabry–Pérot imaging spectrometer to make quantitative measurements of the methane column density, with a focus on resolving enhancements above background values within the 12×12 km² measurement field of view. Since beginning on-orbit measurements in 2016, GHGSat-D has made over 5000 observations including several demonstrated discoveries of industrial methane emissions from space (Varon et al., 2019). This mission has proven that a compact spectrometer on a small satellite can be used to detect and quantify methane plumes from individual facilities with unambiguous attribution. While its detection threshold is estimated to be relatively high at 1000–3000 kg h⁻¹, the experience of designing, manufacturing, and operating GHGSat-D has been a highly fruitful process. Detailed investigations of the retrieval outputs and comparison with simulations under various conditions have helped us understand the limiting sources of error and informed the design of our next satellites, which have much better projected performance.

The first satellite in GHGSat's commercial constellation, GHGSat-C1 or “Iris”, was successfully launched on 2 September 2020. The second commercial satellite, GHGSat-C2 or “Hugo”, is scheduled to launch in December 2020. At the time of this paper's submission, GHGSat-C1 was in its commissioning phase and had detected its first emission plumes from industrial facilities. GHGSat-C1 has an improved design informed by lessons learned from GHGSat-D. Most importantly, GHGSat-C1 has a 100× reduction (approx.) in straylight magnitude, a 5× reduction in ghosting magnitude, increased per-pixel signal levels, an on-board dark and flat-field calibration system, and a re-optimized spectroscopic configuration. GHGSat-C1 has also undergone an intensive test and characterization campaign in which camera and instrument behaviour have been more extensively explored, calibrated, and parameterized than was done for GHGSat-D. The retrieval method has also been advanced, including significant improvements in alignment and spatially resolved column retrievals, tested using aircraft trials ahead of the launch of GHGSat-C1. GHGSat will perform a calibration and validation campaign for GHGSat-C1 that includes, among other activities, several controlled methane release campaigns. We estimate that we will achieve column errors of ∼2 % with GHGSat-C1 – including systematic errors – for a nominal observation (subsequent satellites will have similar or better performance). Given these improvements in the column precision and finer spatial resolution (∼25 m expected), we anticipate satellites in our constellation to achieve detection thresholds at or below 100 kg h⁻¹ for nominal conditions. As we scale up the constellation, this will allow us to provide increasing amounts of high-fidelity, actionable data to industrial operators worldwide, ultimately leading to significant emissions reductions.