Quantitative imaging of carbon dioxide plumes using a ground-based shortwave infrared spectral camera

We present the first results of a ground-based imaging experiment using a shortwave infrared spectral camera to quantify carbon dioxide (CO₂) emissions from a coal-fired power plant in Mannheim, Germany. The power plant emits more than 4.9 Mt CO₂ yr⁻¹ and is a validation opportunity for the emission estimation technique. The camera is a hyperspectral imaging spectrometer that covers the spectral range from 900 to 2500 nm with a spectral resolution of 7 nm. We identify CO₂ enhancements from hourly averaged images using an iterative matched filter retrieval using the 2000 nm absorption band of CO₂. We present 11 plume images from 5 d in 2021 and 2022 covering a variety of ambient conditions. We design a forward model based on a three-dimensional, bent-over Gaussian plume rise simulation and compare our observed emission plumes with the forward model. The model depends on the parameters ambient wind velocity, wind direction, plume dispersion, and emission rate. We retrieve the emission rate by minimizing the least-squares difference between the measured and the simulated images. We find an overall reasonable agreement between the retrieved and expected emissions for power plant emission rates between 223 and 587 t CO₂ h⁻¹. The retrieved emissions average 84 % of the expected emissions and have a mean relative uncertainty of 24 %. The technique works at wind speeds down to 1.4 m s⁻¹ and can follow diurnal emission dynamics. We also include observations with unfavorable ambient conditions, such as background heterogeneity and acute observation angles. These conditions are shown to produce considerable biases in the retrieved emission rates, yet they can be filtered out reliably in most cases. Thus, this emission estimation technique is a promising tool for independently verifying reported emissions from large point sources and provides complementary information to existing monitoring techniques.

1 Introduction

Carbon dioxide (CO₂) is the most important anthropogenic greenhouse gas driving climate change (Masson-Delmotte et al., 2021). The signatory countries of the Paris Agreement set ambitious goals for CO₂ emissions reductions, and they agreed on implementing a stock taking mechanism to supervise mitigation progress. In consequence, the international science community has elaborated a plan to build a monitoring and verification support (MVS) capacity for anthropogenic CO₂ emissions (Janssens-Maenhout et al., 2020). Besides refining inventory-making and atmospheric modeling, the plan includes further developing and implementing atmospheric measurement techniques that can help quantify anthropogenic emissions from global to local scales using in situ and remote sensing observations with the goal to independently verify reported emissions and to monitor the effectiveness of reduction measures.

The local scales have received particular attention recently since spectroscopic imaging techniques are emerging that enable quantification of localized anthropogenic emission sources. Pioneering space missions such as the Greenhouse Gases Observing Satellite (GOSAT) and the Orbiting Carbon Observatory (OCO-2/-3) were able to quantify the urban CO₂ domes of megacities such as Los Angeles (Kort et al., 2012; Schwandner et al., 2017; Kiel et al., 2021). OCO-2 was the first satellite to deliver images of CO₂ plumes from individual coal-fired power plants (Nassar et al., 2017), and since then, manifold activities have demonstrated that spectroscopic imaging with fine spatial resolution from satellites and aircraft can deliver facility-scale emission estimates for CO₂ (Thorpe et al., 2017; Cusworth et al., 2021a; Fujinawa et al., 2021) and methane (CH₄) (Frankenberg et al., 2016; Duren et al., 2019; Guanter et al., 2021). Emission estimation methods either use a mass-balance approach (Liu et al., 2021; Varon et al., 2018, 2020), Gaussian plume models (Nassar et al., 2017; Schwandner et al., 2017; Varon et al., 2018), or machine learning (Jongaramrungruang et al., 2022). Many activities aim at designing next-generation satellite missions with ground resolution on the order of a few tens of meters and enhanced quantification capabilities for CO₂ and CH₄ (Strandgren et al., 2020; Jacob et al., 2022). On the ground, pilot studies and field campaigns were able to constrain city-scale and localized emissions sources using various techniques (Christen, 2014) ranging from eddy-covariance methods (Crawford and Christen, 2015; Christen et al., 2011) to local in situ and remote sensing concentration measurements deployed on mobile platforms and in small ad hoc networks (Hase et al., 2015; Luther et al., 2022). Gålfalk et al. (2016) report on the development of a CH₄ camera that operates in the thermal infrared, enabling imaging of CH₄ released from localized sources. Our precursor study (Knapp et al., 2023) demonstrates the methane imaging capabilities of a stationary spectral camera in the shortwave infrared (SWIR) and derives sub-hourly emission estimates for coal mine ventilation shafts.

Here, we present the first results from our ground-based SWIR camera (HySpex SWIR-384) on imaging of CO₂ plumes from strong point sources. The camera collects sunlight scattered in the sky, and the retrieval exploits the CO₂ absorption structures around 2000 nm wavelength, which is analog to previously reported aircraft and satellite techniques that work on reflected sunlight (Thorpe et al., 2017). Thus, the measurements only work in daytime and require fair weather conditions. Observing skylight in the SWIR further comes with the disadvantage that the sky, while being spectrally smooth, is dark, and thus, our technique currently requires co-adding on hourly timescales. This, however, might be improved by future developments in sensor technology. We develop an emission estimation technique based on bent-over Gaussian plume modeling and use the well-known emissions of a hard-coal power plant to validate our results. Such ground-based CO₂ cameras could be particularly useful to monitor emissions of point-like sources such as coal-fired power plants, cement factories, industrial facilities, or natural sources such as volcanoes. These observations would be a valuable addition to our MVS capacities per se and could be used together with network techniques to disentangle contributions from various sources in complex emission landscapes. We also suggest that imaging CO₂ plumes might be a tool to raise public awareness on the scale of local greenhouse gas pollution and the urgency to implement emission reductions (Jungmann et al., 2022).

For demonstration purposes, we deploy the HySpex hyperspectral camera together with a portable wind lidar at a few kilometers' distance from a medium-sized hard-coal power plant (> 4.9 Mt CO₂ yr⁻¹) located in Mannheim, Germany (Sect. 2 for measurement setup). Repeatedly scanning the sky above the power plant, we retrieve CO₂ column enhancements in the power plant plume using a matched filter algorithm similar to previous satellite (Guanter et al., 2021) and airborne studies (Foote et al., 2020, 2021) (Sect. 3.2 for data analysis). Based on the plume enhancements, we estimate the emissions of the power plant from the observed plumes. To this end, we design a forward model to match our observations based on the plume rise model IBJpluris (Janicke and Janicke, 2001). The best match between the measured and simulated observations feeds our emission estimate (Sect. 4 for the emission estimation method). In total, we collect 11 plumes over 5 d, for which we compare our emission estimates to the emissions of the power plant calculated from the instantaneous power production. Finally, we discuss the capabilities and limitations of the technique (Sect. 5 for the results).

2 Instrumentation and field deployment

2.1 HySpex SWIR-384 camera

The HySpex SWIR-384 camera is a commercially available hyperspectral camera by Norsk Elektro Optikk^® (NEO). The camera optics focuses the incoming light onto a slit, which passes light from a horizontal opening angle of 7.3 mrad width. After collimation, a grating disperses the incoming light, and a mirror focuses the spectrum on a 2D detector array. The detector array samples the vertical dimension of the scene with 384 pixels (384 “lines”) and the spectral dimension with 288 pixels (“288 channels”). The lines of the detector cover a vertical field of view of 16°, and each pixel covers a solid angle of approximately 7.3 mrad × 7.3 mrad. The channels record the spectrum between 950 and 2500 nm with a sampling of 5.45 nm. Each read-out of the 2D detector is called a frame. The camera is mounted on a rotation stage, which turns it in a horizontal (azimuth) direction, thus scanning over the target scene by collecting a sequence of frames. The rotation pattern is clockwise in steps that correspond to the acceptance angle of a single pixel. A shutter closes the camera aperture prior to and after each scene scan in order to take 200 dark spectra. The detector is cooled to 147 K during operation to reduce dark current. We refer to taking a hyperspectral image as a scan. The raw output is in analog–digital converter units with a 16-bit resolution, representing radiances between 1 and 2¹⁶−1, and the datacube has dimensions of number of frames × number of lines × number of channels. A tripod carries the rotation stage with the camera during field deployment. A rugged, field-deployable GETAC^® laptop controls the camera and the rotation stage. The camera and the laptop each run on their own battery packs, which provide enough power for more than 6 h of consecutive measurements. Figure 1 shows the camera deployed in the field on 6 September 2021 ​​​​​​​ observing the sky above the power plant in Mannheim.

