Estimating Upper Silesian coal mine methane emissions from airborne in situ observations and dispersion modeling

. Abundant mining and industrial activities located in the Upper Silesian Coal Basin (USCB) lead to large emissions of the potent greenhouse gas (GHG) methane (CH 4 ). The strong localization of CH 4 emitters (mostly conﬁned to known coal mine ventilation shafts) and the large emissions of 448 / 720 kt CH 4 yr − 1 reported in the European Pollutant Release and Transfer Register (E-PRTR 2017) and the Emissions Database for Global Atmospheric Research (EDGAR v4.3.2) make the USCB a prime research target for validating and improving CH 4 ﬂux estimation techniques. High-precision observations of this GHG 5 were made downwind of local (e.g. single facilities) to regional-scale sources (e.g. agglomerations) in the context of the CoMet 1.0 campaign in early summer 2018. A Quantum Cascade/Interband Cascade Laser (QCL/ICL) based spectrometer adapted for airborne research was deployed aboard the German Aerospace Centers (DLR) Cessna 208B to sample the planetary boundary layer (PBL) in situ. Regional CH 4 emission estimates for the USCB are derived using a model approach including assimilated wind soundings from three ground-based Doppler lidars. Although retrieving estimates for individual emitters is difﬁcult using 10 only single ﬂights, due to sparse data availability, the combination of two ﬂights allows for exploiting different meteorological conditions (analogous to a sparse tomography algorithm) to establish conﬁdence on facility level estimates. Emission rates from individual sources are not only needed for unambiguous comparisons between bottom-up and top-down inventories but become indispensable if (independently veriﬁable) sanctions are to be imposed on individual companies emitting GHGs. An uncertainty analysis is presented for both the regional scale and facility level emission estimates. We ﬁnd instantaneous 15 emission estimates of 452 / 442 ± 78 / 75 kt CH 4 yr − 1 for the morning / afternoon ﬂight of June 6th, 2018. The derived emission rates coincide ( ± 2 % ) with annual-average inventorial data from E-PRTR 2017 albeit they are distinctly lower (-37 % / -40 % ) than values reported in EDGAR v4.3.2. Discrepancies in available emission inventories could potentially be narrowed down with sufﬁcient observations using the method described herein to bridge the gap between instantaneous emission estimates and yearly averaged inventories. ﬂight

