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

Abstract. A reduction of the anthropogenic emissions of CO2 (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 CO2 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 CO2 emissions of a nearby point source from its downwind mass flux without the need for any ancillary data.
Specifically, CO2 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 CO2 measurements have a precision of 3 ppm and that CO2 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−1 and the accuracy to be 0.33 m s−1.
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 CO2 concentrations can be used to infer mass fluxes of atmospheric CO2 related to the emissions of the facility.


15kC including the UAV, which has the highest unit price followed by the CO 2 sensor.

CO 2 sensor characterization
Flux quantification of an individual source requires measurements of the difference between elevated concentrations minus 95 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 100 to temperature, pressure, and humidity and proposed to derive an empirical correction which can be specific for a sensor unit.
We compare roughly one week of continuous measurements of our CO 2 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 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. Average power consumption during flight including camera and for the total take-off weight of 6022g. Considering only the weight of the camera. Including meteorological sensors, aluminum housing, and power converter.

Device
Weight 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 2 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 2 , all parameters vary more than we expect for a typical measurement flight.
Throughout the entire comparison time, the CO 2 measurements are always close together with an average difference of -110 0.33 ppm. We compute the 1h running average of the difference between the GMP343 and the LGR instrument. Deviations from zero can safely be assumed to be not 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 2 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 115 test-flights with very variable CO 2 concentrations near a stack revealed that always two GMP343 measurements (sampling rate 0.5 Hz) are paired, i.e., they exhibit similar CO 2 values within some ppm and Vaisala's technical support confirmed that truly independent measurements are only possible every 4 seconds. In order to account for this, we multiply the derived standard deviation by √ 2 and get 2.49 ppm as estimate for the noise of the individual soundings. This is an upper limit of the measurement noise because the multiplication by √ 2 would imply that both measurements within a 4 seconds interval are 120 identical, which is not the case.
In order to analyze the drifts of the 1h 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 −3 ), and the measured CO 2 concentration (GMP343 in ppm). Similar  Note, the absolute humidity is computed from the relative humidity (h r in %) by: wherein SVD is the saturated vapor density (in g m −3 ), which can be well approximated for temperatures between 0°C and 130 40°C by the following equation (Nave, 2017): 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 135 correction comes from the apparent underestimation of 1.3% of the CO 2 variability as measured by the GMP343.
However, it can be expected that other factors (especially, wind speed, see Sec. 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.

140
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 an 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 145 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 x-and y-direction 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.
Applying the rotation matrix R with γ being the azimuthal orientation (yaw) of the UAV results in the calibrated wind at the anemometer but in geographic coordinates: The wind at the anemometer is a superposition of the wind relative to ground w and the headwind because of the velocity v 155 of the UAV (u = w − v), so that the wind relative to ground becomes A straight forward 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 160 GPS is not available.
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).

165
Specifically, we solve Eq. 7 for the uncalibrated anemometer measurement 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 −2 . The results for the free fit parameters show that the calibration scales the measured wind speed in x-direction with 0.881 and 0.832 in y-direction. A rotation is basically not applied and there is only a very small constant offset correction below 0.1 m s −1 needed.
As shown in Fig. 4 (a, b), the fitted wind measurements (right side of Eq. 8, green) agree well with the actually measured values (left side of Eq. 8, blue). The derived wind relative to ground w ( Fig. 4 c, d, green) scatters around the mean values of 175 1.06 m s −1 in east and west direction. 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 −1 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 180 estimate the noise of the measurements by the standard deviation of the difference of successive measurements divided by √ 2.
As just mentioned, the first method (east component: 0.55 m s −1 , north component: 0.66 m s −1 ) overestimates the true noise.  and pitch as well as corresponding fits. c) Wind speed relative to ground in east direction for the calibrated (green) and uncalibrated (blue) anemometer as well as computed from the attitude of the UAV (orange). d) Wind speed relative to ground in north direction. Histograms of wind speed and direction for the calibrated (e) and uncalibrated (f) anemometer as well as computed from attitude (g).
However, the second estimate (both components: 0.33 m s −1 ) 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. 4 f) and g) shows that the calibration of the anemometer narrows the histogram of the derived wind relative For stationary conditions, pitch and roll could serve as proxies for the wind speed in x-and y-direction because the UAV 190 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 m s −1 and 0.89 m s −1 ) 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 195 acceleration of less than 0.5 m s −2 . 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.  We performed two validation flights with a distance of about 100 m to the tower. In order not to risk signal interference (99% 215 or 700 kW of the broadcasting power is emitted above 200 m), we only visited the four lower measurement levels.
During the first flight, the UAV twice ascended 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 long track so that the UAV was always moving while measuring. In total there are 12 measurement sequences ( Fig. A1 and 6, S1 -S12) where the UAV visited one of the four measurement levels. The flight pattern 220 is shown in detail in Fig. A1.
For each measurement level, we linearly interpolate the ICOS CO 2 , 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 −2 . We here 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.