Figure 1 Photograph of the HySpex SWIR-384 camera deployed in the vicinity of the hard-coal power plant Grosskraftwerk Mannheim on 6 September 2021. In this particular case, the field of view (red frame) covers the sky above two stacks of the power plant. The instrument scans the scene from left to right. The observed chimneys are 200 m (left) and 180 m (right) high and are approximately 3.2 km away from the camera. The inset panel shows a typical sky spectrum in the entire accepted spectral range from roughly 900 to 2500 nm. Notice the absorption features of carbon dioxide at around 2000 nm.

The HySpex camera comes with a radiometric and spectral calibration by the manufacturer. For calibrating the instrument spectral response function (ISRF), we used a tunable diode laser operating around 1600 nm in our lab, similar to the setup of Lenhard et al. (2015). The ISRF has a Gaussian-like profile with a full width at half maximum (FWHM) of 6.7–7 nm; see Fig. 2. This is in agreement with the manufacturer calibration results of approximately 7 nm for the 1694 nm Argon emission line. We use a radiative transfer model similar to Guanter et al. (2021) for the spectral calibration of our fit interval and correct an offset of +2.25 nm to the manufacturer's calibration. Figure 1 shows a typical clear-sky spectrum with the entire recorded spectral range.

Figure 2 The instrument spectral response function of channel 119 at selected lines across the detector. Panel (a) shows the measured radiance and a fitted Gaussian. Each data point represents a mean and standard deviation of 286 frames. Panel (b) shows the residuum of the data points to the Gaussian fit.

2.2 Attitude and heading reference system

We constrain the observation geometry during measurements with an attitude and heading reference system (AHRS). An MTi-7 miniature GNSS/INS module from XSENS^® was mounted on the HySpex camera such that it rotates with the camera during observation. The device performed well in previous campaigns in shipborne applications (Dörner et al., 2018). We conducted performance tests on its inertial navigation system (INS) and found its 10 min precision well below 0.05° for the instrument pitch and roll angle. All data transfer in real time to the GETAC^® laptop controlling the HySpex camera. The sensor does not provide reliable data on the instrument viewing azimuth angle (VAA) in stationary operation. We find the VAA of each image frame by identification of a distinct landmark within the image, e.g., a chimney, and assigning the forward azimuth angle from the camera location to the landmark. The VAA of all other frames in the image follows from the horizontal opening angle of the camera. For typical distances between the camera and the landmark, the uncertainty of the VAA is below 1°. Additionally, the sensor provides us with GPS/GNSS-based geolocation data. We use the software of Holmgren et al. (2018) to calculate the solar zenith angle (SZA) and solar azimuth angle (SAA) from the geolocation and time of the observation.

2.3 Windranger 200 lidar

The Windranger 200 is a compact and lightweight wind lidar produced by METEK^®. Adler et al. (2021) used it to successfully study boundary layer turbulence. The instrument weighs approximately 50 kg and has dimensions of 840 mm × 540 mm × 580 mm. For absolute reference of wind direction, the lidar needs to point to the north or another known azimuth reference point. The lidar works with frequency-modulated continuous-wave (FMCW) technology (Peters, 2018). It uses the Doppler-shifted back-reflection of a 1545 nm laser for the measurements. The laser beam rotates by 360° once per second, producing a data point of wind velocity and direction at one height layer. It measures at six height levels up to 200 m. In our setting, the layers are at 10, 20, 50, 100, 150, and 200 m. The method is sensitive to turbulent wind fluctuations within one rotation, but it provides a complete profile approximately every 10 s, such that statistical fluctuations can be reduced by averaging over several rotations.

2.4 Field measurements

The HySpex operates in a 3 to 4 km distance from the target; thus an image has a spatial resolution of approximately 2 to 3 m. This provides us with a plume image that is spatially well resolved, while we measure enough background-sky pixels for the retrieval of enhancements (Sect. 3.2). We level the camera based on the INS real-time data prior to each observation by collecting 10 min means of the INS pitch and roll and adjusting the camera's alignment for potential offsets (Sect. 2.2). This ensures a roll angle of 0.0 ± 0.1°. Subsequently, the Viewing Elevation Angle (VEA) is adjusted by tilting the camera such that the lower edge of the image contains the upper part of the CO₂ source, e.g., the chimney tip, and the upper part of the image shows the sky (Fig. 1). Typical single exposures of the detector range between 8 and 20 ms. Each frame adds 5–10 of these exposures before the camera rotates a step in the azimuth direction. Scans that cover azimuth angles of 10 to 15° have typical scan times of 90 s. We record consecutive scans as long as atmospheric conditions like cloud cover and solar illumination are favorable, i.e., a bright and clear sky. To achieve a sufficient signal-to-noise ratio, the hyperspectral datacubes submitted to the retrieval are co-additions of scans of at least 55 min.

The lidar is positioned next to the camera. Its power supply is a Jackery^® Explorer 500 battery with a Jackery^® SolarSaga 100W solar panel. During field observations, we monitor the wind data quality and re-calibrate the lidar height levels if data quality decreases. We measure wind velocities and directions in six levels from 10 to 200 m height.

Table 1 The first four columns list the date, the time interval of the observations, the unit number of the power plant that was observed, and if the Windranger 200 lidar was available. Columns 5 and 6 list the mean viewing azimuth angle (VAA) of the camera and the distance d to the observed unit. Columns 7 and 8 list the AERONET aerosol optical thickness (AOT) at 2000 nm and the asymmetry parameter (g) of the scattering phase function; see Sect. 3.3. Column 9 lists if condensation was observed in the plume. For more detailed atmospheric conditions, see Table S2.

Here, we present observations taken on 5 d in 2021 and 2022. We headed to the field when weather forecasts predicted clear skies and targeted the largest section of the power plant, unit 9, preferentially. Table 1 lists the observation periods alongside a short description of atmospheric conditions and the instrumental setup. Atmospheric conditions varied over the day and between days concerning cloud cover, wind conditions, and aerosol load. For more detailed information on the atmospheric conditions, see Sect. 4.3 and Table S2 in the Supplement. We process scans that pass a visual quality filter for heavy clouds or sporadic obstruction by pedestrians or cars. Time intervals with favorable conditions vary between 100 and 258 min d⁻¹, and we retrieve plume images for averaged images of approximately 1 h, as described in Sect. 3.1.

We observed unit 9 on all days except 13 May 2022, on which unit 9 malfunctioned, and we observed unit 6 instead.

3 Carbon dioxide retrieval

3.1 Preprocessing the images

For each scan, we term the raw detector output of the camera DN_ijk, where i labels the channel, j the line, and k the frame. We take 200 background frames prior to and after each scan. Background frames are taken with a closed shutter in front of the detector. We calculate a mean detector background BG_ij before and after the scan as the average of the background frames. The backgrounds are linearly interpolated to the frames k during observation for each detector element and subtracted from each frame. The background-corrected image L_ijk is therefore given by

$$\begin{array}{}\text{(1)}& L_{ijk}={\mathrm{DN}}_{ijk}-{\mathrm{BG}}_{ij}^{\mathrm{before}}-{\displaystyle \frac{{\mathrm{BG}}_{ij}^{\mathrm{after}}-{\mathrm{BG}}_{ij}^{\mathrm{before}}}{\mathrm{number}\mathrm{of}\mathrm{frames}}}\cdot k.\end{array}$$

