Development of a small unmanned aircraft system to derive CO₂ emissions of anthropogenic point sources

A reduction of the anthropogenic emissions of CO₂ (carbon dioxide) is necessary to stop or slow down man-made climate change. To verify mitigation strategies, a global monitoring system such as the envisaged European Copernicus anthropogenic CO₂ monitoring mission (CO2M) is required. Those satellite data are going to be complemented and validated with airborne measurements. Unmanned aerial vehicle (UAV)-based measurements can provide a cost-effective way to contribute to these activities. Here, we present the development of an sUAS (small unmanned aircraft system) to quantify the CO₂ emissions of a nearby point source from its downwind mass flux without the need for any ancillary data. Specifically, CO₂ is measured by an NDIR (non-dispersive infrared) detector, and the wind speed and direction are measured with a 2-D ultrasonic acoustic resonance anemometer. By means of laboratory measurements and an in-flight validation at the ICOS (Integrated Carbon Observation System) atmospheric station Steinkimmen (STE) near Bremen, Germany, we estimate that the individual CO₂ measurements have a precision of 3 ppm and that CO₂ enhancements can be determined with an accuracy of 1.3 % or 0.9 ppm, whichever is larger. We introduce an anemometer calibration method to minimize the effect of rotor downwash on the wind measurements. This method derives the fit parameters of a linear calibration model accounting for scaling, rotation, and a potential constant bias. For this purpose, it analyzes wind measurements taken while following a suitable flight pattern and assuming stationary wind conditions. From the calibration and validation experiments, we estimate the single measurement precision of the horizontal wind speed to be 0.40 m s⁻¹ and the accuracy to be 0.33 m s⁻¹. By means of two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten about 40 km west of Bremen, Germany, we demonstrate how the measurements of elevated CO₂ concentrations can be used to infer mass fluxes of atmospheric CO₂ related to the emissions of the facility.

1 Introduction

CO₂ (carbon dioxide) emissions are the primary cause of man-made climate change, and in order to limit this, a reduction of emissions is necessary (IPCC, 2013). Large parts of the anthropogenic CO₂ emissions originate from point sources such as coal- or gas-fired power plants and observing systems are needed to verify mitigation strategies through independent measurements (Pinty et al., 2017). On a global level, it is planned to monitor CO₂ emissions remotely by means of satellites (e.g., Nassar et al., 2017; Reuter et al., 2019) such as the envisaged European Copernicus anthropogenic CO₂ monitoring mission (CO2M; Janssens-Maenhout et al., 2020) building upon the heritage of CarbonSat (Bovensmann et al., 2010; Buchwitz et al., 2013), which was a candidate for European Space Agency (ESA)'s Earth Explorer 8 mission. At the regional level, this can be achieved by a dense network of ground-based observations or by airplane-based measurements, which can also be used for smaller point sources and for the validation of the satellite data (e.g., Krings et al., 2011, 2018; Carotenuto et al., 2018). Unmanned aerial vehicle (UAV)-based measurements can complement these activities, reduce flight costs, and enable measurements under conditions not suitable for satellite measurements, e.g., in cloud-contaminated scenes, during nighttime, or at facilities with emissions below the satellite's detection limit.

The quantification of emissions with the mass balance approach using atmospheric measurements downwind of the source requires knowledge of the atmospheric concentration of CO₂, the wind speed, and the density of air. As done, e.g., by Krings et al. (2018), the CO₂ enhancement from the source can be computed from the measurements of a cross-section of its plume by subtracting the CO₂ background concentration derived from measurements outside the plume. Additionally, measurements of the wind speed and the density of air are needed to compute the CO₂ cross-sectional flux by integrating over the cross-sectional flux density.

Earlier studies showed that UAVs are suitable to serve as platforms for in situ CO₂ sensors (e.g., Berman et al., 2012; Khan et al., 2012; Kunz et al., 2018; Allen et al., 2019; Chiba et al., 2019), anemometers (e.g., Palomaki et al., 2017; Hollenbeck et al., 2018; Shimura et al., 2018; Barbieri et al., 2019), and/or sensors for meteorological parameters (e.g., Barbieri et al., 2019). Our goal was to create an sUAS (small unmanned aircraft system) which is capable of determining atmospheric CO₂ mass fluxes from its own sensor data independently from external data sources. The system was required to be reliable and affordable, and to be realized within a relatively short development phase. For this reason, we decided to use only commercially available and mature components for the UAV as well as the sensor equipment.

In the following section, we introduce the hardware setup and the instrumentation. In Sect. 3, we characterize the measurements of our in situ CO₂ sensor by comparing them with those from an accurate laboratory instrument, and in Sect. 4 we introduce a calibration method for the anemometer correcting, e.g., for the influence of rotor downwash. In Sect. 5, we quantify the performance of our atmospheric measurements with a focus on CO₂ and wind by analyzing two validation flights at the ICOS (Integrated Carbon Observation System) atmospheric station Steinkimmen (STE) near Bremen, Germany. In Sect. 6, we demonstrate how the sUAS can be used to measure elevated CO₂ concentrations from a nearby source and how these measurements can be used to infer the CO₂ mass flux related to the emissions of the source. Finally, we summarize and conclude our results.

2 Hardware setup and instrumentation