225
Our CO 2 measurements (Fig. 6, a) 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 2 measurements vary negligibly during both flights and from height to height (Fig. 6, a) 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 2 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 230 conditions caused by the rotors during the accelerations when flying back and forth in sequence S9 -S12. As discussed in Sec. 3, we multiply the derived scatter with √ 2 and get 3.00 ppm as estimate for the in-flight precision of our CO 2 measurements. This is similar to the noise error estimated in the laboratory (Sec. 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 235 dependencies and we also find no significant systematic offsets between both flights. According to Sec. 3, we estimate that CO 2 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. 6, c) is only slightly larger (0.1 m s −1 ) compared to the ICOS measurements and the standard deviation of the difference amounts to 0.90 m s −1 . The relative large scatter of the difference to ICOS is mainly driven by a poor agreement in S12. As visible in Fig. 5 (a and b), the tower is located close to a small piece of 240 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 −1 and the standard deviation of the difference reduces to 0.67 m s −1 .
The east component of the wind (Fig. 6, b) is on average 0.9 m s −1 larger compared to the ICOS measurements and the In order to estimate the systematic uncertainties of our wind component measurements, we compute the average differences in sequence S1-S11 of the north component. The standard deviation of these values amounts to 0.34 m s −1 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 260 than the ICOS observations. On average, the difference amounts to 0.84 m s −1 and has a standard deviation of 0.51 m s −1 (when excluding S12). We estimate from the latter value and from the lower limit of the scatter (Sec. 4) that the precision of the individual measurements of the total horizontal wind speed is 0.40 m s −1 , which is consistent with the manufacturer's specification of ±0.3 m s −1 for wind speeds below 16 m s −1 . The accuracy of the total horizontal wind speed is 0.33 m s −1 and has been estimated as for the wind components.

265
The temperature measured with the internal temperature sensor of the GMP343 CO 2 sonde is typically about 1 K larger compared to the ICOS measurements (Fig. 6, d). However, it can reasonably well reproduce the temperature profile with smaller values at higher altitudes, even though, especially within the first 3-4 minutes 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 2 sonde, which is intended to reduce the possibility of condensation 270 on the optical components but which could also slightly warm the temperature sensor. Relative humidity (Fig. 6, e) is typically up to 10% lower than measured by ICOS and the sensor appears to have a relatively long response time in the order of some minutes.

Elevated CO 2 concentrations downwind of an industrial facility
In order to demonstrate how the sUAS can be used to measure elevated CO 2 concentrations of a nearby source, we performed 275 two flights on 10.04.2020 near the ExxonMobil natural gas processing facility in Großenkneten about 40 km east of Bremen, Germany (Fig. 7). In this facility, hydrogen sulfide (H 2 S) is removed from natural gas (i.e., sour gas is "sweetened"), which is an energy intense process. It can be expected that most of the CO 2 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 285 (https://www.erdgas-aus-deutschland.de, last access 11.05.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 2 have been released via the main stacks in 2019. For later emissions or the inner-annual variability, we do not have an estimate.
Targeting at the emissions of the main stacks, we flew two vertical cross-sections roughly perpendicular to the average  The CO 2 mass flux density perpendicular through the cross-section resulting from the facility's emissions can be computed by: Here w ⊥ is the wind speed normal (perpendicular to the cross-section), ρ air is the molar density of air, ∆CO 2 is the CO 2 300 enhancement caused by the facility's emissions, and M CO2 is the molar mass of CO 2 . 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 −1 (Sharan et al., 1996). The sUAS measures all quantities needed to compute the CO 2 flux density for each measurement interval and Fig. 8 shows the quantities needed to compute Eq. 10.
In order to derive ∆CO 2 from the measured CO 2 mole fraction, we manually define measurement intervals which we 305 consider to represent background concentrations not disturbed by the emissions of the facility. For each of these intervals (gray background in Fig. 8, a), we compute the average CO 2 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 2 enhancement caused by the facility's emissions by the difference between the actually measured and the estimated background concentration (light blue area in Fig. 8, a). The CO 2 enhancement is often in the order of 100 ppm and reaches maximum values 310 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 −2 . 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. 8, b). 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 ground. The 315 average wind speed normal amounts to 4.0 m s −1 .
By applying the ideal gas law, we compute the molar air density, which is shown in Fig. 8 (c). 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 feature visible in the air density plot (Fig. 8, c) are the increasing values for lower heights due to the larger pressure. 320 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 −1 and an accuracy of 0.34 m s −1 . This compares to an accuracy of 0.5 m s −1 of the wind components measured aboard a manned aircraft by Krings et al. (2018).

375
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 base 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. Of course, quantitative comparisons require that complete cross-sections of the plume could be examined. Whilst we measured mostly background concentrations at the horizontal edges of our flight pattern, we measured some significant CO 2 enhancements at the vertical edges, especially at the top in a height of 200 m, which was the maximum altitude allowed by 390 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 395 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 −1 . It is relatively insensitive to the magnitude of the flux and will reduce for moderately larger wind speeds, even though the CO 2 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 in 400 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 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 However, it shall be noted that there are also several weather conditions under which flights with our sUAS are not possible 415 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. Beside this, many countries have put in place strict regulations for UAV flights (for good reasons), which also impacts the possible areas of application.
For the future, we envisage many interesting potential applications, improvements to be made, and scientific questions to 420 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 investigate 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 425 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 2 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.