For later convenience, we chain the spectral dimension i into a vector ${\mathit{L}}_{jk}=[L_{\mathrm{1}jk},L_{\mathrm{2}jk}\mathrm{\cdots }{L}_{\mathrm{288}jk}]$.

Malfunctioning detector elements (pixels on the detector) are identified in laboratory measurements using halogen lamps fed into an integrating sphere (40 cm diameter). Elements are removed from further processing if they exhibit an exceptionally high, low, or variable response to the broadband illumination. The missing values are interpolated in the spatial (line) direction. Field deployment of the HySpex camera typically produces several hundred single scans. Each scan is corrected according to Eq. (1) and visually inspected for quality control. We remove images that are corrupted by obstacles obscuring the field of view in single frames or by significant changes in the overall atmospheric conditions, e.g., cloud formation. The scans remaining after quality control are averaged. Some observed emission plumes show condensation from the water vapor co-emitted with the CO₂. The light paths of photons under clear sky are substantially different from those under condensate conditions. To avoid intermingling light path differences with concentration differences, we identify and exclude pixels with condensate from the averaged image. The saturation of each pixel in a scan is defined as the maximum of the spectrum in the pixel $S_{jk}=\underset{i}{max}\left(L_{ijk}\right)$. We calculate the background saturation BS_jk of the sky for each scan by

$$\begin{array}{}\text{(2)}& {\mathrm{BS}}_{jk}=\underset{j}{\mathrm{med}}\left(S_{jk}\right)\cdot {\left({\displaystyle \frac{\underset{k}{\mathrm{med}}\left(S_{jk}\right)}{{〈\underset{k}{\mathrm{med}}\left(S_{jk}\right)〉}_j}}\right)}^T,\end{array}$$

where med is the median operator. Equation (2) calculates the clear-sky saturation for each pixel in a scan from the scan itself, assuming that the plume condensate covers a small part of the sky. We subtract BS_jk from the scan saturation, such that the residual image scatters around zero where a clear sky was observed. The use of the median assures that the high reflectivity of the condensed water does not affect BS_jk. Pixels are masked in scans if their saturation deviates more than +3σ from the residual distribution since this shows an exceptionally high reflectivity caused by plume condensation or clouds. We build the averaged image from the spectra of the single scans, which were not removed by the saturation mask. Therefore, the averaged image contains pixels to which not all scans have contributed. We exclude pixels from the retrieval of the averaged image if less than 90 % of the scans contribute to the average image. Figure 3 shows how the mask based on Eq. (2) identifies over-saturated pixels in single scans and how the averaged image is constructed from there.

Figure 3 (a) Saturation of a single scan focused on the part of the image with a condensate plume. The image was taken on 23 March 2022 in Mannheim. The red border surrounds all pixels that are excluded from averaging. (b) Average image with colored lines marking the areas in which a certain fraction of scans contributed to the spectra in the pixels. The black line marks a mask (generated by eye), which excludes the chimney from further processing. Note that the image dimensions are VAA and VEA.

We use the geolocation data of the camera and the target to calculate distance d and the viewing azimuth angle of the instrument based on the WGS84 reference system (Slater and Malys, 1998). The observed chimney also provides us with a line of known height h in the image; thus the viewing elevation angle is found by

$$\begin{array}{}\text{(3)}& \mathrm{VEA}=\arctan \left({\displaystyle \frac{h}{d}}\right).\end{array}$$

The geometric area of each image pixel A_j at the distance d to the target source is

$$\begin{array}{}\text{(4)}& A_j=\left(d\cdot \tan \left(\mathrm{\Delta }{\mathrm{VEA}}_j\right)\right)\cdot \left(d\cdot \tan \left(\mathrm{\Delta }\mathrm{VAA}\right)\right),\end{array}$$

where ΔVAA=0.73 mrad is the horizontal opening angle of the camera, and ΔVEA_j≈0.73 mrad is the vertical opening angle of each line.

Inhomogeneities in detector response or the optical setup cause striping patterns in imaging spectrometer data (Borsdorff et al., 2019). We compensate for this by dividing the spectral vectors L_jk element-wise by a reference spectrum ${\widehat{\mathit{L}}}_j$ from the same line. Thus, our normalized spectral vectors l_jk are given by

$$\begin{array}{}\text{(5)}& {\mathit{l}}_{\mathit{j}\mathit{k}}={\displaystyle \frac{{\mathit{L}}_{jk}}{{\widehat{\mathit{L}}}_j}},\end{array}$$

and scatter around unity if the scene was homogeneous. The reference spectrum is taken from each scan itself as the mean of 20 frames upwind of the source.

3.2 Matched filter retrieval

The matched filter is a data-driven statistical approach for signal identification and quantification in noisy data, e.g., hyperspectral images (Dennison et al., 2013; Manolakis et al., 2014; Zhang et al., 2020; Cusworth et al., 2021b; Guanter et al., 2021). It estimates a spectral background variability from two-dimensional spatial tiles of absorption spectra and identifies a pre-defined spectral signature exceeding this variability, i.e., in our case, the spectral absorption signature of CO₂ around 2000 nm wavelength. Foote et al. (2020) improve on the classical matched filter by introducing an albedo correction and a sparsity constraint on the enhancements, which we also used in our precursor study for CH₄ (Knapp et al., 2023). The sparsity constraint sets enhancements below the detection limit to zero, which enables the matched filter to iteratively remove the target signal from the background clutter estimation. This is particularly important for the detection of weak signals in the presence of strong background signals (Foote et al., 2020; Schaum, 2021; Pei et al., 2023). Heterogeneities in the observed scene are particularly challenging for statistical retrievals (Ayasse et al., 2018). Thus, the homogeneous reflectivity of the clear sky is an advantage for the matched filter, although the sky is rather dark in the shortwave infrared spectral range.

The classic matched filter (CMF) calculates the enhancement α in each pixel as

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{jk}={\displaystyle \frac{{\left({\mathit{l}}_{jk}-\mathit{\mu }\right)}^T{\mathbf{C}}^{-\mathrm{1}}\mathit{t}}{{\mathit{t}}^T{\mathbf{C}}^{-\mathrm{1}}\mathit{t}}},\end{array}$$

where l_jk is the normalized spectral vectors defined in Eq. (5), μ=〈l〉_jk is the mean normalized spectral vector across the image, and C=cov(l_jk) is the covariance matrix of all the normalized spectral vectors across the image. The target signature t contains the spectral feature the matched filter retrieval is sensitive to. In our case, it follows from the absorption cross-section of CO₂. The target signature is commonly derived from the assumption that an enhancement α in a pixel acts on a spectrum L according to Beer–Lambert's law:

$$\begin{array}{}\text{(7)}& \mathit{L}={\mathit{L}}_{\mathrm{0}}\cdot \exp \left(-\mathit{\alpha }\cdot \mathit{s}\right),\end{array}$$

where s is the unit absorption spectrum of the gas (Sect. 3.3). We rearrange Eq. (7) for referenced spectra and linearize the exponential function via Taylor expansion to first order, which yields

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{0}}}}\approx \mathrm{1}-\mathit{\alpha }\cdot \mathit{s}.\end{array}$$

The target signature t of the matched filter is given by the derivative of Eq. (8) with respect to α, which is the unit absorption spectrum s.

We adapt our retrieval from the MAG1C (Matched filter with Albedo correction and reweiGhted L1 sparsity Code) algorithm of Foote et al. (2020), which decreases enhancement uncertainty and background noise. The albedo correction reduces the effect of background heterogeneity, which accounts for the gradient in the brightness of the clear sky in our case. We calculate a scalar factor

$$\begin{array}{}\text{(9)}& r_{jk}={\displaystyle \frac{{\mathit{L}}_{jk}^T\widehat{{\mathit{L}}_j}}{{\widehat{{\mathit{L}}_j}}^T\widehat{{\mathit{L}}_j}}}\end{array}$$