We use a DJI Matrice 210v2 (https://www.dji.com, last access: 2 June 2020) as UAV airborne platform, which weighs about 4.8 kg (including batteries) and can carry a maximum payload of 1.34 kg. An overview of the sUAS including UAV and instrumentation is shown in Fig. 1. The DJI Matrice 210v2 is powered by two LiPo (lithium polymer) batteries with a capacity of 7660 mAh at 22.8 V allowing for flight times typically in the range of 24 to 34 min, depending on payload. According to the technical specifications of the UAV, it can operate in temperatures between −20 and 50 ^∘C; its maximum flight altitude is reached at 3 km above mean sea level (with special propellers) or at 500 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (limited by its firmware); under optimal conditions, its radio system can operate over distances of up to 8 km; its maximum wind resistance is 12 m s⁻¹. The DJI Matrice 210v2 is often used for industrial applications and emergency services. Its top cover has four mounting threads designed to attach additional payloads and we use it to mount an aluminum base plate.

Figure 1 UAV and instrumentation. (a) Overview including UAV coordinate system as defined by us. (b) Zoom on base plate.

The Vaisala GMP343 in situ CO₂ sensor (https://www.vaisala.com, last access: 2 June 2020) is the primary measurement device aboard the UAV and the only part of the equipment not attached to the base plate but to the secondary gimbal connector of the UAV. The GMP343 has an NDIR (non-dispersive infrared) sensor that uses an electrically tunable Fabry–Pérot interferometer to switch back and forth between a CO₂ absorption band at around 4.26 µm and the neighboring continuum at around 3.9 µm serving as reference. It is specified to measure the mole fraction of atmospheric CO₂ with a single sounding precision of ±3 ppm (1σ). In order to achieve a fast response time, we operate the CO₂ sonde with an open measuring cell through which the ambient air flows.

The second most important sensor aboard the UAV is the FT Technologies FT205 anemometer (https://fttechnologies.com, last access: 2 June 2020). It is a lightweight 2-D ultrasonic acoustic resonance anemometer specifically designed for UAV applications. It is specified to deliver measurements of the horizontal wind speed with a precision of ±0.3 m s⁻¹ for wind speeds below 16 m s⁻¹. In order to reduce the influence of the rotor downwash on the wind measurements, we mounted the anemometer 36 cm above the rotor plane on a carbon fiber pole with a diameter of 7 mm.

Deriving the CO₂ mass flux from the atmospheric mole fraction and the wind speed requires knowledge of the molar air density, which we compute from the measured meteorological parameters. For this purpose, we use two Adafruit (https://www.adafruit.com, last access: 2 June 2020) breakout boards equipped with a Bosch BMP388 pressure sensor and a Sensirion SHT31-DIS temperature and humidity sensor. The BMP388 pressure sensor has a specified absolute accuracy of ±0.5 hPa (±0.08 hPa relative accuracy). The Sensirion SHT31-DIS temperature and relative humidity measurements have a specified typical accuracy of ±0.3 ^∘C and ±2 %, respectively. GNSS (Global Navigation Satellite System) position as well as attitude information is provided in flight directly by the UAV via DJI's OSDK API (onboard software development kit application programming interface).

With a maximum rate of 0.5 Hz, the CO₂ sensor offers the slowest readout rates of our sensor equipment, and we decided to synchronize all instrument readouts with those of the GMP343 in order to generate a convenient output data stream with the same time stamp for all sensors.

A Dronee Zoon v2 telemetry module (https://dronee.aero, last access: 2 June 2020) that allows transmission ranges up to 3 km establishes an independent bidirectional data link. It is used to transmit the most important measurement data (CO₂, wind, position, etc.) to the ground enabling, e.g., to optimize the flight pattern according to the wind measurements during flight.

All instruments are connected to a Raspberry Pi 3B+ via an RS-232 to USB converter (CO₂), an RS-485 to USB converter (wind), a transistor–transistor logic (TTL) UART to USB converter (UAV), USB cable (telemetry), and I²C-bus (pressure, temperature, and relative humidity). The Raspberry Pi as well as the entire science payload is powered by the batteries of the UAV.

Beside the science payload, the UAV carries also a DJI X4S camera for documentary purposes and for high-quality, gimbal-stabilized first-person view during flight.

The entire payload weighs about 1.2 kg and consumes 6 W, reducing the maximum flight time of 34 min to about 24.7 min. The reduction in flight time results mainly from the weight of the payload and only marginally from its power consumption. The individual contributions are listed in Table 1. All components of the sUAS have a total retail price of about EUR 15 000 including the UAV, which has the highest unit price followed by the CO₂ sensor.

Table 1 UAV and instrumentation: weight, power consumption, and contribution to flight time. The reduction in flight time results from the weight and the power consumption of the payload.

^a Average power consumption during flight including camera and for the total take-off weight of 6022 g. ^b Considering only the weight of the camera. ^c Including meteorological sensors, aluminum housing, and power converter.

3 CO₂ sensor characterization

Flux quantification of an individual source requires measurements of the difference between elevated concentrations minus undisturbed background concentrations. This means that a constant bias of the measurements is less critical than drifts during flight. Unfortunately, some NDIR gas sensors can experience such drifts, e.g., related to temperature (e.g., van Leeuwen et al., 2013; Shusterman et al., 2016; Kunz et al., 2018; Ouchi et al., 2019). The GMP343 corrects for changes in temperature, pressure, or humidity. For this purpose, it measures the temperature with an internal sensor and accepts user inputs for pressure and humidity. Shusterman et al. (2016) and van Leeuwen et al. (2013) found some residual dependencies of the sensor output to temperature, pressure, and humidity and proposed to derive an empirical correction which can be specific for a sensor unit.

We compare roughly 1 week of continuous measurements of our CO₂ sensor with simultaneous measurements of a highly precise (better than ±0.3 ppm without averaging as done here) ABB LGR-ICOS ultra-portable greenhouse gas analyzer (https://new.abb.com, last access: 2 June 2020) in the laboratory under conditions of varying ambient pressure, temperature, and humidity (Fig. 2).

Figure 2 (a) CO₂ concentration measured in the laboratory with the Vaisala GMP343 CO₂ sonde with and without linear correction (Eq. 1) as well as highly accurate reference measurements performed with an ABB LGR-ICOS ultra-portable greenhouse gas analyzer. (b) Difference between the Vaisala GMP343 CO₂ (with and without linear correction) and the reference measurements. Pale colors represent instantaneous differences and intense colors 1 h running averages. (c) Temperature measured with the Vaisala GMP343. (d) Relative humidity measured with a Sensirion SHT31-DIS sensor. (e) Pressure measured with a Bosch BMP388 sensor.

Air pressure is measured with the BMP388, temperature with the GMP343 internal, and humidity with the SHT31-DIS sensor. During the measurement period, the CO₂ concentration changes by more than 200 ppm, the air pressure varies by about 17 hPa, the relative humidity is in the range between about 35 % and 60 %, and the temperature changes by approximately 6 ^∘C. This means, except for CO₂, all parameters vary more than we expect for a typical measurement flight.

Throughout the entire comparison time, the CO₂ measurements are always close together with an average difference of −0.33 ppm. We compute the 1 h running average of the difference between the GMP343 and the LGR instrument. Deviations from zero can safely be assumed to not be dominated by instrumental noise or short-term fluctuations. The standard deviation of the running average amounts to 0.89 ppm, and we consider it an estimate for potential systematic CO₂ drifts during a measurement flight.

The individual soundings scatter around the running average with a standard deviation of 1.76 ppm. However, one of the first test flights with very variable CO₂ concentrations near a stack revealed that always two GMP343 measurements (sampling rate 0.5 Hz) are paired; i.e., they exhibit similar CO₂ values within some ppm, and Vaisala's technical support confirmed that truly independent measurements are only possible every 4 s. In order to account for this, we multiply the derived standard deviation by $\sqrt{\mathrm{2}}$ and get 2.49 ppm as an estimate for the noise of the individual soundings. This is an upper limit of the measurement noise because the multiplication by $\sqrt{\mathrm{2}}$ would imply that both measurements within a 4 s interval are identical, which is not the case.

In order to analyze the drifts of the 1 h running average in more detail, we searched for correlations with temperature (T in ^∘C), pressure (p in hPa), absolute humidity (h_a in g m⁻³), and the measured CO₂ concentration (GMP343 in ppm). Similar to the approach used by van Leeuwen et al. (2013), we fit a linear model with these parameters as explanatory variables and the CO₂ concentration of the LGR instrument as response variable and get the following potential correction (in ppm) for our GMP343 sensor:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }{\text{GMP343}}^{\prime }& {\displaystyle }=\mathrm{1.013}\text{GMP343}+\mathrm{0.434}T{\displaystyle }& {\displaystyle }-\mathrm{0.013}{h}_{\mathrm{a}}+\mathrm{0.021}p-\mathrm{38.952}.\end{array}$$

Note that the absolute humidity is computed from the relative humidity (h_r in %) by

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{h}_{\mathrm{a}}=h_{\mathrm{r}}/\mathrm{100}\text{SVD},\end{array}$$

wherein SVD is the saturated vapor density (in g m⁻³), which can be well approximated for temperatures between 0 and 40 ^∘C by the following equation (Nave, 2017):

$$\begin{array}{ll}\text{(3)}& {\displaystyle }\text{SVD}& {\displaystyle }\approx \mathrm{5.018}+\mathrm{0.32321}T+\mathrm{8.1847}\cdot {\mathrm{10}}^{-\mathrm{3}}{T}^{\mathrm{2}}{\displaystyle }& {\displaystyle }+\mathrm{3.1243}\cdot {\mathrm{10}}^{-\mathrm{4}}{T}^{\mathrm{3}}.\end{array}$$

In principle, Eq. (1) could be used as refinement of the internal corrections of our GMP343 sensor unit, and indeed, Fig. 2 suggests that it could reduce potential drifts to 0.45 ppm. However, in terms of expected variability of temperature, pressure, and humidity during a measurement flight, the variability of the correction is small (<1 ppm). The dominating part of the correction comes from the apparent underestimation of 1.3 % of the CO₂ variability as measured by the GMP343.

However, it can be expected that other factors (especially wind speed; see Sect. 4) add larger uncertainties to the flux estimates. Therefore, and because more tests would be needed to confirm the robustness of the derived linear model, we decided to use the unmodified, only internally corrected GMP343 measurements as read from the instrument.

4 Anemometer calibration

Operating an anemometer near an object that can disturb the local wind field (e.g., a tower) requires calibration to derive the wind speed of the undisturbed wind field. This is particularly the case for the operation aboard a UAV influencing the local wind field not only by its static parts but especially by the downwash of the rotors. We tried to reduce both effects by mounting the anemometer on a 7 mm diameter carbon fiber pole about 36 cm above the rotor plane. A larger distance to the rotor plane would of course reduce the influence of the rotors, but the flight behavior could suffer if the distance becomes too large. In addition, pitch and roll maneuvers of the UAV would result in larger relative velocities of the anemometer.

Let m^′ be the horizontal wind vector (in the x and y directions of the UAV coordinate system; see Fig. 1) as measured by the anemometer. We define the calibrated wind vector u^′ (also in UAV coordinates) by the following calibration function which allows for scaling and rotating the original measurements by matrix A and for adding a constant vector b.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathit{u}}^{\prime }=\mathbf{A}{\mathit{m}}^{\prime }+\mathit{b}\end{array}$$

Applying the rotation matrix R with γ being the azimuthal orientation (yaw) of the UAV,

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathbf{R}=\left(\begin{array}{cc}\cos \mathit{\gamma }& -\sin \mathit{\gamma }\\ \sin \mathit{\gamma }& \cos \mathit{\gamma }\end{array}\right),\end{array}$$

results in the calibrated wind at the anemometer but in geographic coordinates:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\mathit{u}={\mathbf{R}}^{-\mathrm{1}}(\mathbf{A}{\mathit{m}}^{\prime }+\mathit{b}).\end{array}$$

The wind at the anemometer is a superposition of the wind relative to ground w and the headwind because of the velocity v of the UAV ($\mathit{u}=\mathit{w}-\mathit{v}$), so that the wind relative to ground becomes

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\mathit{w}={\mathbf{R}}^{-\mathrm{1}}(\mathbf{A}{\mathit{m}}^{\prime }+\mathit{b})+\mathit{v}.\end{array}$$

A straightforward approach to derive the coefficients of A and b would be to make measurements in a wind tunnel at different wind speeds w and orientations γ of the UAV while continuously hovering at the same position (v=0). However, this procedure is impractical because it requires a relatively large wind tunnel and hovering in place becomes more difficult if GPS is not available.

Figure 3 Flight pattern of the anemometer calibration experiment. From top to bottom: location (a), distance to the (arbitrarily chosen) origin in the east direction (b), distance to the (arbitrarily chosen) origin in the north direction (c), ground speed (d), and yaw (e) of the UAV. Valid measurements are shown in dark blue; filtered measurements are in light blue. The filtering removes measurements that do not lie within the calibration pattern or with UAV accelerations larger than 0.5 m s⁻².

Consequently, we selected a different calibration strategy: under the assumption of a stationary wind field, we performed a calibration flight with a flight pattern that enables w and the coefficients of A and b to be derived simultaneously. The flight pattern has to include back and forth movements of the UAV in two ideally perpendicular directions and has to be repeated with varying velocities and orientations (see Fig. 3).

Specifically, we solve Eq. (7) for the uncalibrated anemometer measurement:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{m}}^{\prime }={\mathbf{A}}^{-\mathrm{1}}\left[\mathbf{R}(\mathit{w}-\mathit{v})-\mathit{b}\right],\end{array}$$

and perform a least squares fit to find the coefficients of A, b, and w. In order to best fulfill the stationarity assumption of the wind field, we performed the calibration flight on a day with relatively calm conditions and low gusts. Additionally, we use only measurements with a horizontal total acceleration of less than 0.5 m s⁻². The results for the free fit parameters,

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \mathbf{A}=\left(\begin{array}{cc}\mathrm{0.881}\pm \mathrm{0.010}& \mathrm{0.018}\pm \mathrm{0.010}\\ -\mathrm{0.021}\pm \mathrm{0.020}& \mathrm{0.832}\pm \mathrm{0.010}\end{array}\right),\\ & \mathit{b}=\left(\begin{array}{c}\mathrm{0.074}\pm \mathrm{0.040}\\ -\mathrm{0.058}\pm \mathrm{0.038}\end{array}\right),\mathit{v}=\left(\begin{array}{c}\mathrm{1.064}\pm \mathrm{0.041}\\ \mathrm{1.063}\pm \mathrm{0.041}\end{array}\right),\end{array}\end{array}$$

show that the calibration scales the measured wind speed in the x direction with 0.881 and 0.832 in the y direction. A rotation is basically not applied and there is only a very small constant offset correction below 0.1 m s⁻¹ needed.

Figure 4 Results of the anemometer calibration experiment. (a) X component of the uncalibrated wind measurement and roll (blue and red, left side of Eq. 8) as well as corresponding fits (green and orange, right side of Eq. 8). (b) Y component of the uncalibrated wind measurement and pitch as well as corresponding fits. (c) Wind speed relative to ground in the east direction for the calibrated (green) and uncalibrated (blue) anemometer, as well as that computed from the attitude of the UAV (orange). (d) Wind speed relative to ground in the north direction. Histograms of wind speed and direction for the calibrated (g) and uncalibrated (f) anemometer, as well as that computed from attitude (e).

As shown in Fig. 4a and b, the fitted wind measurements (right side of Eq. 8, green lines in Fig. 4a and b) agree well with the actually measured values (left side of Eq. 8, blue). The derived wind relative to ground w (Fig. 4c and d, green) scatters around the mean values of 1.06 m s⁻¹ in the east and west directions. We repeated the calibration experiment and obtained similar coefficients even though the conditions were less ideal with a slightly larger average wind speed of 2.40 m s⁻¹ and more gusts (not shown).

We estimate the scatter of the wind components in two different ways. (i) We compute the standard deviation of the wind components and consider it an upper limit of the scatter as it includes not only the noise of the measurements (including calibration) but also the actual variability of the wind. (ii) Under the assumption that the actual wind varies only smoothly, we estimate the noise of the measurements by the standard deviation of the difference of successive measurements divided by $\sqrt{\mathrm{2}}$. As just mentioned, the first method (east component: 0.55 m s⁻¹, north component: 0.66 m s⁻¹) overestimates the true noise. However, the second estimate (both components: 0.33 m s⁻¹) is probably too optimistic because it assumes that slow changes of the wind only come from the actual variability of the wind but not from the instrument or the calibration method.

Comparing Fig. 4f and g shows that the calibration of the anemometer narrows the histogram of the derived wind relative to ground. The average calibrated total wind speed amounts to 1.5 m s⁻¹. The scatter of the total wind speed has been computed in the same manner as done for the components and amounts to 0.29 and 0.60 m s⁻¹ for the optimistic and pessimistic computation, respectively. Using the uncalibrated measurements would result in a higher average total wind speed of 1.77 m s⁻¹ and a reduced precision between 0.40 and 0.86 m s⁻¹.

For stationary conditions, pitch and roll could serve as proxies for the wind speed in the x and y directions because the UAV has to tilt in order to resist the wind. We use the same calibration method to find calibration coefficients for pitch and roll as explanatory variables instead of the measured uncalibrated anemometer signal. Overall, this relatively simple method to get wind information from the attitude of the UAV instead of an anemometer gives also reasonable results but with the drawbacks of a lower precision (between 0.74 and 0.89 m s⁻¹) and that a stricter filtering for non-static conditions has to be applied because of pitch and roll due to acceleration. Specifically, this was the driver to use only measurements with a horizontal total acceleration of less than 0.5 m s⁻². Choosing considerably larger threshold values results in much larger spikes in the derived attitude-based wind components.

Note that the anemometer-based wind measurements are also somewhat influenced by accelerations of the UAV but to a far lesser extent. This comes from pitch-and-roll maneuvers of the UAV, resulting in movements of the anemometer at the top of the pole and from variations of the rotor downwash.

5 Validation using ICOS measurements

In order to assure that all instruments of the sUAS work as expected in flight, we performed two validation flights on 9 April 2020 at the ICOS atmospheric station Steinkimmen (STE) near Bremen, Germany (Fig. 5), which is hosted on a 285 m broadcasting tower of the NDR (Norddeutscher Rundfunk). During the flights, the average wind speed was 2.79 m s⁻¹ mainly in the east direction (Fig. 5c).

Figure 5 Overview of the validation flights on 9 April 2020 at the ICOS atmospheric station Steinkimmen (STE) near Bremen hosted on a broadcasting tower of the NDR. (a) Positions of measurements of both flights projected onto an aerial photograph. (b) Top view of the measurement site including launch point (H), position of measurements (stars and circles), and position of the lattice tower (M). (c) Histogram of wind speed and direction. (d) ICOS instruments and orientation of the lattice tower; image courtesy of the German meteorological service (DWD).

At this ICOS station, greenhouse gas concentrations (CO₂, CH₄, and N₂O) and meteorological parameters (e.g., wind speed and direction, temperature, and humidity) are continuously measured at five different levels with altitudes of 32, 82, 127, 187, and 252 m (https://icos-atc.lsce.ipsl.fr/panelboard/ste, last access: 5 January 2021). We use ICOS level 1 data with a temporal resolution of 1 min (Kubistin et al., 2020).

The CO₂ concentration is measured with a Picarro G2301 cavity ring-down spectroscopy analyzer, and each measurement level is sampled for 5 min so that a full measurement cycle takes 25 min. The 1 min averaged CO₂ concentrations are precise to about ±0.1 ppm (1 σ). The wind in each height is measured with Thies Clima ultrasonic 2-D anemometers (accuracy: 0.1 m s⁻¹ for wind speeds up to 5 m s⁻¹) and is also averaged to 1 min intervals. Temperature and humidity at each height are measured with Vaisala HMP155 sensors.

We performed two validation flights with a distance of about 100 m to the tower. In order not to risk signal interference (99 % or 700 kW of the broadcasting power is emitted above 200 m), we only visited the four lower measurement levels.

Figure 6 Measurements and results from the validation flights at the ICOS atmospheric station Steinkimmen (STE), Germany, 9 April 2020. From top to bottom: CO₂ mole fraction (a), wind speed in the east direction (b), wind speed in the north direction (c), temperature (d), relative humidity (e). Left: flight 1. Right: flight 2. Gray areas correspond to measurement sequences S1–S12 during which the UAV visited one of the ICOS measurement levels. Light colors represent filtered measurement with a UAV acceleration larger than 0.1 m s⁻². Bold colors represent measurements used for the validation.

During the first flight, the UAV ascended twice to 187 m and descended level by level to 32 m, hovering about 45 s at each measurement level. During the second flight, we visited each measurement level only once, but instead of hovering, we flew three times back and forth on a 40 m track so that the UAV was always moving while measuring. In total, there are 12 measurement sequences (Figs. A1 and 6, S1–S12) where the UAV visited one of the four measurement levels. The flight pattern is shown in detail in Fig. A1.

For each measurement level, we linearly interpolate the ICOS CO₂, wind, temperature, and humidity to the times of our measurements. We only use those measurements for comparisons which fall into one of the 12 measurement sequences and for which the UAV acceleration is below 0.1 m s⁻². Here, we use a rather strict threshold for the maximum accepted acceleration because we wanted the results to be dominated by steady-state conditions, as expected when flying at a constant speed.

Our CO₂ measurements (Fig. 6a) have basically no systematic offset relative to the ICOS observations (0.07 ppm) and the standard deviation of the difference amounts 2.12 ppm. The ICOS CO₂ measurements vary negligibly during both flights and from height to height (Fig. 6a), which indicates well-mixed atmospheric conditions so that no significant representation errors, i.e., differences due to atmospheric variability, have to be expected. The measured CO₂ scatters slightly more in the second flight compared to the first (especially when the filtered measurements are included). This could be due to fluctuating pressure conditions caused by the rotors during the accelerations when flying back and forth in sequence S9–S12. As discussed in Sect. 3, we multiply the derived scatter with $\sqrt{\mathrm{2}}$ and get 3.00 ppm as the estimate for the in-flight precision of our CO₂ measurements. This is similar to the noise error estimated in the laboratory (Sect. 3) and agrees with the precision as specified by the manufacturer.

We find no significant drifts from height to height which would hint, e.g., at strong uncorrected pressure or temperature dependencies and we also find no significant systematic offsets between both flights. According to Sect. 3, we estimate that CO₂ enhancements can be measured in flight with an accuracy (or trueness) of 1.3 % or 0.9 ppm, whichever is larger.

The average of the north component of the wind speed (Fig. 6c) is only slightly larger (0.1 m s⁻¹) compared to the ICOS measurements, and the standard deviation of the difference amounts to 0.90 m s⁻¹. The relative large scatter of the difference to ICOS is mainly driven by a poor agreement in S12. As visible in Fig. 5a and b, the tower is located close to a small piece of forest, while the flight track is above a relatively free field which can lead to significant differences of the measurements at 32 m height. Additionally, it shall be noted that the second flight has a larger distance to the tower compared to the first flight. When considering only S1–S11, the average difference becomes 0.02 m s⁻¹ and the standard deviation of the difference is reduced to 0.67 m s⁻¹.

The east component of the wind (Fig. 6b) is on average 0.9 m s⁻¹ larger compared to the ICOS measurements and the standard deviation of the difference is 0.57 m s⁻¹. The ICOS anemometers are mounted on short outriggers at the west side of the tower, leading to a possible underestimation of the ICOS measurements especially in the east direction due blockage of eastbound winds as prevailing during the validation flights (Fig. 5b and c). This would be consistent with the generally larger east component of the wind speed measured with the sUAS compared to the ICOS measurements, independent of the orientation of the sUAS so that a systematic overestimation in one direction in UAV coordinates can be excluded as explanation.

The average of the scatter of the north and east component is 0.62 m s⁻¹ (excluding S12 of the north component). In addition to the noise of our measurements, this includes the representation error and the noise of the ICOS measurements. Therefore, it can be considered an upper limit of the scatter of our measurements of a horizontal wind component. This agrees well with the upper limit of the scatter determined in Sect. 4 and compares to the lower limit of 0.33 m s⁻¹ computed in the same section. For convenience, we estimate that the precision of the individual measurements of a horizontal wind component lies in the middle of both values (0.48 m s⁻¹).

In order to estimate the systematic uncertainties of our wind component measurements, we compute the average differences in sequences S1–S11 of the north component. The standard deviation of these values amounts to 0.34 m s⁻¹ and is considered an estimate of the accuracy of our wind component measurements.

Due to the overestimation of the east component (or underestimation by ICOS), also the total horizontal wind speed is larger than the ICOS observations. On average, the difference amounts to 0.84 m s⁻¹ and has a standard deviation of 0.51 m s⁻¹ (when excluding S12). We estimate from the latter value and from the lower limit of the scatter (Sect. 4) that the precision of the individual measurements of the total horizontal wind speed is 0.40 m s⁻¹, which is consistent with the manufacturer's specification of ±0.3 m s⁻¹ for wind speeds below 16 m s⁻¹. The accuracy of the total horizontal wind speed is 0.33 m s⁻¹ and has been estimated in the same way as for the wind components.

The temperature measured with the internal temperature sensor of the GMP343 CO₂ sonde is typically about 1 K larger compared to the ICOS measurements (Fig. 6d). However, it can reasonably well reproduce the temperature profile with smaller values at higher altitudes, even though, especially within the first 3–4 min after launch (S1, S2, and S9), one can clearly see the hysteresis of the sensor resulting in an overestimation of up to 2 K. A potential explanation of the general overestimation could be the heating of the GMP343 CO₂ sonde, which is intended to reduce the possibility of condensation on the optical components but which could also slightly warm the temperature sensor. Relative humidity (Fig. 6e) is typically up to 10 % lower than measured by ICOS, and the sensor appears to have a relatively long response time on the order of some minutes.

6 Elevated CO₂ concentrations downwind of an industrial facility

In order to demonstrate how the sUAS can be used to measure elevated CO₂ concentrations of a nearby source, we performed two flights on 10 April 2020 near the ExxonMobil natural gas processing facility in Großenkneten about 40 km west of Bremen, Germany (Fig. 7). In this facility, hydrogen sulfide (H₂S) is removed from natural gas (i.e., sour gas is “sweetened”), which is an energy-intense process. According to the European Pollutant Release and Transfer Register (E-PRTR), the annual CO₂ emissions of the facility declined from 1260 kt in 2014 to 1010 kt in 2017 which is the most recent value in the E-PRTR database (https://prtr.eea.europa.eu, last access: 11 June 2020). Similar values for the CO₂ equivalent emissions can be found in the most recent (2019) list of stationary facilities in Germany that are subject to emissions trading published by the German emissions trading authority (DEHSt; https://www.dehst.de, last access: 11 June 2020) of the German Environment Agency, but in addition to the E-PRTR database, this list also includes values for 2018 (855 kt) and 2019 (892 kt).

Figure 7 Overview of two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten 40 km west of Bremen, Germany, 10 April 2020. (a) CO₂ measurements of both flights projected onto an aerial photograph taken from a height of 80 m. Note that only for the purpose of plotting, we added 3 m to the altitudes of the first flight. (b) Top view of the measurement site including positions of the sources S1–S3, flight track, and histogram of wind speed and direction (c).

It can be expected that most of the CO₂ emissions are released by the two 150 m high stacks, prominently visible in Fig. 7 (S1). In 2014, a combined heat and power plant was put into operation and from information on the ExxonMobil website (https://www.erdgas-aus-deutschland.de, last access: 11 May 2020), one can roughly estimate that about 17 % of the total emissions may be emitted via its 34 m high stack (Fig. 7, S2). Moreover, an unknown but likely small proportion is emitted by gas flaring (Fig. 7, S3). We roughly estimate that 740 kt CO₂ have been released via the main stacks in 2019. For later emissions or the interannual variability, we do not have an estimate.

Targeting the emissions of the main stacks, we flew two vertical cross-sections roughly perpendicular to the average wind direction with a length of 272 m and altitudes ranging from 120 to 200 m. The second cross-section has been shifted by roughly 70 m westwards relative to the first one in order to get a denser sampling in the region with the largest expected enhancements. The number of tracks and velocity of the UAV relative to ground have been chosen to realize a similar sampling density in horizontal (9.1 m) and vertical (8.9 m) direction for 2 s measurement intervals. Due to flight regulations, we limited the maximum altitude to 200 m, although we anticipated that parts of the plume would probably have risen above that. The details of the flight pattern are shown in Fig. A2.

Figure 8 Measurements and results from two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten 40 km west of Bremen, Germany, 10 April 2020. (a) CO₂ mole fraction as well as background values (gray line) estimated from measurements within light gray areas and CO₂ enhancement (light blue area). (b) Wind speed normal including 100 s running average (gray line), filtered measurements with UAV accelerations larger than 0.1 m s⁻² (light blue), and valid measurements (bold blue) used to compute the running average. (c) Air density including a 10 s running average and filtering information (light and bold blue) as for the wind speed normal. Left: flight 1. Right: flight 2.

The CO₂ mass flux density perpendicular through the cross-section resulting from the facility's emissions can be computed by

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }F=w_{\perp }{\mathit{\rho }}_{\text{air}}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}{M}_{{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

Here, w_⊥ is the wind speed normal (perpendicular to the cross-section), ρ_air is the molar density of air, ΔCO₂ is the CO₂ enhancement caused by the facility's emissions, and $M_{{\mathrm{CO}}_{\mathrm{2}}}$ is the molar mass of CO₂. Strictly speaking, Eq. (10) represents only the mass flux density due to advection, and we assume that diffusion can be neglected, which is the case for wind speed normals larger than about 2 m s⁻¹ (Sharan et al., 1996). The sUAS measures all quantities needed to compute the CO₂ flux density for each measurement interval, and Fig. 8 shows the quantities needed to compute Eq. (10).

In order to derive ΔCO₂ from the measured CO₂ mole fraction, we manually define measurement intervals which we consider to represent background concentrations not disturbed by the emissions of the facility. For each of these intervals (gray background in Fig. 8a), we compute the average CO₂ mole fraction and estimate the background concentration for each sounding by linear interpolation. Under the assumption of no additional nearby upwind sources, we compute the CO₂ enhancement caused by the facility's emissions by the difference between the actually measured and the estimated background concentration (light blue area in Fig. 8a). The CO₂ enhancement is often on the order of 100 ppm and reaches maximum values of almost 400 ppm.

For each sounding, we compute the wind speed normal from the calibrated wind measurements and filter out wind measurements with UAV accelerations larger than 0.1 m s⁻². The resulting gaps are filled with linearly interpolated values, and in order to reduce the noise of the wind measurements, we compute a 100 s running average (Fig. 8b). The wind speed reduces with height superimposed by a refreshing of the wind towards the end of the second flight, which we also noticed on the ground. The average wind speed normal amounts to 4.0 m s⁻¹.

By applying the ideal gas law, we compute the molar air density, which is shown in Fig. 8c. We use the same filtering and interpolation as for the wind measurements, but because of less relative noise, we compute a running average with a smoothing kernel of only 10 s instead of 100 s. The most noticeable features visible in the air density plot (Fig. 8c) are the increasing values for lower heights due to the larger pressure.

Figure 9 CO₂ mass flux density perpendicular to the plume cross-section computed from the measurements of two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten 40 km west of Bremen, Germany, 10 April 2020. The x axis corresponds to the distance from the west-most point of the flight pattern.

We compute the flux density according to Eq. (10) and grid it with a spacing of 10 m horizontally and 8.89 m vertically, ensuring that each grid box includes at least one data point (Fig. 9).

Integrating the flux density field results in a total flux of 196 kt CO₂ yr⁻¹. In order to estimate the uncertainty of this value due to instrumental noise, we compute the total flux from measurements to which we add Gaussian noise with a standard deviation according to the instrument characteristics and the results of the ICOS validation (Sect. 5: 3 ppm for CO₂, 0.48 m s⁻¹ for the wind speed normal component, 1 ^∘C for temperature, and 1 hPa for pressure). We repeat this 100 times and get a standard deviation of 2.5 kt CO₂ yr⁻¹, which corresponds to the stochastic 1σ uncertainty.

We also estimate the influence of potential systematic instrumental errors summarized in Sect. 5. The systematic CO₂ uncertainty (1.3 %) corresponds to a systematic total flux uncertainty of 2.5 kt CO₂ yr⁻¹. The systematic wind component uncertainty (0.34 m s⁻¹ = 8.3 %) translates to a flux uncertainty of 16.3 kt CO₂ yr⁻¹. Assuming 1 ^∘C to be the systematic uncertainties of the temperature measurement results in total flux uncertainties of 0.7 kt CO₂ yr⁻¹ (3.5 ‰). Systematic uncertainties introduced by the pressure measurements can safely be neglected.

Adding up the variances of the systematic uncertainties plus the stochastic uncertainty results in a total uncertainty of 16.7 kt CO₂ yr⁻¹ or 8.5 %. It shall be noted that this uncertainty estimate only includes systematic and random uncertainties introduced by the instruments. There are other effects which also can result in significant differences between the inferred total flux and the facility's emissions such as turbulence (i.e., short-term fluctuations due to non-steady meteorology) or undersampling of the plume morphology. Most of such effects will average out with increasing number of flights or increasing flight time. How well this works and how large the corresponding uncertainties actually are can be assessed, e.g., by analyzing the variability of results from multiple successive flights or by high-resolution plume simulations. However, such assessments are beyond the scope of this publication.

The inferred total flux of 196±17 kt CO₂ yr⁻¹ is significantly lower than our rough estimate of the annual emissions of the main stacks for the year 2019 based on DEHSt (740 kt CO₂ yr⁻¹). The uppermost tracks of both flights contain some significantly elevated CO₂ measurements (Figs. 7, 8, 9), strongly indicating that we have only seen parts of the plume, and we assume that most of the plume has risen above 200 m. In addition, we see some non-background values in the lowermost tracks which could either come from the nearby gas flare but also indicate that parts of the plume may have been below 120 m.

We use a Gaussian plume model (Beychok, 2005) to estimate the expected plume extent for moderately unstable conditions (Pasquill stability class B), resulting in a full width half maximum of 197 m horizontally and 124 m vertically. The expected corresponding plume rise can be estimated with Briggs' equations for bent-over, hot buoyant plumes (Beychok, 2005). Most input parameters to the Briggs equations such as temperature and wind speed at stack height have been measured, but other parameters require ad hoc assumptions. We assume that the exiting flue gas consists of 21 % CO₂ and that it is 50 ^∘C warmer than the ambient air. By applying the ideal gas law, these values are used to estimate that the annual emissions through the main stacks have an average volumetric flow rate of roughly 78 m³ s⁻¹. For this scenario, Briggs' equations estimate that the expected center of the plume 500 m downwind of the source has risen to 234 m.

In the case of a Gaussian plume morphology, this would mean that roughly 74 % of the emitted CO₂ has risen above 200 m. However, it shall be noted that the Gaussian plume shape and a plume rise according to Briggs' equations is only on average a good estimate for reality, but on short timescales, turbulence can result in large deviations form that. Also, the results of our simple simulations are relatively sensitive to the ad hoc assumptions. Nevertheless, they indicate that the width of the flight pattern (about 340 m) was sufficient to sample the expected plume width but that large parts of the plume may indeed have risen above the maximum flight altitude.

The fact that we use kt CO₂ yr⁻¹ as the unit for the total flux does not imply that we actually determine annual emissions; our measurements still only correspond to a snapshot of the current situation and the unit has only been chosen because it is commonly used. The estimation of annual emissions usually requires repeated flights throughout the year, and assumptions on the variability of the actual emissions have to be made to estimate the corresponding uncertainties (Velazco et al., 2011).

7 Summary and conclusion

We introduced a small sUAS measuring all atmospheric quantities needed to quantify the CO₂ emissions of a nearby point source from its downwind mass flux. The entire sUAS weighs about 6 kg and has been built from commercially available components, which allowed us to realize an affordable but reliable system in a relatively short development phase. In order to quantify the atmospheric CO₂ mass flux without the need for any ancillary data, the payload includes a CO₂ probe, an anemometer, and sensors for temperature, pressure, and relative humidity. The onboard computer uses a separate radio data link, which allows the flight pattern to be adapted to the CO₂ and wind conditions in flight.

With special emphasis on CO₂ and wind, we performed calibration, comparison, and validation experiments in the laboratory and in flight. We introduced a method to calibrate the anemometer under flight conditions when the headwind has to be accounted for and rotor downwash can influence the local wind field. We validated our CO₂ and wind measurements by comparison with ICOS measurements at the 285 m high NDR broadcasting tower in Steinkimmen near Bremen, Germany.

According to the validation experiment, an upper limit for the 1σ single sounding precision of our CO₂ measurements during a flight is 3 ppm (measurement interval of 2 s), and from a linear analysis of correlated errors in the laboratory, we conclude that CO₂ enhancements can be measured with an accuracy of 1.3 % or 0.9 ppm, whichever is larger.

Our anemometer calibration method derives the free fit parameters of a linear calibration model accounting for scaling, rotation, and a potential constant bias. For this purpose, it analyzes wind measurements taken while following a suitable flight pattern and assuming stationary wind conditions. From the calibration and validation experiments, we estimated that we can retrieve the total horizontal wind speed relative to ground during a flight with a 1σ single measurement precision of 0.40 m s⁻¹ and an accuracy of 0.33 m s⁻¹. The precision is similar to the root mean square errors of 0.5 and 0.4 m s⁻¹ found by Palomaki et al. (2017) and Shimura et al. (2018), respectively, but in both cases a bias of 0.5 m s⁻¹ has previously been subtracted.

For the purpose of flux estimation, the uncertainty of the wind speed components is more relevant than the uncertainty of the total wind speed. We estimated that these can be retrieved with a 1σ single measurement precision of 0.48 m s⁻¹ and an accuracy of 0.34 m s⁻¹. This compares to an accuracy of 0.5 m s⁻¹ of the wind components measured aboard a manned aircraft by Krings et al. (2018).

Additionally, we showed that our calibration method can also be used to derive wind information from the attitude of the UAV. The inferred wind speeds are reasonably consistent with those of the calibrated anemometer but feature a reduced quality. It cannot be excluded that the attitude-based wind retrieval might be significantly improved in the future, but for the time being, our findings are consistent with those of Palomaki et al. (2017) and Barbieri et al. (2019) who concluded that sonic anemometers provide the most accurate wind information from multi-rotor platforms.

During two flights, we measured the CO₂ enhancement downwind of the ExxonMobil natural gas processing facility in Großenkneten about 40 km west of Bremen, Germany, and demonstrated how the measurements of the sUAS can be used to infer mass fluxes of atmospheric CO₂ related to the emissions of the facility. We aimed at the emission of the two 150 m high main stacks of the facility and centered the cross-sectional flight patterns according to the wind direction and found enhancements of up to almost 400 ppm. Integration over the inferred flux densities resulted in a total cross-sectional flux of 196±17 kt CO₂ yr⁻¹. This is significantly lower than our rough estimate of 740 kt CO₂ yr⁻¹ for the annual emissions of the main stacks for the year 2019 based on DEHSt.

Quantitative comparisons usually require that complete cross-sections of the plume can be examined. Whilst we measured mostly background concentrations at the horizontal edges of our flight pattern, we measured some significant CO₂ enhancements at the vertical edges, especially at the top in a height of 200 m, which was the maximum altitude allowed by flight regulations for that day. This strongly indicates that we have only seen parts of the plume, and we assume that most of the plume has risen above the maximum flight altitude. Additionally, our result represents only a snapshot of the current situation which does not necessarily agree with the annual average emissions. Ideally, emissions should be compared only on the basis of instantaneous or short-term average values.

Our uncertainty estimate is not affected by the fact that our measurements sampled only parts of the plume. It suggests that it is possible to determine the cross-sectional flux of a point source with this type of sUAS with an uncertainty of about 8.5 % when considering only instrumental effects and neglecting, e.g., the influence of turbulence. The flux uncertainty is dominated by the uncertainty of the wind speed and valid for average wind speeds of about 4 m s⁻¹. It is relatively insensitive to the magnitude of the flux and will reduce for moderately larger wind speeds, even though the CO₂ enhancements will become smaller. Our uncertainty estimate is consistent with that of Krings et al. (2018) for in situ measurements taken aboard a manned aircraft with a larger distance to the source. They estimated the flux error due to the uncertainties of the primary measurements to be well below 10 %.

It shall be noted that our uncertainty estimate only accounts for instrumental effects, and we agree with Krings et al. (2018) that there are other factors which can result in significant differences between the inferred total flux and the facility's actual emissions. These are, for example, turbulence (i.e., short-term fluctuations due to non-steady meteorology) or undersampling of the plume morphology. However, such effects can average out when averaging over a large enough number of cross-sections, i.e., long-enough flight times so that the uncertainty converges towards our estimate of systematic instrumental effects.

The introduced sUAS can also be used for weaker or stronger sources as long as a suitable distance to the source and extent of the flight pattern can be used. However, flight regulations often prescribe a minimum distance to be kept and, additionally, flights are usually only permitted within visual line of sight, which limits the maximum extension of the flight pattern. One main target of future satellite missions such as CO2M is the quantification of anthropogenic CO₂ point source emissions and with the estimated accuracy, the sUAS is suitable to complement other, e.g., airplane-based validation activities. Moreover, it may also be useful for the sample verification of reported emissions and can provide emission estimates under conditions not suitable for satellite measurements, e.g., in cloud-contaminated scenes or during night (if flight regulations allow).

However, it shall be noted that there are also several weather conditions under which flights with our sUAS are not possible or meaningful. For example, strong winds may prevent a save operation, rain may destroy parts of the payload or prevent its functioning, and too-low wind speeds may render the results hard to interpret because the mass flux becomes dominated by diffusion instead of advection. In addition to this, many countries have put in place strict legal regulations for UAV flights (for good reasons), which also impact the possible areas of application. These regulations vary significantly from country to country and are subject to frequent adjustments, and they often also depend on which organization is responsible for the flights. For example, some government agencies can have far-reaching permissions. In addition, the local aviation authority may grant special permissions, as was the case with our flights.

For the future, we envisage many interesting potential applications, improvements to be made, and scientific questions to be answered. A promising next step would be to use the sUAS for the quantification of the emissions of known sources by measuring complete plume cross-sections and investigating how much averaging has to be applied on the cross-sectional flux in order to converge to the actual emissions given short-term fluctuations due to turbulence and potential undersampling of the complex plume morphology. As the wind measurement is the dominating source of uncertainty, it should be investigated how these measurements can be further improved. This could, e.g., include experiments with an anemometer with reduced noise or the implementation of more complex calibration models. In addition, the capabilities of the onboard computer should be exploited to allow fully autonomous flights with a flight pattern automatically adapted to the wind and CO₂ measurements, resume after battery replacement, complex no-fly zones, etc. The change to a more advanced UAV could enable longer flight times, heavier payloads (e.g., for simultaneous measurements of other gases), and flights at higher wind speeds.

Appendix A

Figure A1 Detailed flight pattern of the validation flights at the ICOS atmospheric station Steinkimmen (STE), Germany, 9 April 2020. From top to bottom: altitude (a), horizontal distance from the first valid measurement within each flight (b), and yaw of the UAV (c). Left: flight 1. Right: flight 2.

Figure A2 Detailed flight pattern of two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten 40 km west of Bremen, Germany, 10 April 2020. (a) Altitude. (b) Horizontal distance from the west-most point of the flight pattern. Left: flight 1. Right: flight 2.

Figure A3 Temperature (a), relative humidity (b), and pressure measured (c) during two flights downwind of the ExxonMobil natural gas processing facility in Großenkneten 40 km west of Bremen, Germany, 10 April 2020. Left: flight 1. Right: flight 2.