emission estimates of 436 ± 115 kt CH 4 yr −1 and 477 ± 101 kt CH 4 yr −1 from two research flights along with a detailed 55 uncertainty analysis. The present study is based on the very same research flights and aims at contributing an advanced model approach. Previous studies have used Lagrangian models to simulate the dispersion (Tuccella et al. (2017); Raut et al. (2017)) of plumes emanating from oil and gas platforms or identification of CH 4 sources (Platt et al., 2018). Atmospheric transport models have been used to infer CH 4 emissions from the oil and natural gas industry (Barkley et al., 2017). Here, a combination of a Eulerian atmospheric transport model and a Lagrangian particle dispersion model is used in conjunction with assimilated 60 Doppler lidar soundings to infer instantaneous CH 4 emissions for Europe's largest coal extraction region, the USCB.
Section 2 provides an introduction on the USCB as the region of prime interest followed by a research flight overview in Sect. 3. Section 4 details a model based flux estimation approach. CH 4 emission estimates will be given in the form of a case study in Sect. 4.1 for two research flights on June 6th, 2018 along with an estimate of the uncertainties involved. Section 5 summarizes our findings and concludes the study.  zinc and lead ore exploitation. Coal mining activities make up for the largest part, with an approximate total of 10 billion metric tonnes extracted since the industrial revolution, where over 70 % of this exploitation took place after 1945. To date an approximate 75 million tonnes of coal are extracted every year from 27 active mines. It is these figures and the large area of approximately 7400 km 2 covered, that make the USCB the largest coal extraction region in Europe (Dulias, 2016). The intensive coal mining activities and the heavy industry spread around the city of Katowice, Poland, located in the north of the 75 USCB, lead to significant amounts of GHG emissions in this area. Fugitive CH 4 emanating from the coal mine shafts reaching several hundred meters into the ground is either actively ventilated (active mines) or degasses passively from abandoned mines.
Mines located in the north of the USCB are mostly abandoned and partially flooded, while intensive, active coal exploitation is located in the southern USCB, both in Poland and the Czech Republic (Gzyl et al., 2017).  based on a subset of the publicly available EDGAR v4.3.2 CH 4 (https://edgar.jrc.ec.europa.eu/) emission inventory (Janssens-Maenhout et al., 2017). It shows CH 4 emissions range up to approx. 100 kt yr −1 on a 0.1 × 0.1 degree grid with source strengths increasing towards the southern USCB. Accordingly, the strongest sources are located near the Czech border mid ways between the cities of Bielsko-Biala, Poland and Ostrava, Czech Republic. According to EDGAR v4.3.2, these CH 4 sources are among the strongest in Europe. The total CH 4 emissions from this inventory amount to approximately 720 kt yr −1 for the USCB 85 region, where ∼ 620 kt yr −1 are attributed to the fuel exploitation sector. Facility level emission data of CH 4 are provided by E-PRTR 2017. The locations of 74 documented coal mine ventilation shafts (active and inactive) have been added to Fig. 1 for reference. These locations were visually identified from satellite imagery and emission values from E-PRTR 2017 were evenly distributed among the ventilation shafts for each company (see also Nickl et al. (2019)

Research flight overview
The CoMet mission in early summer 2018 primarily aimed at providing observations of GHG (mainly CO 2 and CH 4 ) gradients along large-scale latitudinal transects over Europe from co-ordinated operation of several state-of-the-art instruments on ground 95 and aboard five research aircraft. Aboard the Cessna 208B, a rich dataset of simultaneous airborne observations of CH 4 , C 2 H 6 , CO 2 , CO, N 2 O and H 2 O was collected using the QCLS instrument (see Fig. 2 and Kostinek et al. (2019) for details) during ∼30 flight hours. In the following, a subset of these data from two research flights undertaken on June 6th, 2018 was used to retrieve CH 4 fluxes emanating from the USCB region. Both flight tracks are shown in Fig. 1 along with the locations of three co-deployed Leosphere Windcube 200S Doppler wind lidars (Wildmann et al., 2020). The morning (black line) and afternoon 100 flights (red line) circumvent all known ventilation shafts in the area (gray triangles) and are in fact very similar (congruent) from the top-down perspective. This is well intended to enhance confidence on retrieved GHG fluxes. Moderate (3 -6 m s −1 ) winds throughout the day from north-easterly directions drive advection of the CH 4 plumes towards the Czech border and into the Ostrava region.  The model-based approach developed in this work employs a combination of Eulerian and Lagrangian particle dispersion models. Due to the known locations of the coal mine ventilation shafts, their emissions are modeled forward in time with constant emission rates. Modeled data are then extracted at the aircraft position in space and time and compared to actual airborne in situ observations. This comparison depends on the quality of the a-priori emission data, the quality of the measurements and 125 the quality of the transport model simulation, which in turn depends on the quality of the meteorological data (winds, PBL heights, etc.). Validating the meteorological data is therefore important to enable regional emission estimates based on particle dispersion models. Here, meteorological driver data is generated using the Weather Research and Forecasting (WRF) v4.0 model (Powers et al., 2017)    Soundings from three Doppler lidars (marked DLR85, DLR86 and DLR89) deployed in the USCB area during the CoMet mission (see Fig. 1 for respective positions) have been used to augment the model output (Wildmann et al., 2020). These data are available on a regular, continuous basis throughout the campaign period at 10 min time intervals with soundings typically reaching up to ∼2.5 km a.M.S.L depending on the atmospheric condition. Domains D1 and D2 are both nudged towards the Doppler soundings using the WRF-FDDA subsystem (Deng et al., 2008). Sensitivity of the model output on three key 150 parameters of the observational data assimilation subsystem, namely the radius of influence r xy , time window ∆ t and horizontal wind coefficient c uv were analyzed through numerous runs with the goal of finding an appropriate configuration. Figure 4 shows ensemble runs with varying r xy in comparison to interpolated NCEP GDAS/FNL and the actual lidar soundings for June 6th, 2018 at 0900 UTC. The shaded gray area beneath the orange-colored lidar soundings shows the maximum variability including soundings timed 20 min before and 20 min after the observations used. Figure 4 demonstrates that modeled data are in good agreement to observed Doppler soundings when using WRF-FDDA. It also shows discrepancies between NCEP GDAS/FNL driver data and observations in wind direction and more importantly on the wind speed in the lower troposphere and the PBL depth. To further enhance compatibility between model and observations, the WRFDA submodule (Barker et al., 2012) was used in 3DVar cycling mode similar to Liu et al. (2013) using the NCAR CV3 background error covariance (Barker et al., 2004). The required observational error covariances are taken from the measurement uncertainties.
The total emission estimate Φ in units kg s −1 follows from the scaled sum of the individual contributions ϕ i For use cases where no a-priori information on emitters is available, the above equations directly yield the best total emission estimate Φ. In the present case however, a-priori information on the emissions of the individual shafts can be used to compute a maximum a posteriori (MAP) solution x based on the observation vector y and the a-priori vector x a (Rodgers, 2000).
with later defined a-priori and observational error covariance matrices S a and S , respectively. The MAP solution can be found by solving for ∇ x J (x) = 0 and is given bŷ with the gain matrix By exploiting the averaging kernel A = GK the number of degrees of freedom for signal d s can be computed as This number describes the reduction in the normalized error on x introduced by the available observations and hence provides a measure for the improvement in knowledge of x, relative to the a-priori, due to the observations.

210
Here, the Non-Negative Least Squares (NNLS) algorithm (Lawson and Hanson, 1995) has been used to minimize the MAP cost function subject to the constraint x > 0. This constraint is equivalent to the absence of negative sources. The NNLS algorithm solves the constrained least squares problem by splitting into active and passive subsets, where active and passive refer to the state of the constraint. The algorithm subsequently solves the unconstrained least squares problem for the passive set.

Estimating total uncertainty
The outlined approach is based on assumptions, of which the most important ones are: a constant emission rate over the timescale of transport from the source to the aircraft, an appropriate atmospheric background vector b (used to compute CH 4 enhancements y = ρ−b from measured mass densities ρ), long lifetime of the species of interest, i.e. no chemical and physical removal on the timescale of a flight and the model being able to adequately represent the meteorological state variables. To 220 assess uncertainty on the retrieved emission rates, several variables have been selected as most influencing systematic error sources: wind speed, wind direction, PBL height, source dislocation and an error in sensed mole fractions, that is further intended to include an error due to chosen background. Individual contributions of these error sources to total uncertainty can be identified from ensemble model runs with systematically perturbed parameters.
In addition to derived systematic uncertainties, statistical errors related to the MAP fit are to be acknowledged. The statistical 225 uncertainty i in the retrieved parameters x i can be expressed in terms of the parameter covariance matrixŜ as The parameter covariance matrixŜ is computed from the m-by-n dimensional forward model K usinĝ σ 2 on its main diagonal inflated with a transport model error χ m obtained from ensemble runs with perturbed parameters (see Sect. 4.6).
The measurement uncertainties σ i were obtained via standard error propagation from the uncertainties associated with different instruments aboard the aircraft needed for the computation of the CH 4 mass density observations (see Eq. 10). Static 235 air temperature can be probed with an uncertainty of σ T = 0.15 K, static air pressure with σ p = 1 hPa and wind speed with σ u = 0.3 ms −1 (1s-1σ, Mallaun et al. (2015)). CH 4 mole fractions were sampled with a total uncertainty better than 1.85 ppb (1s-1σ). The a-priori error covariance matrix S a is a diagonal matrix with the squared a-priori uncertainties, estimated with 50 % of the nominal value on its main diagonal. Finally, the Jacobian ∇ x J (x) is the first derivative of the cost function with respect to the parameters x i . As such it describes the change in residuals introduced by a change in parameter x i . In situ E-PRTR a-priori Optimized model where m x denotes the total mass of the species of interest. The unit-less coefficient m x m −1 air = c x M x M −1 air is obtained from the sensed CH 4 mole fractions c x in units mol mol −1 .
Prior to conversion, the atmospheric background has to be subtracted from the observed c x . The choice of background is to some extent ambiguous, because there is no clear edge between background and in-plume sampling. This contributes to total 250 flux estimation uncertainty, as will be discussed later. Here, a piecewise linear interpolation between the outermost boundaries of each of the 4 flight legs (see Fig. 7) has been considered as the best guess of atmospheric background. The mean value of 20 samples has been used on both edges of each flight leg. Using this approach, latitudinal and longitudinal gradients in background CH 4 mole fractions are accounted for by using both edges of each flight leg. Vertical gradients in background CH 4 levels are accounted for by treating each constant-altitude flight leg separately.

255
From Fig. 6 a good overall match between model (solid black line) and in situ observations (gray) is apparent with a mean bias of 2.5 × 10 −9 kg m −3 and a root mean square error of 1.6 × 10 −8 kg m −3 . Some of the minor structure is not reproduced in detail by the model, which is expected due to the model's 3 km horizontal grid resolution. The reason for the discrepancy between model and the first few hundred observations in Fig. 6 becomes more obvious when looking at the 2D scene shown in Fig. 7. The left panel of Fig. 7 shows a cross section of the model output along the downwind wall including the in situ   Fig. 1. Figure 9 shows the corresponding time series of the measured and modeled CH 4 mass density as a function of flight time during the downwind wall phase with source coefficients x i already optimized. Alike for the morning flight a good overall match between model and in situ observations can be observed with a mean bias of 3 × 10 −10 kg m −3 and a root mean square error of 1 × 10 −8 kg m −3 .

285
The sensed mixing ratios are lower compared to the morning flight due to a further developed and hence more diluted boundary layer in the afternoon. The corresponding snapshot 2D scene is depicted in Fig. 8 The model based approach provides a unique advantage over established mass balance techniques in terms of spatial information, as it enables attributing sensed CH 4 mole fractions to remote sources at distances of tens to hundreds of kilometers.
The total emission estimate has been introduced in Sect. 4.2 as the sum over n sources ϕ i that are individually scaled with a coefficient x i . The emission rate Φ i corresponding to the i-th source is thus given by x i ϕ i . Here, the availability of data from both research flights on June 6th, 2018 was exploited to estimate emission rates Φ i for individual sources. By combining data 305 from both flights the number of degrees of freedom for signal increases to d s = 48 with 42 sources actively emitting. in key sources of uncertainty. Systematic uncertainties for each source are directly obtained from this ensemble run. Differences in estimated and reported (E-PRTR 2017) Φ i are evident. This is however expected due to the comparison of instantaneous emission estimates and yearly averages.

Uncertainty analysis 315
The influence of several variables on the total flux estimate Φ has been computed following Sect. 4.3. Figure 11 shows the influence of an error in wind speed (σ u = 0.9 ms −1 ), wind direction (σ d = 5 • ), PBL height (σ pbl = 100 m) and a source dislocation (σ sd = 1 km) to total uncertainty for the flights detailed in the previous section. An assumed error in sensed mole fractions (σ c = 10 ppb) is intended to include an error due to wrongly chosen background. The error on wind speed σ u is taken as the standard deviation of the difference between WRF modeled wind and non-assimilated in situ observations from the data 320 depicted in Fig. 5. The same holds for the wind direction. The difference between modeled data and observations should therefore reflect overall uncertainty in these variables. Two spiral-up soundings out of the PBL revealed a boundary layer height of 1150 m at 0937 UTC and 1300 m at 1145 UTC. Based on these two soundings the uncertainty on boundary layer depth is estimated with σ pbl = 100 m for the downwind wall phase between 1000 UTC and 1100 UTC. For this sensitivity analysis, the WRF fields were perturbed systematically during the FLEXPART read phase in "readwind.f90". The source dislocation was 325 implemented in the FLEXPART configuration file. It is evident from Fig. 11 that all selected error sources contribute on a similar level to total systematic uncertainty, which is ultimately computed as the standard deviation of the ensemble.
In addition to the derived systematic uncertainties, statistical errors related to the least squares fit have been computed following Sect. 4.3. Figure 12 depicts the Jacobian with respect to x i and the observations of the morning flight on June 6th,    Figure 11. Ensemble runs to assess uncertainty in the flux estimates derived using a model based approach. All selected error sources contribute to total uncertainty on a similar level.

Conclusions
A modified Aerodyne Dual QCLS instrument has been deployed aboard the DLR Cessna 208B in the context of the CoMet 1.0 campaign in early summer 2018 with the goal of estimating hard coal mine CH 4 emissions emanating from the USCB 340 area -Europe's largest coal extraction region. Intensive mining activities and the heavy industry spread around the city of Katowice lead to significant amounts of GHGs emitted into the atmosphere. The reported inventorial CH 4 emission rates for the entire USCB region amount to 720 kt yr −1 (EDGAR v4.3.2) and 448 kt yr −1 (E-PRTR 2017). The latter corresponds to 12.5 MtCO 2 -eq yr −1 using a CH 4 -GWP 100 =28 from the IPCC Fifth Assessment report (Pachauri et al., 2014). Assuming an average carbon content of 75 %, a net calorific value of 29 MJ kg −1 , an emission factor of 94 tCO 2 (TJ) −1 and the approxi- reported emissions are observed. This is expected due to several reasons: limited amount of measurements relative to the yearly averages provided in the inventories, wind directions do not differ by much between the two flights and the evenly distributed emissions among the ventilation shafts for each mining company. In general, the approach described herein delivers more information compared to the conventional mass balance, albeit at increased effort: wind lidars need to be deployed during the 370 measurement campaign, models need to be run, wind lidar data needs to be assimilated and inverse estimation techniques need to be applied. The additional possibility of remote source attribution however, coupled with the results obtained in Sect. 4.2 for the regional USCB anthropogenic CH 4 emissions make this approach a potent alternative to the mass balance technique.
Although retrieving estimates for individual emitters is not possible using only single flights, due to sparse data availability, the combination of two or more flights allows for exploiting different meteorological conditions to enhance confidence on facility 375 level estimates.