for each pixel, which gives its brightness relative to the line reference spectrum. The sparsity prior adds an L0 regularization to the cost function of the retrieval to minimize the number of detected enhancements since we expect enhancements in less than 1 ‰ of the pixels. Foote et al. (2020) further introduce a positivity constraint on α since gas enhancements in the atmosphere are non-negative. Thus, Eq. (6) transforms to an iterative retrieval given by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\alpha }}_{jkn}=max\left({\displaystyle \frac{{\left({\mathit{l}}_{jk}-{\mathit{\mu }}_n\right)}^T{\mathbf{C}}_n^{-\mathrm{1}}{\mathit{s}}_j-w_n}{r_{jk}{\mathit{s}}_j^T{\mathbf{C}}_n^{-\mathrm{1}}{\mathit{s}}_j}},\mathrm{0}\right),\\ & w_n={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{n-\mathrm{1}}+\mathit{\kappa }}},\end{array}\end{array}$$

where w_n is the pixel's regularization weight of the sparsity prior, n the iteration step, and κ a small number for numerical robustness. The target signature is given by the unit absorption spectrum s_j of the gas, which depends on the image line j (Sect. 3.3). The iteration is necessary due to the sparsity constraint (Candès et al., 2008). It allows for iteratively improving on the estimate of the background mean and covariance, μ and C. We perform 30 iteration steps for each scene to get an accurate estimate of the background distribution. Pixels with low abundances are artificially forced to zero in the MAG1C algorithm by Foote et al. (2020). To avoid a systematic bias from these pixels in the emission inversion (Sect. 4.4), we perform the final iterative retrieval step without the positivity and sparsity constraint from Eq. (10). This ensures that the background estimation (μ,C) benefits from the iterative retrieval, while the retrieved enhancements are non-zero in the background pixels of the image. Thus, we obtain a final estimate of the enhancement α_jk in every pixel. The uncertainty of the matched filter retrieval is given by the retrieval error covariance matrix (Köhler et al., 2015) and depends on s and C. It can be calculated by

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_j^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{{\mathit{s}}_j^T{\mathbf{C}}^{-\mathrm{1}}{\mathit{s}}_j}},\end{array}$$

which gives an enhancement variance ${\mathit{\sigma }}_j^{\mathrm{2}}$ for each line of the image.

3.3 Unit absorption spectrum

In our units, the unit absorption spectrum s defined in Eq. (8) is the relative change of observed radiance due to a 1 ppm increase in the atmospheric mixing ratio of CO₂ over a path length of 1 m. The target spectrum can be calculated from the absorbing molecules' absorption cross-section by radiative transfer calculations and convolution with the instrument function (Thorpe et al., 2013; Thompson et al., 2015). Foote et al. (2021) show that a scene-specific spectral signature greatly improves the quality of the retrieved data. Thus, we choose to simulate the change of radiance with the atmospheric enhancement of CO₂ for a ground-based observer. An enhancement α of the atmospheric column of an absorber increases the optical thickness such that the at sensor radiance (ASR), I changes according to Beer's law:

$$\begin{array}{}\text{(12)}& I(\mathit{\lambda };\mathit{\alpha })=I_{\mathrm{0}}\left(\mathit{\lambda }\right)e^{-\mathit{\tau }\left(\mathit{\alpha }\right)},\end{array}$$

where τ is the optical thickness due to α and I₀ the observed radiance of an atmosphere without an enhancement. This assumes that there is no change in the photon light path induced by the absorber or a co-emitted substance like water vapor. To simulate the effect of an additional column of carbon dioxide on the spectrum, we calculate the optical thickness τ as

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\tau }& {\displaystyle }=n_{{\mathrm{CO}}_{\mathrm{2}}}\cdot {\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}}\cdot \mathrm{\Delta }z\text{(14)}& {\displaystyle }& {\displaystyle }=n_{\mathrm{air}}\cdot {\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}}\cdot \mathit{\alpha }\cdot {\mathrm{10}}^{-\mathrm{4}},\end{array}$$

where $n_{{\mathrm{CO}}_{\mathrm{2}}}$ and n_air are the number densities of carbon dioxide and ambient air (in molec. cm⁻³), respectively; ${\mathit{\sigma }}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the absorption cross-section of carbon dioxide (in cm² molec.⁻¹); Δz is the path length (in cm); α is the enhancement (in ppm m); and 10⁻⁴ is a unit conversion factor (from cm² to m²). Both I₀(λ) and I(λ;α), assume infinite spectral resolution of the instrument i.e., they are line-by-line quantities. Thus, we need to convolve the ASR with the ISRF of the instrument to get simulated measurements F(λ;α),

$$\begin{array}{}\text{(15)}& F(\mathit{\lambda };\mathit{\alpha })=\int I({\mathit{\lambda }}^{\prime };\mathit{\alpha })\cdot \mathcal{N}(\mathit{\lambda }-{\mathit{\lambda }}^{\prime };\mathrm{FWHM})\mathrm{d}{\mathit{\lambda }}^{\prime },\end{array}$$

where 𝒩 is the Gaussian kernel with an FWHM of 7 nm. From Eq. (7), it follows that

$$\begin{array}{}\text{(16)}& s\left(\mathit{\lambda }\right)=-{\displaystyle \frac{\partial }{\partial \mathit{\alpha }}}\ln \left({\displaystyle \frac{F(\mathit{\lambda };\mathit{\alpha })}{F_{\mathrm{0}}\left(\mathit{\lambda }\right)}}\right),\end{array}$$

where s is the unit absorption spectrum in units of (ppm m)⁻¹. It is binned to each channel i according to the spectral calibration of the instrument.

We calculate the ASR using the single scattering approximation of the radiative transfer equation (RTE) for an upward looking observer. Since we are only interested in a relative change of radiance, we simplify the RTE by neglecting thermal emission and multiple scattering. We consider scattering on molecules using Rayleigh theory and scattering on aerosols via a Henyey–Greenstein phase function (Henyey and Greenstein, 1941). These assumptions lead to an analytically solvable RTE of the form

$$\begin{array}{}\text{(17)}& \begin{array}{rl}I\left(\mathit{\lambda }\right)=& S_{\mathrm{0}}\left(\mathit{\lambda }\right)\cdot {\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }\left(\right|\mathit{\gamma }|-|{\mathit{\gamma }}_{\mathrm{0}}\left|\right)}}\\ \cdot & \sum _k^{\mathrm{layer}}{\stackrel{\mathrm{̃}}{\mathit{\omega }}}_k\left(\mathit{\lambda }\right)p_k(\mathit{\lambda };\mathrm{\Omega },{\mathrm{\Omega }}_{\mathrm{0}})e^{-\frac{{\mathit{\tau }}_k\left(\mathit{\lambda }\right)}{\left|{\mathit{\gamma }}_{\mathrm{0}}\right|-\left|\mathit{\gamma }\right|}}\\ & \left(\mathrm{1}-e^{-\mathrm{\Delta }{\mathit{\tau }}_k\left(\frac{\mathrm{1}}{\left|{\mathit{\gamma }}_{\mathrm{0}}\right|}-\frac{\mathrm{1}}{\left|\mathit{\gamma }\right|}\right)}\right),\end{array}\end{array}$$

where S₀ is the top of the atmosphere (TOA) radiance, ${\mathit{\gamma }}_{\mathrm{0}}=\frac{\mathrm{1}}{\cos \left(\mathrm{SZA}\right)}$, $\mathit{\gamma }=\frac{\mathrm{1}}{\cos \left(\mathrm{VZA}\right)}$, $\stackrel{\mathrm{̃}}{\mathit{\omega }}$ is the single scattering albedo, and p is the scattering phase function for the beam direction before (Ω₀) and after (Ω) the single scattering. Equation (17) holds for an atmosphere consisting of plane-parallel, horizontally homogeneous layers with an optical thickness Δτ_k per layer k and a total atmospheric optical thickness down to a layer k of τ_k.

We solve Eq. (17) between 1900 and 2300 nm at a resolution of 0.001 nm to calculate I₀(λ). The model uses absorption cross-sections of CO₂, CH₄, and water vapor from the HITRAN2016 database (Gordon et al., 2017) and the US standard atmosphere from Anderson et al. (1986), in which the mean background concentrations of CO₂ and CH₄ are scaled to 420 and 1.81 ppm, respectively. The simulation of the radiance vector builds on 100 layers of 400 m thickness. Since we observe sky-scattered sunlight in shallow viewing elevation angles, we use an empirical adjustment by Kasten and Young (1989) for SZAs above 70° and VEAs below 20° to correct for the sphericity of the atmosphere.

Foote et al. (2021) show the importance of scene-specific unit enhancement spectra in a matched filter retrieval for significant parameters of airborne instruments, i.e., the atmospheric water column, surface elevation, solar zenith angle, and sensor altitude. We adapt their technique and compute the unit absorption spectra specifically for each observation period. Specific quantities are the mean SAA and VAA during observation as well as the aerosol optical thickness and scattering phase function. We assign an SZA and SAA to each observation from astronomical calculations based on the observation location and time and VEA and VAA based on geometrical calculations (Sect. 3.1). Aerosol information is taken from the closest AERONET station in Karlsruhe, Germany. We use the daily mean of the aerosol optical thickness (AOT) at 870 nm and the Ångström exponent (AE) to calculate the AOT at 2000 nm. Aerosol concentrations are assumed to be horizontally homogeneous and vertically decreasing exponentially with the atmospheric scale height. We use the asymmetry parameter of AERONET at 1020 nm, which is the largest wavelength available, to compute the Henyey–Greenstein phase function. Furthermore, we compute a look-up table of unit absorption spectra for a set of viewing geometries, using a specific unit absorption spectrum for each line in an image since the VEA covers a 16° range. Thus, our unit absorption spectrum is specific to the changing geometry within a single observation. The look-up table contains unit absorption spectra for five SZAs (10, 30, 50, 70°) and seven VEAs (2, 5, 8, 11, 14, 17, 20°). We interpolate a specific unit absorption spectrum for each line from the look-up table with the image mean SZA and the lines' VEA. Figure 4 shows CO₂ unit enhancement spectra at an SZA of 50° (right panel) at an VEA of 8°. Both plots show the same expected behavior from the observation geometry. The longer the light path in the atmosphere, the more light is absorbed by background CO₂. Thus, an additional enhancement has a smaller effect on the ASR, and the absolute values of the unit absorption spectrum decrease.

Figure 4 CO₂ unit absorption spectra for the HySpex SWIR-384 camera. Panel (a) shows the dependency on the solar zenith angle (SZA) and panel (b) the dependency on the viewing elevation angle (VEA).

4 Estimating emission rates

Various methods have been developed to estimate emission rates from images of CO₂ and CH₄ enhancement plumes of localized sources. CH₄ has received more attention recently, but the methods work mostly analogously for both gases. Gaussian plume models (GPMs) are widely used to estimate emission rates by fitting the model to snapshot observations of top-down observed plume images (Bovensmann et al., 2010; Krings et al., 2011; Rayner et al., 2014; Matheou and Bowman, 2016; Nassar et al., 2017; Schwandner et al., 2017). GPMs simulate the plume spread from point sources along a horizontal propagation direction. Typically, these models utilize stability classes to parameterize the turbulent dispersion properties. Fitting a model to the image can theoretically account for plume mass below the detection limit, which is an advantage over mass-balance methods. The GPM is only valid for an ensemble of plumes due to the stochastic nature of turbulence (Varon et al., 2018; Jongaramrungruang et al., 2019). Since we use temporal averages of plume measurements, we average over short-term turbulent structures and therefore meet the requirements for applying the GPM (Fig. 9). As we also observe the vertical plume rise with our observations, we use the bent-over Gaussian plume model by Janicke and Janicke (2001), which also accounts for the plume rise driven by buoyancy and initial vertical velocity after the release (Fig. 5). This allows us to simulate an hourly averaged observation of our camera from a set of parameters. The model input contains ambient conditions and source parameters, including the emission rate. We use the model to simulate an ensemble of observations and fit the simulation to the observations to estimate the emission rate. Section 4.1 describes the model for three-dimensional, bent-over Gaussian plume shapes, and Sect. 4.2 explains how we simulate an observation from the model. Section 4.3 describes the a priori data used for the emission estimation and Sect. 4.4 how we find emission estimates and uncertainties from the inversion.

Figure 5 Example for the output of the Gaussian plume model IBJpluris with a chimney of 180 m height. The example input parameters are the ambient wind velocity (7.0 m s⁻¹), ambient temperature (27 °C), and relative humidity (40 %) at 200 m height. Furthermore, the initial velocity (13.4 m s⁻¹), temperature (63.0 °C), and initial concentration (189 g m⁻³) of the exhaust gas are given. IBJpluris calculates the central plume axis (blue) and the plume boundaries (orange). Output parameters like the concentration of carbon dioxide are provided along the central plume axis.

4.1 Gaussian plume model

Since we observe the plumes in a horizontal viewing geometry, the model needs to account for the bent-over plume shape (Fig. 5). We use the plume rise model of Janicke and Janicke (2001) to calculate the plume properties along the central plume travel axis. Figure 5 shows the simulated plume shape for an ensemble of initial conditions. The properties are given as a set of discrete points P, which contain the spatial coordinates $x,y,z$; the plume radius R; the distance along the plume axis s; and the mass concentration c of a carried quantity, e.g., CO₂. The model requires ambient wind velocity u_a, temperature T_a, pressure p_a, and relative humidity RH_a as well as exhaust gas initial velocity u_e, temperature T_e, and concentration c₀ as input parameters (see Sect. 4.3). The total gas enhancement in the plume depends linearly on the emission rate Q. The concentration c₀ in the plume right above the chimney follows from

$$\begin{array}{}\text{(18)}& c_{\mathrm{0}}={\displaystyle \frac{Q}{\dot{V}}},\end{array}$$

where $\dot{V}$ is the air volume flux from the chimney. The concentration c₀ holds for a homogeneous plume segment of cylindrical form and radius R. We transfer each plume segment from a cylindrical concentration profile to the Gaussian profile,

$$\begin{array}{}\text{(19)}& c\left(r\right)=c^{\ast }\exp \left(-{\displaystyle \frac{r^{\mathrm{2}}}{b^{\mathrm{2}}}}\right),\end{array}$$

where c^∗ is the core concentration, r the distance to the central axis, and b the plume width. For each segment of the plume, c^∗ follows from the conservation of mass compared to a cylindrical plume segment with concentration c from

$$\begin{array}{}\text{(20)}& {\displaystyle }c\cdot \mathit{\pi }{R}^{\mathrm{2}}\mathrm{\Delta }s& {\displaystyle }=\underset{s_{\mathrm{0}}}{\overset{s_{\mathrm{1}}}{\int }}\mathrm{d}s\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\varphi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}r\mathrm{d}rc\left(r\right)\text{(21)}& {\displaystyle }& {\displaystyle }=\mathrm{\Delta }s\cdot \mathrm{2}\mathit{\pi }\cdot {\displaystyle \frac{c^{\ast }{b}^{\mathrm{2}}}{\mathrm{2}}}\text{(22)}& {\displaystyle }\Rightarrow c^{\ast }& {\displaystyle }=c\cdot {\displaystyle \frac{R^{\mathrm{2}}}{b^{\mathrm{2}}}},\end{array}$$

where s is the distance along the plume axis, i.e., Δs the plume cross-section segment thickness.

We create a three-dimensional domain around the central plume axis as our model domain, which covers at least a 4b radius around each point on the plume axis. The spatial resolution of the domain cells is a third of the HySpex pixel width (approximately 1 m). For each domain grid cell in the vicinity of the central axis, we find its distance r to the plume axis and the mass concentration c(r) at this distance from Eq. (19). We introduce two more parameters to Eq. (19), which are later used as fitting parameters for the inversion (Sect. 4.4). The parameters are a scaling factor k_c of the concentration c^∗ and a scaling factor k_b of the plume width b. We introduce them in a way that the total mass of the plume depends linearly on k_c and is independent of k_b. Thus, the concentration scaling represents the source strength, while the width scaling accounts for the turbulent diffusion during the time of observation and needs to be fitted (Carhart and Policastro, 1991). The total mass M_s in each slice of the plume is given by

$$\begin{array}{}\text{(23)}& {\displaystyle }{M}_s& {\displaystyle }=\underset{s_{\mathrm{0}}}{\overset{s_{\mathrm{1}}}{\int }}\mathrm{d}s\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\varphi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}r\mathrm{d}r{\displaystyle \frac{k_c}{(k_b{)}^{\mathrm{2}}}}{c}_s^{\ast }\exp \left(-{\displaystyle \frac{r^{\mathrm{2}}}{(k_b{b}_s{)}^{\mathrm{2}}}}\right)\text{(24)}& {\displaystyle }& {\displaystyle }=\mathit{\pi }\cdot \mathrm{\Delta }s\cdot b_s^{\mathrm{2}}\cdot k_c{c}_s^{\ast },\end{array}$$

where b_s and $c_s^{\ast }$ are the radial width and core concentration of the segment from s₀ to s₁, respectively. Thus, the parameter k_b scales the plume width without changing its mass, while k_c scales the mass and in the plume, which is linearly related to the source emission rate Q.

4.2 Observation forward model

We simulate a plume observation from the three-dimensional Gaussian model output by projecting the plume on our plane of observation and aggregating the mass according to the image pixels. A plume cell at ($x,y,z$) is projected on the observation plane that is perpendicular to the viewing direction. Figure 6 shows a sketch which explains the projection. The observation angle ϕ is the angle between the plume travel direction and the viewing direction, the plume cell angle $\mathit{\theta }=\arctan (x/y)$ is the angle between the location of a plume parcel location and the travel direction, and the projection angle β is the angle between the plume parcel and the observation plane. These three angles add up to 90°. Thus, given the cell location ($x,y,z$) and the observation angle ϕ, the projection angle β is given by $\mathit{\beta }=\mathrm{90}\mathit{°}-\mathit{\varphi }-\mathit{\theta }$. The observation angle is defined between −180 and +180°, where a negative angle denotes a plume moving to the left, a positive angle a plume moving to the right, and a zero angle a plume moving straight away from camera. The projection of the cell location ($x,y,z$) to the observation plane ($x^{\prime },z^{\prime }$) is given by

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{x}^{\prime }& =\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}\cdot \sin (\mathit{\varphi }+\arctan (x/y\left)\right),\\ z^{\prime }& =z,\end{array}\end{array}$$

where ϕ is known from the measurement geometry and the ambient wind direction, while ($x,y,z$) denotes the plume grid cells. We calculate a CO₂ mass distribution $m(x^{\prime },z^{\prime })$ in the observation plane by multiplying each concentration by the cell volume and projecting it according to Eq. (25). The simulated image is then given by aggregating the mass points m_jk in the pixels of the observation and converting them to a column enhancement ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{jk}$ using

$$\begin{array}{}\text{(26)}& {\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{jk}=m_{jk}\cdot {\displaystyle \frac{{\mathit{\nu }}_{{\mathrm{CO}}_{\mathrm{2}}}}{A_{jk}}},\end{array}$$

where ${\mathit{\nu }}_{{\mathrm{CO}}_{\mathrm{2}}}\approx \mathrm{0.546}$ m³ kg⁻¹ is the specific volume of CO₂ at normal conditions, and A_jk is the pixel area.

Figure 6 Sketch of the projection of a plume cell at ($x,y,z$) to the observation plane ($x^{\prime },z^{\prime }$). The gray shading shows a conceptual top-down view on a horizontal plume cross-section – the triangular shape is chosen for simplicity and does not represent the actual plume shape. The camera viewing direction is marked by the green line, such that the angle ϕ is the angle between the viewing direction and the y axis. The orange line points to an arbitrary plume cell at ($x,y,z$), with a plume cell angle of $\mathit{\theta }=\arctan (x/y)$. The blue line marks the projection plane which is perpendicular to the viewing direction. Thus, the projection angle β is given by $\mathit{\beta }=\mathrm{90}\mathit{°}-\mathit{\varphi }-\mathit{\theta }$. Using the distance $\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}$ of the plume cell to the origin, the projected length coordinate follows from elementary geometry to $x^{\prime }=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}\cdot \sin (\mathit{\varphi }+\arctan (x/y\left)\right)$.

We choose all input parameters of the model to be constant for an observation except for the ambient wind speed (Sect. 4.3). Thus, a simulated observation has a total of four independent parameters, which are ambient wind speed u_a, observation angle ϕ, emission scaling k_c, and the plume width scaling k_b.

4.3 A priori data

Plume formation in the atmosphere depends on ambient meteorological conditions like temperature, pressure, humidity, and wind velocity. The exhaust gas temperature, initial velocity, and source shape and diameter also affect the plume formation. The operator of Grosskraftwerk Mannheim (GKM) provided us with their operational data for all days on which we have been observing the GKM exhaust plume. The power plant operates four different sections, namely units 6, 7, 8, and 9. The stack height of the power plant is 180 m for unit 9 and 200 m for units 6 to 8. Their data contain, for each unit, the coal consumption, the exhaust volume flow rate at the chimney top, the chimney tip diameter, and the exhaust gas temperature. Unit 9 is the youngest facility, having been in operation since 2015. It contributes a maximum gross power production of 911 MW to the power plant total of 2146 MW. GKM also provided us with the instantaneous gross power production of unit 9. Furthermore, the GKM dataset includes 10 m meteorological data at the power plant location. These data consist of wind velocity and direction, temperature, pressure, and relative humidity. Both the meteorological and operational datasets have a resolution of 1 min.

We calculate the expected emission rate of CO₂ from power production using the official German emission factor for hard coal of 93.1 t CO₂ TJ⁻¹ (Sandau et al., 2021). Unit 9 is a modern unit with an efficiency of 46.4 % (Grosskraftwerk, 2015), giving an emission factor of 722 g CO₂ kWh⁻¹. We calculate the power plant's mean emission factor from its annual reported power production and the CO₂ emissions reported to the European Pollutant Release and Transfer Register (E-PRTR; Table S1). The total power plant has an average emission factor of 955 g CO₂ kWh⁻¹, leading to an estimated average emission factor of 1127 g CO₂ kWh⁻¹ for the older units 6, 7, and 8, for which efficiencies are unavailable. These factors allow for estimating emission rates from the instantaneous power production. Thus, the GKM operator provided us with all the necessary source input parameters for the Gaussian plume model. On 13 May 2022, unit 9 of GKM was malfunctioning, and the power production was shifted to unit 6. The Fraunhofer-Institut für Solare Energiesysteme (ISE) allocates the power production of German power plants with 15 min resolution. We use their data to estimate the power production of unit 6 on this day. The ambient input parameters were partly provided by the meteorological data ($p_{\mathrm{a}},T_{\mathrm{a}},{\mathrm{RH}}_{\mathrm{a}}$) of GKM and partly by our wind lidar (u_a,ϕ). We measure wind speed and direction in six height levels, the two uppermost being 150 and 200 m above ground level, from which we derive the a priori wind data at plume height. A 20 min running mean removes high-frequency fluctuations from the lidar data due to turbulence within a laser rotation. We average all input parameters over the same period used for the HySpex measurements, i.e., approximately 1 h. Typical wind speeds during observation are 5 m s⁻¹, for which we find hourly standard deviations of 0.9 m s⁻¹. Depending on meteorological stability during observation, the uncertainty of the wind direction varies between 18 and 44°. One exception to the a priori wind data is the observation on 8 September 2021. Since we do not have lidar data for this day, we scaled the 10 m winds of the GKM data to shaft height using data from the Copernicus Atmospheric Monitoring Service (CAMS) reanalysis (Inness et al., 2019). Figure S1 shows both the scaled wind speed and typical lidar observations. All input parameters for the Gaussian plume model are listed in Table S2.

4.4 Inverse estimate

As described in Sect. 4.2, there are four parameters that we can vary to fit the simulated image to the observed one. These are the plume width scaling k_b, the emission scaling k_c, the plume velocity u_a, and the observation angle ϕ. Despite wind lidar measurements being available, we use u_a and ϕ as free parameters since the wind lidar measurements are performed at kilometer distance from the source, and typical plume heights are above the lidar top height of 200 m.

We use a brute force method to scan over a space of parameter sets and compare the simulated observation with the observation. The plume mask is defined as the largest continuous patch of enhancements above twice the noise level σ_j given in Eq. (11). Figure 7a shows an example of an observed plume from GKM. Measurements in the same frames as the chimney are excluded since the comparatively high brightness of the chimney affects the other spectra in the frame. Figure 7b shows an example of a simulated plume. The translucent area is below the 2σ_j noise level, while the colored area is above 2σ_j. Enhancements above 2σ_j noise level are well above the detection limit; thus, the matched filter retrieval should detect them. We use the union of the plume mask and all pixels with an enhancement above 2σ_j of the noise level as our fit mask. Including background pixels in the fit mask keeps the fit close to the observation even outside the observed plume mask. While computationally costly, the brute force method provides further insights into the parameter space. For each parameter set, we calculate the reduced chi-square ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ by

$$\begin{array}{}\text{(27)}& {\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N-\mathrm{4}}}\sum _{jk}^{\mathrm{fitmask}}{\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{jk}(k_b,k_c,u_{\mathrm{a}},\mathit{\varphi })-{\mathit{\alpha }}_{jk}}{{\mathit{\sigma }}_j}}\right)}^{\mathrm{2}},\end{array}$$

where N is the number of pixels in the fit mask, ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{jk}(k_b,k_c,u_{\mathrm{a}},\mathit{\varphi })$ is the simulated column enhancement in pixel jk for the parameter set $(k_b,k_c,u_{\mathrm{a}},\mathit{\varphi })$, α_jk is the observed column enhancement in pixel jk, and σ_j is the uncertainty of the observed column enhancement for the pixels in line j. Figure 7c shows the residual of the observation and the best fit simulation.

Figure 7 Observed (a) and simulated (b) CO₂ plume from GKM Mannheim. The simulations correspond to the optimal parameter set. Residuals are shown in panel (c). The black contour in each panel marks the fitmask, which is the union of the plumemask from panel (a) and all pixels with a simulated enhancement larger than 2⋅σ_j (bright colors in panel (b)). The part above the chimney is excluded since the chimney is so bright that it affects the retrievals for the respective entire frames. The example is the observation from 26 March 2022, 15:56–17:36 UTC.

Figure 8 Cross-sections through the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ hypersurface (color coded) for the four-dimensional parameter space: emission scaling factor k_c, plume width factor k_b, wind speed u_a, and observation angle ϕ. The contour lines mark where the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ increases by 1, 2, and 3. The blue dot marks the minimal ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$. Note that all cross-sections involving the observation angle are symmetric around −90°, since a plume moving away from or toward us looks identical. The example is the observation from 26 March 2022, 15:56–17:36 UTC.

Since we scan the parameter space in a brute force manner, we calculate the reduced ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ for each parameter set in a wide parameter range. We use the optimal parameter set (${\widehat{k}}_b,{\widehat{k}}_c,{\widehat{u}}_{\mathrm{a}},\widehat{\mathit{\varphi }}$) with the lowest ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ for our emission estimate. Figure 8 shows an illustrative ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ hypersurface of the four-dimensional parameter space along the cross-sections through the optimal parameter set for the plume on 26 March 2022, retrieved from the averaging period between 15:56–17:36 UTC. The ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ surfaces are smooth, which indicates that ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ is a continuous function of the parameters. There is a unique minimum on each hypersurface, marked by a blue dot, around which the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ increases monotonically. For purely statistical errors, an increase of 1 in ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ corresponds to a mean deviation of 1 standard deviation (1σ) between the simulated and observed image due to parameter changes (Bevington et al., 1993). Thus, we consider all parameter sets with

$$\begin{array}{}\text{(28)}& {\displaystyle }{\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}<min\left({\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\right)+\mathrm{1}\end{array}$$

for the uncertainty of the emission estimate. Fixing (${\widehat{k}}_b,{\widehat{u}}_{\mathrm{a}},\widehat{\mathit{\varphi }}$) to the optimal parameter set, we can estimate the uncertainty of the emission estimate by varying k_c within the range of Eq. (28). The minimum and maximum emission estimates are used as uncertainty ranges. This estimate neglects systematic errors between observation and simulation, which are challenging to account for. Therefore, we use the criterion as our best approximation of the uncertainty range.

The shape of the well in the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ hypersurface reveals correlations between the parameters. A circular well indicates that the parameters are uncorrelated, such as for the well spanned by the wind speed parameter u_a and the plume width parameter k_b. Flat wells suggest that the observations do not constrain the parameter strongly since the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ does not change much with the parameter. Figure 8 shows that the observation angle is the least constrained parameter and that there is a correlation between the plume width parameter k_b and emission scaling factor k_c. This correlation is expected since the plume width scaling distributes the total mass over a larger area, while the emission scaling can increase the overall mass to match the observation again. Observations under favorable conditions constrain the ambiguity well since they provide enough information on plume width. If the observation is made under challenging conditions, such as clouds or low emission rates, the plume width is not well constrained, which leads to a flat slant well in the plume width–emission scaling (k_b–k_c) plane. Some observations show a correlation between wind speed and wind direction (e.g., Fig. S5), which is explained by an ambiguity in the observed plume shape for horizontally viewing observers. A plume traveling perpendicular to the viewing angle at a slow wind speed will look the same as a plume traveling at a higher wind speed in an acute observation angle. In theory, the emission scaling is unaffected by this ambiguity. This becomes clear using a simple mass-balance argument. The observation gives the plume mass, while the travel time of the plume is approximately given by plume length divided by wind speed. A geometric observation factor will act on both speed and direction similarly, thus canceling out in the travel time and emission estimate.

The ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ surfaces alongside a plot of the best fit for each observation are given in the Supplement. We sample the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ surfaces with step sizes of 0.2–0.3 m s⁻¹ for the ambient wind speed, 5° for the wind direction, 0.1 for k_c and 0.1–0.4 for k_b. We sample the k_c dimension a final time with 0.01 step size at the optimal parameters (${\widehat{k}}_b,{\widehat{u}}_{\mathrm{a}},\widehat{\mathit{\varphi }}$) to improve the emission estimate given the other three parameters.

5 Results

We present CO₂ plume images from GKM for 5 d in 2021 and 2022. On each day, we measured hyperspectral images for several hours, as listed in Table 1. Figure 9 shows the total of 11 retrieved plumes identified by the iterative matched filter algorithm described in Sect. 3.2. In every observation, the plume is clearly visible and can be attributed to a power plant chimney. For each of the plume images, we estimate the emission rates according to Sect. 4.4.

Figure 9 Plume images for the periods listed in Table 2. The yellow to red color bar is shared over all plumes and shows the atmospheric CO₂ enhancement (in ppm m). The blue (dark) to white (bright) color in the background shows the spectrum saturation. The bright rectangular shapes in the lower part of the images are the power plant chimneys.

The observations on 8 September 2021 were made under favorable conditions with clear skies and a priori observation angles of $-\mathrm{53}\pm \mathrm{13}$°. The sky, though, became more heterogeneous in the afternoon, causing the retrieval noise to increase (Fig. S3). On 23 March 2022, there was significant condensation in the early stages of the plume, and an acute observation angle 33 ± 28° poses further challenges for the retrieval of the plume images collected on that day. On 26 March 2022, we again had favorable conditions allowing for measuring two plume images. The power plant ramped up its power production during our observation period, so we can compare estimates over a range from 223 to 455 t h⁻¹. A Sahara dust event increased the aerosol load on 28 March 2022, which caused the sky to be hazy and bright in the shortwave infrared. On 13 May 2022, we observed power plant unit 6 instead of unit 9.

The plume images collected under unfavorable conditions on 23 March and 13 May 2022 are particularly well-suited to illustrate the limitations of the method. The acute observation angles on 23 March 2022 reduce the apparent plume size due to the unfavorable projection. Furthermore, condensation in the early stages of the plume removes a considerable part of the plume. Thus, we consider these observations challenging for the retrieval. The plume images on 13 May 2022 are taken under comparably high aerosol load, with AOT in the range of 0.065 (at 2000 nm) and favorable wind conditions of approximately 5.5 m s⁻¹ at a 62° observation angle. Yet, the plume images show unexpected shapes that are not reproducible by the Gaussian plume model (Fig. S2). We observe two enhancement patches, one above and below the expected plume shape. One stripe-like patch is around 290 m height, while another patch is below the chimney. A possible reason for these patterns is spectral artifacts due to background heterogeneity, e.g., thin clouds. Irrespective of the origin of these artifacts, the Gaussian model is not capable of representing these patterns, and thus, our approach is ill-suited for emissions estimates from these images.

Figure 10 Correlation of retrieved emissions and expected emissions. Every color represents a different observation day, while each symbol represents a different time interval. Filled symbols mark observations under favorable conditions, while open symbols mark observations under unfavorable conditions. The dashed black line marks the 1 : 1 line.

Table 2 Column one lists the observation time period and column two the expected emissions E_exp (Sect. 4.3). The retrieved emissions E_ret (Sect. 4.4) are listed in column three, with the uncertainty range in parentheses. Columns four to seven give the optimal inversion parameters ${\widehat{k}}_c$, ${\widehat{k}}_b$, $\widehat{\mathit{\varphi }}$, and ${\widehat{u}}_{\mathrm{a}}$. Note that ${\widehat{k}}_c$ is the relative scaling between retrieved and expected emissions; i.e., it represents the relative deviation. The last column shows the minimum ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ for each observation.

For each of the observations, we calculate the CO₂ emission rates as described in Sect. 4. Since we have precise knowledge of the power plant's power production from the company itself, we can compare our retrieved emissions to the correlative bottom-up calculations. These expected emissions are considered errorless in comparison. Figure 10 shows how our retrieved emissions compare to the expected ones calculated as described in Sect. 4.3. Table 2 lists the retrieved emissions with the optimal parameter sets for each observation.

We consider 7 out of the 11 plume images to be taken under favorable conditions. For these observations on 8 September 2021, 26 March 2022, and 28 March 2022, we find an overall reasonable agreement between retrieved and expected emissions. The retrieved emissions average 84 % of the expected emissions and have a mean relative uncertainty of 24 %. Five of the seven observations agree with the expected emissions within their uncertainties. Notably, the retrieved emissions agree well with the variability of the expected emissions from 223 to 455 t h⁻¹ on 26 March 2022, indicating that the method can observe diurnal changes in emission rates. The plume on 28 March 2022 was observed during a high aerosol load and small wind speeds (${\widehat{u}}_{\mathrm{a}}=\mathrm{1.4}$ m s⁻¹), which gives a preliminary lower limit for the wind speed necessary for the method to work. Observations on 8 September 2021 between 14:24 and 16:35 UTC agree with the expected emissions within the uncertainty range, while observations between 12:13 and 14:23 UTC underestimate the expected emissions. The estimated observation angle $\widehat{\mathit{\varphi }}$ is between −35 and −30° for these observations. This is significantly steeper than the a priori observation angle of −53 ± 13°, which was derived from ERA5 data (Sect. 4.3) since there are no wind lidar data on 8 September 2021. Observations on 23 March 2022 indicate that acute observation angles may cause emission underestimation, which might also apply to the observations on 8 September 2021. Potential sources of systematic errors in the retrievals are background heterogeneity of the scene, CO₂ features in the image region of the reference spectrum, or assumptions in the unit absorption spectrum calculations like aerosol content.

Observations on 23 March 2022 have been made under challenging conditions, as described above. The measurement between 14:51 and 16:13 UTC underestimates the expected emissions significantly by 61 % to 77 %. We find a retrieved observation angle $\widehat{\mathit{\varphi }}$ of 15° for this period, which agrees reasonably well with the a priori observation angle of 33 ± 30°. The measurement from 15:35 to 17:36 UTC agrees with the expected emissions within the uncertainty, but its observation angle $\widehat{\mathit{\varphi }}$ is 65° is inconsistent with the lidar observation of 33 ± 25°. Thus, acute observation angles and plume condensation are considered factors that limit the method's applicability and may be used as filter criteria. The observations on 13 May 2022 overestimate the power plant emissions. The enhancements outside the expected plume cause an increase in the width scaling factor in the GPM inversion since the simulated enhancements need to spread out (Fig. S22). While the width scaling factors range typically between 1.2 and 2.0, we find width scaling factors of 7.0 for these observations. The additional mass is attributed to the power plant emissions, leading to overestimation of the emission scaling factor. For the observation on 13 May 2022, such biased observations can be identified by non-compliance with a Gaussian plume shape and excluded from the emission estimate. However, in other cases, the effect might be too small to be identified visually but still large enough to propagate in the emission estimate.

6 Conclusions

We report on a proof of concept for estimating CO₂ emission rates for a coal-fired power plant using a ground-based hyperspectral camera. The power plant is located in Mannheim, Germany, and has yearly CO₂ emissions of more than 4.9 Mt CO₂ yr⁻¹. We demonstrate our capability to reliably image CO₂ plumes from individual chimneys with 11 observations over 5 d in 2021 and 2022. The camera observes sky-scattered sunlight in the shortwave infrared spectral range above the chimneys of the power plant. We use an iterative matched filter algorithm to retrieve the CO₂ enhancements from the observed spectra. Scattering on molecules and aerosols is inefficient in the infrared spectral range. Thus, the observed signal is small, and we need to average over more than 50 min to reach a sufficient signal-to-noise ratio in typical cases. Averaging over such a period typically causes the observed plume to be of Gaussian shape in good approximation. Therefore, we estimate emissions by fitting a Gaussian plume model to the observed plume. The forward model is based on the plume rise model by Janicke and Janicke (2001). We estimate the ambient wind velocity and direction, the plume width, and the source emissions by minimizing the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ between the observed and simulated plume. Therefore, we sample the ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$ space using a brute force approach, which, for all cases, reveals an unambiguous minimum in ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}$. For validation, we calculate the expected emissions based on the power plant's power production during the observation.

Favorable observation conditions are homogeneous skies, stable wind speeds, and a wind direction perpendicular to the viewing direction. We present seven observations taken on 3 d with ambient conditions matching such favorable conditions. For these seven observations, the estimated emission rates average 84 % of the expected emissions with a mean relative uncertainty of 24 %; thus, they agree reasonably well with the expected emissions. Observations indicate that we can follow the diurnal trend of the power plant emissions under such conditions and that our technique works in wind speeds down to ∼ 1.4 m s⁻¹.

We also present plumes that show enhancement artifacts, plume condensation, and unfavorable wind conditions to demonstrate the limitations of the technique. Obvious non-compliance with the Gaussian plume shape causes inaccurate emission estimates, but these cases are easily identified in the image and neglected for further analysis. Furthermore, we find that plume condensation and steep observation angles pose challenges for the technique. More observations are needed to quantify the impact of these effects on the emission estimates and to develop quantitative and suitable filter criteria when measuring under non-favorable conditions.

Our spectral imaging technique adds to the pool of tools to verify CO₂ emission rates of localized sources. In that context, our ground-based setup allows for monitoring individual sources over prolonged periods, which, for example, is complementary to the snapshot images provided by satellites. We envision that further development of our technique can provide independent data for emission inventories and can be used to verify bottom-up emission estimates. Our instruments fit in a car, and the ground-based observation geometry enables us to choose the targets flexibly. Potential future targets are less well-known anthropogenic sources such as facilities in the chemical, metallurgy, or cement industry and natural CO₂ sources like volcanoes.