Observational constraints on methane emissions from Polish coal mines using a ground-based remote sensing network

. Given its abundant coal mining activities, the Upper Silesian Coal Basin (USCB) in southern Poland is one of the largest sources of anthropogenic methane (CH 4 ) emissions in Europe. Here, we report on CH 4 emission estimates for coal mine ventilation facilities in the USCB. Our estimates are driven by pairwise upwind– downwind observations of the column-average dry-air mole fractions of CH 4 (XCH 4 ) by a network of four portable, ground-based, sun-viewing Fourier transform spectrometers of the type EM27/SUN operated during the CoMet campaign in May–June 2018. The EM27/SUN instruments were deployed in the four cardinal directions around the USCB approximately 50 km from the center of the basin. We report on six case studies for which we inferred emissions by evaluating the mismatch between the observed downwind enhancements and simulations based on trajectory calculations releasing particles out of the ventilation shafts using the Lagrangian particle dispersion model FLEXPART. The latter was driven by wind ﬁelds calculated by WRF (Weather Research and Forecasting model) under assimilation of vertical wind proﬁle measurements of three co-deployed wind lidars. For emission estimation, we use a Phillips–Tikhonov regularization scheme with the L-curve criterion. Diagnosed by the emissions averaging kernels, we ﬁnd that, depending on the catchment area of the downwind measurements, our ad hoc network can resolve individual facilities or groups of ventilation facilities but that inspecting the emissions averaging kernels is essential to detect correlated estimates. Generally, our instantaneous emission estimates range between 80 and 133 kt CH 4 a − 1 for the southeastern part of the USCB


Introduction
The atmospheric abundance of methane (CH 4 ) has increased by a factor of 2.6 since pre-industrial times from roughly 720 ppb (parts per billion) to about 1879 ppb in 2020 (Dlugokencky, 2021), mainly driven by anthropogenic influences (e.g., Bousquet et al., 2006;Loulergue et al., 2008;Kirschke et al., 2013;IPCC, 2013;Nisbet et al., 2014;Conley et al., 2016;Schwietzke et al., 2016;Worden et al., 2017;Alvarez et al., 2018;Saunois et al., 2020;Hmiel et al., 2020). Roughly 20 % of the total global anthropogenic CH 4 emissions are caused by the fossil fuel industry (Bousquet et al., 2006;Schwietzke et al., 2016;Saunois et al., 2020) and an extensive source of CH 4 is hard coal mining. Poland is the largest hard coal producer in the European Union with the Upper Silesian Coal Basin (USCB) as one of the largest hard coal producing regions in Europe. Several bottom-up inventories report on the total CH 4 emissions for the USCB: according to the GESAPU database, the USCB emitted a total of 405 kt CH 4 in 2010 (Bun et al., 2019). The E-PRTR (http: //prtr.ec.europa.eu/, European Pollutant Release and Transfer Register, 2018) collects emission reports of every individual mine, yielding an aggregated total of 507 kt CH 4 a −1 for the USCB. Dreger (2021) reports hard coal mining emissions of 530 kt CH 4 for the USCB in 2018, and the Copernicus Atmosphere Monitoring Service regional emission inventory (CAMS-REG-GHG/AP) lists 632 kt CH 4 a −1 (Granier et al., 2019;Fiehn et al., 2020). EDGAR v4.3.2 (Emission Database for Global Atmospheric Research) gives emissions of 675 kt CH 4 a −1 (Janssens-Maenhout et al., 2017) for fuel exploitation emissions (EDGAR abbreviation PRO) for the USCB in the year 2012. EDGAR v6.0 (Crippa et al., 2020) states 454 kt CH 4 a −1 for the sector abbreviation PRO COAL for the USCB in 2018. The extrapolated seasonal breakdown provided by EDGAR for the relevant months May and June 2018 matches the total 2018 EDGAR emission estimate for the USCB.
In addition to the bottom-up inventories, top-down approaches have examined USCB emissions. During the CoMet mission (Carbon dioxide and Methane mission 2018, from 23 May to 12 June 2018), several ground-based instruments and aircraft measured the atmospheric CH 4 abundance in the USCB. In our precursor study (Luther et al., 2019), we used stop-and-go measurements of the columnaverage dry-air mole fractions of CH 4 (XCH 4 ) by a mobile, ground-based Fourier transform spectrometer (FTS) to evaluate the mining emissions of individual ventilation facilities and found similar emissions as suggested by the E-PRTR in-ventory. The total USCB emission estimates of Fiehn et al. (2020) and Kostinek et al. (2021), based on airborne in situ measurements, are in broad agreement with the E-PRTR data for single flights. Using airborne imager data, Krautwurst et al. (2021) found some discrepancies between their estimates and the E-PRTR inventory for small groups of ventilation facilities. The isotopic CH 4 composition was measured by Menoud et al. (2021) with ground-based in situ instruments. Swolkień (2020) discusses the short-term, shaft-wise CH 4 release in the USCB.
Given the magnitude of emissions and the range of estimates, CH 4 in the USCB warrants further investigation. Here, we report on CH 4 emission estimates derived from measurements of four stationary, sun-viewing FTSs of the type EM27/SUN arranged in a network-like pattern enclosing the USCB during the CoMet campaign activities. The setup largely mimics previous network deployments for quantifying urban greenhouse gas emissions in Berlin , Paris , St. Petersburg (Makarova et al., 2021), Munich (Dietrich et al., 2021), Indianapolis (Jones et al., 2021), and other places. Our four EM27/SUN instruments were positioned in the four cardinal directions at a distance of a few tens of kilometers to the center of the USCB. We calculate differences between pairs of upwind and downwind observations to determine enhancements ( XCH 4 ) in our XCH 4 records attributable to sources within the USCB. Further developing the model setup of Kostinek et al. (2021), we use the flexible Lagrangian particle dispersion model (FLEXPART) together with wind fields from the Weather Research and Forecasting model (WRF) constrained by three wind lidars to translate the observed XCH 4 enhancements into emission estimates for groups of coal mine ventilation shafts. To this end, we set up an inverse estimation scheme based on Phillips-Tikhonov regularization. This setup enables careful information content analysis while not relying on estimates of a priori emission uncertainties, which are often inaccessible. This work demonstrates emission estimation of a methane-emitting hotspot based on an FTS network combined with a Phillips-Tikhonov regularized inversion approach.
Our paper first summarizes the campaign setup (Sect. 2). Then, Sects. 3 and 4 detail the modeling and regularization methods. Section 5 reports on emission estimates for six case studies, and Sect. 6 discusses our results in terms of compatibility with other emission estimates and methodological strengths and weaknesses.

Campaign deployment and XCH 4 measurements
As part of the CoMet activities in May and June 2018, our campaign deployment of the EM27/SUN focused on the USCB, extending roughly 80 × 80 km 2 in the southwest of Poland. Figure 1, adopted from our precursor study on mobile measurements (Luther et al., 2019), illustrates the network pattern together with the locations of the most important coal mine ventilation shafts and the wind lidars. The four stationary EM27/SUN spectrometers are positioned roughly in the four cardinal directions around the center of the basin, ensuring that at least one instrument measures upwind and one downwind of the USCB for most wind situations. Due to mainly easterly wind conditions prevailing during our deployment, the eastern station (The Glade) functions as the upwind station for all case studies discussed here. The respective downwind stations are located in the south (Pustelnik) and west (Raciborz) of the USCB. The northern station (Za Miastem) was identified as neither an upwind nor downwind station for any of the cases.
The functioning of the EM27/SUN FTS and the data reduction techniques are described in detail by Gisi et al. (2012), Hase et al. (2015Hase et al. ( , 2016, and with a particular focus on our setup by Luther et al. (2019). Generally, the EM27/SUN instruments observe spectra of direct sunlight in the shortwave infrared spectral range from which the total column concentrations of CH 4 , carbon dioxide (CO 2 ), water vapor (H 2 O), molecular oxygen (O 2 ), carbon monoxide (CO), and other substances (Butz et al., 2017) can be retrieved. CH 4 and O 2 are of relevance here. All EM27/SUN spectrometers participating in the campaign successfully underwent the instrumental quality assurance tests required by the Collaborative Carbon Column Observing Network (COCCON) and presented in Frey et al. (2019) before field deployment.
We run the software package PROFFIT (Hase et al., 2004) to retrieve the column concentrations [O 2 ] and [CH 4 ] from the 7765 to 8005 and 5897 to 6145 cm −1 spectral windows, respectively, using absorption line parameters by Toon (2017) and Rothman et al. (2009). The respective columnaverage dry-air mole fractions of methane, XCH 4 , are calculated via normalization through [CH 4 ] [O 2 ] × 0.2094, where the atmospheric O 2 mole fraction is assumed to be 0.2094. The normalization is generally recommended to lessen spurious artifacts due to pressure and solar zenith angle (SZA) dependencies. The SZAs during our measurements did not exceed 56 • , which is within the range that does not require air-massdependent bias correction . Slightly deviating from the processing recipe  and in line with our precursor study (Luther et al., 2019), our CH 4 retrievals only scale the CH 4 concentrations in the layers below 1700 m a.g.l. where we expect elevated methane due to the localized sources at the ground. Further, we extract the a priori CH 4 profiles from a dedicated ECHAM/MESSy Atmospheric Chemistry simulation described by Jöckel et al. (2016) and Nickl et al. (2020).
For network deployments such as undertaken here, it is common practice to cross-calibrate the network nodes by side-by-side measurements Chen et al., 2016;Jones et al., 2021;Dietrich et al., 2021) in order to exclude spurious gradients when conducting upwinddownwind analyses. We calibrated the four instruments through side-by-side measurements on 23 and 26 May 2018 at the southern location Pustelnik. Figure 2 shows the raw and calibrated XCH 4 records, and Table 1 lists the respective calibration factors. These instrument-specific empirical calibration factors have been applied to the XCH 4 campaign data discussed in the following. After calibration, the mean absolute instrument-to-instrument difference is roughly 1 ppb, which is compatible with Frey et al. (2019), who state the instrument-to-instrument difference for an EM27/SUN ensemble as 0.8 ppb for XCH 4 . We estimate the precision for XCH 4 as 0.6 ppb for 1 min integration time from measurements sampled during the campaign in the rather variable methane field of the USCB, which is consistent with the precision stated by Chen et al. (2016). Based on measurements of the eastern and northern instruments (The Glade -E and Za Miastem -N) observed from 07:00 to 10:00 UTC on 28 May 2018 (upper panel Fig. 7), we calculated the precision as the standard deviation of the observations averaged between the two instruments during this period. We chose this period, since the two timelines are not affected by any strong methane sources in the vicinity.  In addition to the EM27/SUN network, we also operated three Leosphere Windcube 200S Doppler wind lidars (Vasiljević et al., 2016;Wildmann et al., 2018Wildmann et al., , 2020 marked as red stars in Fig. 1 and as detailed in Luther et al. (2019) and Kostinek et al. (2021). The measured wind profiles (10 min time interval) reach up to 4 km a.g.l. and are assimilated into the WRF simulations to improve the modeled wind fields as discussed in Sect. 3.

Dispersion modeling of methane
Our simulations of methane dispersion in the USCB are partitioned into two steps, largely adopting the basic setup reported by Kostinek et al. (2021): first (Sect. 3.1), the wind fields are modeled by a two-domain WRF setup including assimilation of the wind lidar observations. Second, the WRF wind fields drive the particle dispersion in the Lagrangian particle dispersion model FLEXPART (Sect. 3.2).

WRF wind fields
The WRF V4 (Skamarock et al., 2019) setup is driven by 3-hourly GFS data (NCEP, 2017) in two domains ( Fig. 3) focusing on central Europe and the USCB. The outer domain has a spatial resolution of roughly 15 km, and the inner domain has a spatial resolution of roughly 3 km. The simulations start at 00:00 UTC on the day of interest to allow for several hours of spin-up. Details are explained in Kostinek et al. (2021).
WRF has the possibility to assimilate observational data via four-dimensional data assimilation (FDDA), which we use for our wind lidar measurements. At a 10 min time interval, the wind profile observations are fed into the WRF calculations to constrain the simulated wind fields. The assimilation process is adjustable via several parameters, e.g., radius of influence r xy , time window t, and horizontal wind coefficient c uv , in the WRF input file directly. Kostinek et al. (2021) have identified a selection of parameter settings to obtain the best-guess parameter combinations for the same region, WRF domains, and time periods discussed here. There-fore, we adopt the same setup here and report on the overall WRF performance by comparing the simulations with the wind lidar observations. Figure 4 displays a comparison between modeled and observed wind speed (upper panel) and wind direction (lower panel) for altitude levels in the planetary boundary layer (PBL). PBL height is estimated by means of eddy dissipation rate gradients calculated from the wind lidar observations directly. The observation levels of the wind lidar do not match the WRF grid levels exactly. Thus, the comparison includes modeled and observed data if the WRF wind level is within the range of ±25 m around the level of observation. This represents all of the WRF wind levels but dismisses some of the wind lidar levels as the lidar data have a finer vertical resolution. The comparison is restricted to the three days of the case studies to be discussed below (28 May; 6, 7 June 2018). The prevailing wind directions were northeast to southeast, and prevailing wind speeds range from 2 to 9 ms −1 . The root mean square error (RMSE) for the wind speed comparison calculates to 1.7 ms −1 . When eliminating the top two levels, the wind direction RMSE is reduced from 38 to 26 • , indicating model inaccuracies towards the PBL top where wind shear effects are expected. Low wind speeds on 7 June possibly deteriorate the wind direction bias, as wind direction uncertainties are generally larger for low wind speeds. Further challenges for the wind direction estimation are the onset of convection during the observational morning hours with subsequent PBL rise and the calming winds towards the end of the observational days. In general we observed most outliers for the top comparison levels for wind direction, which could be related to a significant number of conspicuous low-wind-speed simulations and observations for these levels. This might be related to model uncertainties when estimating the PBL height, leading to misinterpretation of actual above-PBL observations with inside-PBL simulations and vice versa.

Lagrangian methane dispersion via FLEXPART
WRF wind fields drive the trajectory calculations in FLEX-PART. The model simulates trajectories for 50 000 particles for every USCB coal mining shaft reported by the E-PRTR and the CoMet database  with a total mass of 10 5 kgĊH 4 . The simulations do not consider background CH 4 . The model releases particles in a 10 m ×10 m ×10 m box on the ground. The modeling period starts at 00:10 UTC the day of interest and continues until 17:50 UTC, which results in 17.7 h simulation time. We chose the grid output option in FLEXPART with 100 × 100 boxes and a spatial resolution of roughly 1.3 km stacked in 24 layers up to 3 km altitude. The simulated XCH 4 measurements are the sums of all boxes above each pixel enclosing an EM27/SUN location. The 6 min FLEXPART output is interpolated to the observational time interval, which is generally one measurement per minute. After unit conversion, the sim- Comparison of WRF wind estimates and wind lidar observations for the observational period between 07:00 and 17:00 UTC for the three days of interest 28 May (blue), 6 June (orange), and 7 June (green). Shaded areas include roughly 80 % of the data points for ±2 ms −1 for wind speed in (a) and include roughly 50 % of wind direction measurements in a ±15 • range in (b). Wind direction simulations differ most from the observations for 7 June (green), a day with generally lower wind speeds than the other two discussed cases. This indicates wind information uncertainties when it comes to low wind speeds. For the wind direction comparison, note that values close to 0 and 360 • represent virtually the same wind direction but may introduce an error to the root mean square error (RMSE) calculations. ulated methane enhancement is compared to the measured upwind-downwind difference XCH 4 .
The FLEXPART simulations are iterated with slightly different meteorological parameters to provide an uncertainty analysis. There are seven ensemble runs: the CONTROL run with best-guess input, the WINDp5 run with +5 • wind direction change of the whole wind field, the WINDm5 run with −5 • wind direction change of the whole wind field, the SPEEDp06 run with the wind speed increased by 0.9 m s −1 , Figure 5. Lagrangian time lag sketch. Instrument A measures background methane (upwind). The methane enhancements induced by the emissions of the shaft are calculated by subtracting the background (upwind) measurements from the downwind observations of instruments B and C. Depending on wind speed, the air mass measured by instrument A will be at instrument B and C after different travel times.
the SPEEDm06 run with the wind speed decreased by 0.9 m s −1 , the PBLp100 run with the PBL height increased by 100 m, and the PBLm100 run with the PBL height decreased by 100 m. We use the same ensemble setup as discussed by Kostinek et al. (2021).
We further used the FLEXPART simulations for modeling the air mass travel time from the upwind to the downwind instruments (Fig. 5). To this end, we implemented a virtual methane source at the upwind measurement location at the beginning of each observational period. Then, we estimated the air mass travel time by recording the time required for the virtual tracer plume to cross the tangent point of the downwind measurement location, i.e., the point that was closest to the downwind location along the PBL-averaged wind direction. Given that the distance between upwind and downwind locations exceeds 50 km for some cases, estimated travel times range between 1.5 and 4.6 h. Thus, in order for the upwind measurement to be representative of the background conditions for the later downwind measurement, we consider the air mass travel time between the two measurements by subtracting the respective time-shifted upwind measurements when calculating XCH 4 . The travel time is calculated once for every case study beginning with the first upwind measurement, which is part of the inversion process, and ending when the virtual tracer reaches the respective downwind instrument. To account for possible background methane variability, we altered the simulated travel time by ±30 min and calculated the average methane enhancement for each individually time-shifted period. The difference between minimum and maximum averaged XCH 4 of these periods is included in the emission estimation routine to represent the error related to background methane uncertainties.

Emission estimates via Phillips-Tikhonov regularization
In order to estimate the shaft emissions from the mismatch between measured and simulated methane enhancements XCH 4 , we use a Phillips-Tikhonov inverse method. We set it up such that the state vector x consists of m dimensionless factors that scale the emissions of each coal mine ventilation shaft considered by the FLEXPART simulations. We assume the scaling factors to be constant for each day of measurement; i.e., we impose the assumption of the source strength being constant over the time of a day. FLEXPART is the forward model K (m × n) relating the emissions of the n shafts to m XCH 4 measurements. The measurement vector y contains the XCH 4 enhancements, observed at 1 min intervals, translated from units of parts per billion (ppb) into total mass column enhancements (kg m −2 ). The Phillips-Tikhonov inverse method then delivers the estimated state vector x λ by minimizing a two-term cost function consisting of a measurement term and an a priori term (Phillips, 1962;Tikhonov, 1963;Twomey, 1963): with S the error covariance matrix, λ the regularization parameter, W the weighting operator, x a the a priori state vector, and || · || 2 representing the L 2 norm. S contains the averaged standard deviation of the FLEXPART simulation ensemble summed in quadrature with the XCH 4 background variability and the measurement noise. Note that, for simplicity, we did not consider correlations in S . The estimated background error ranges between 0.3 and 2.2 ppb for the individual case studies based on the differences of minimum and maximum average XCH 4 calculated under consideration of ±30 min time shifts of the simulated methane travel times. The measurement noise amounts to 0.6 ppb, calculated as the standard deviation of the averaged measurements of The Glade (E) and Za Miastem (N) from 07:00 to 10:00 UTC on 28 May 2018 (see top panel in Fig. 7). W is a diagonal matrix with elements 1 x a,j , j = 1, . . ., m, which renders the second cost term dimensionless (Butz et al., 2012). Technically, we transform the Phillips-Tikhonov regularization problem into a plain least-squares fit using the definitions (Hansen and O'Leary, 1993;Hansen, 1999;Golub and Von Matt, 1997) which transforms Eq.
(1) to and is treatable by a standard least-squares solver (e.g., python module scipy.optimize.lsq_linear). The a priori information x a for each ventilation shaft is taken from the annual E-PRTR emission inventory updated by Gałkowski et al. (2021). The a priori generally guides the minimization process towards physically reasonable solutions if the inverse problem tends to be ill-posed, e.g., when there is insufficient measurement information on some of the state vector elements. However, the solution is dependent on the regularization parameter λ, which has to be found by trading the propagation of measurement errors against influence of the a priori. Here, we use the L-curve criterion to determine the regularization parameter λ for each individual case study (e.g., Hansen, 1999). The L-curve is a graphical representation of the mismatch between measurements and simulations K x λ − y 2 plotted against the norm of the state vector x λ 2 evaluated for a range of λ (see Fig. 9). The plot typically looks like an L. For small λ, the measurement term dominates the cost function, and the estimate x λ becomes noisy and drives large deviations from the a priori. For large λ, the a priori term dominates the cost function, and the estimate x λ ignores the measurements and produces a large norm of the measurement term. The corner of the L indicates a reasonable regularization parameter for the given minimization problem and is graphically chosen. In addition, we found that the shape of the L-curve is sensitive to forward model errors, e.g., when errors in the FLEXPART trajectories and the driving wind fields suggest a spurious, erroneous link between emissions and methane enhancements. Obvious distortions of the L-curve shape are used as a criterion to filter out periods when the forward model does not represent the actual dispersion conditions.
Besides the L-curve, the regularized inversion approach holds another diagnostic measure: the emissions averaging kernel matrix A λ with dimensions m × m. It is defined via the gain matrix G λ (Rodgers, 2000;Butz et al., 2012;Borsdorff et al., 2014).
The averaging kernel matrix, for a given regularization strength λ, diagnoses how information propagates from the true and a priori states, x true and x a , into the emission estimate.
The rows of A λ are called the averaging kernels, quantifying how an estimated state vector element calculates from the other state elements and what portion comes from the prior. For our purposes, the averaging kernels quantify how the emission estimate for a ventilation shaft is affected by the neighboring shafts and whether there is sufficient measurement information. In the perfect case, the emissions averaging kernel is unity for the shaft under consideration and zero for all other shafts, indicating that the shaft can be perfectly resolved and discriminated from neighboring sources and that it is well constrained by measurement information.
In reality, groups of neighboring sources and sources behind each other along the trajectory are not resolvable, and some shafts only marginally affect our measurements, implying broader and smaller emissions averaging kernels. The errors due to measurement noise, background methane variability, and atmospheric transport incorporated in S are propagated into the a posteriori error covariance for the emission estimates via where we report the square root of the diagonal as the error bars of the shaft-wise emission estimates and the square root of the sum of the entire covariance matrix as the error of the total emissions aggregated over all shafts. For the case studies discussed in Sect. 5, the emission errors due to measurement noise range between 0.62 and 4.46 kt a −1 , which is small compared to the errors introduced by the dispersion modeling ensemble that range between 27 and 143 kt a −1 . Errors related to background methane variability introduced by a ±30 min time shift of the air mass travel time range between 0.83 and 8.5 kt a −1 .

Case studies
We report on six case studies on three different days of the CoMet campaign. For the days 28 May and 6 and 7 June 2018, Fig. 6 illustrates typical FLEXPART trajectories of air masses around midday dispersing out of the USCB coal mine ventilation shafts. For all cases, easterly winds led to the southern station Pustelnik being influenced by a few southern shafts (red trajectories) and the western station Raciborz being influenced by many shafts in various parts of the basin (blue trajectories). The eastern station The Glade provides the background measurements, and the northern station Za Miastem was not used here since The Glade (E) was the better background station given the prevailing easterly winds. Figure 7 depicts the corresponding XCH 4 measurements for all stations, indeed pointing at significantly elevated concentrations at the downwind sites typically amounting to XCH 4 on the order of 10 ppb with some diurnal and day-today variability. For 28 May, the maxima during the morning hours at Pustelnik (S) and Raciborz (W) are most likely connected to nighttime methane accumulation and subsequent transport with rising convection and mixing during the morning hours. When considering the Lagrangian travel time of air masses, these morning hours are excluded from the period of investigation since we have to wait for the air masses from the upwind site The Glade (E) to arrive downwind at Pustelnik (S) and Raciborz (W). Therefore, the period of investigation starts later than the Pustelnik (S) and Raciborz (W) measurement records. For Pustelnik (S), travel times ranged between 1.3 and 2 h and for Raciborz (W) between 3 and 5.5 h; 7 June had the longest travel times linked to the slow wind speeds on that day (see Fig. 4). In particular, for 7 June, the consideration of the travel time implies that only a small fraction of the measurements at the end of the day is considered for the emission estimates. The stations The Glade (E) and Za Miastem (N) show roughly consistent background XCH 4 for 28 May and 6 June, but, on 7 June, significant differences indicate that the background concentration field was not homogeneous. Nonetheless, the FLEXPART trajectories point to The Glade (E) being representative of the background conditions for the downwind stations Raciborz (W) and Pustelnik (S). Following the procedures described in Sects. 3.2 and 4, we calculate observed XCH 4 from the XCH 4 observations by subtracting the background record of The Glade (E) from the Pustelnik (S) and Raciborz (W) records, shifting the time series by the calculated Lagrangian travel times. Then, we estimate the shaft-wise emissions by the Phillips-Tikhonov technique.

Southern station Pustelnik
Focusing on the Pustelnik (S) records first, Fig. 8 illustrates the measured and simulated XCH 4 enhancements (translated into units of kg m −2 ) and the time periods used for estimating emission rates (blue shading). Generally, FLEX-PART simulations using the a priori emissions based on E-PRTR underestimate the enhancements substantially for all cases. The FLEXPART simulations with optimized emission rates fit the measurements well. Residual discrepancies range between 2.35 × 10 −8 and 4.56 × 10 −8 kg m −2 (2.62 × 10 −8 and 4.72 × 10 −8 kg m −2 ) in terms of mean bias (root mean square error). The regularization strength for optimizing the emissions was determined via the L-curves depicted in Fig. 9, for which the L-shape is clearly recognizable. The corner of the L-curve and the corresponding regularization parameter λ is identified by visual inspection. Beside excluding data at the start of the daily time series due to the travel time, we also exclude data that we diagnosed to be affected by systematic forward model errors. The periods before 10:00 and after 13:30 UTC on 28 May (Fig. 8, upper panel) are an illustrative example. For these periods, we find substantial deviations between the measured and modeled time series that the optimization cannot resolve by adjusting emission rates. We expect that such patterns originate from systematic errors in the wind fields that drive FLEXPART. If we include these periods in our optimization scheme, we end up with distorted L-curves (Fig. 9, upper panel, grey dashed line). Thus, inspection of the L-curve provides us with a diagnostic tool to identify periods that are affected by forward model errors, which we then exclude. The other excluded periods are around noon on 6 June and the beginning of 7 June (faint colors in Fig. 8).
After selecting the data periods well represented by the FLEXPART simulations and choosing a suitable regulariza- Figure 9. The L-curves for the three Pustelnik (S) case studies. Panel (a) additionally depicts the L-curve of the same case study (dashed grey line) but under consideration of the full data set including the morning and afternoon hours, which suffered from forward model errors and are omitted in the final analysis (blue curve). The regularization parameters, λ, range from 0 to 100. The red star marker depicts the respective λ used for the emission estimation. tion parameter, we estimate the compatible shaft-wise emission rates illustrated in Fig. 10. Before evaluating the estimates, it is mandatory to inspect the emissions averaging kernels (upper sub-panels in Fig. 10), which encode the information on whether individual shafts or groups of shafts can be resolved. For 28 May (upper panel in Fig. 10), the emissions averaging kernels are greater than 0.7 and well positioned for the two Silesia shafts, while, for the Brzeszcze shafts, the emissions averaging kernels are small, indicating that there is very little measurement sensitivity to these shafts. Low sensitivity might be caused by air mass trajectories going over these shafts only for a short period of time. For the Silesia shafts, the estimated emission rates are substantially greater than the E-PRTR-based prior. On 6 June (middle panel in Fig. 10), the emissions averaging kernels indicate substantial sensitivity to Silesia V and Brzeszcze II, VI, and IX, as well as some sensitivity to Brzeszcze IV, but the emissions averaging kernels for Brzeszcze II and VI are double-peaked at each of the shafts, indicating the shafts cannot be resolved individually, e.g., because trajectories went over both shafts most of the time. A similar double-peaked emissions averaging kernel occurs for Silesia V and Brzeszcze VI. The shaft that can be best resolved is Brzeszcze IX. Generally, the opti- Figure 10. Shaft-wise emission estimates (lower sub-panels) for the three case studies at the southern station Pustelnik and the corresponding emissions averaging kernels (upper sub-panels). Colors of the shaft-wise emission estimates resemble colors of the emissions averaging kernel. Error bars contain atmospheric variability and observational uncertainty (measurement noise and background variability). Grey bars illustrate the a priori emission estimates. mized emissions are again significantly greater than the priors. For Brzeszcze VI, the prior even indicates zero emissions, while our estimate exceeds 10 kt a −1 . On 7 June (middle panel in Fig. 10), sensitivity is considerable for Silesia I and V, Brzeszcze IV, and to a lesser degree Brzeszcze II, but generally the emissions averaging kernels are not singlepeaked but distributed among several shafts, which fits the rather unstable wind conditions on that day. It is noteworthy that we find small emissions for Silesia I on 7 June, while for the other days Silesia I showed emissions of 15-20 kt a −1 . This might point at day-to-day variability of the ventilation. Aggregating the emissions over all shafts (Table 2), we find that on all three days the optimized emissions are substantially greater than the prior. On 28 May and 7 June, the total optimized and prior emissions are compatible within the error bars, but on 6 June, our estimates indicate emissions twice as high as the prior: 133 ± 31 kt a −1 compared to 63 kt a −1 . The difference of the total emission estimates from the other days is largely due to better sensitivity to all the Brzeszcze shafts and to the larger respective emission estimates.

Western station Raciborz
In contrast to the case studies for Pustelnik (S), the FLEX-PART simulations indicate that the western station Raciborz is influenced by a large and varying number of shafts (between 30 and 50 shafts for the three days discussed here). Figures 11 and 12 show the fits to the XCH 4 observations for the prior and optimized emission estimates and the Lcurves for the selection of the regularization parameter, respectively. For 6 June, we again find that the FLEXPART simulations could not represent our measurements at the beginning and end of the time series, and therefore we excluded these periods based on visual inspection of distortions of the L-curve. Overall, after optimization the simulated XCH 4 records fit the observations well, while simulations with prior emissions show substantially smaller-thanobserved XCH 4 . The residual differences between simulations with optimized emissions and observations range between 7.13×10 −8 and 1.45×10 −7 kg m −2 (7.49×10 −8 and 1.46×10 −7 kg m −2 ) in terms of mean bias (root mean square error). Figure 13 illustrates the shaft-wise emission estimates and the corresponding emissions averaging kernels (upper subpanels). Due to the large number of shafts involved, the picture is less clear than for the Pustelnik (S) station. The emissions averaging kernels show that it is typically groups of shafts to which the measurements have good sensitivity but that, mostly, individual shafts cannot be resolved. Overall, the optimization points at emissions greater than the prior. When aggregated over all the shafts (Table 2), the optimized emissions are a factor of 2.4 (28 May), 1.7 (6 June), and 2.4 (7 June) greater than the prior emissions, and the differences are greater than the error bars estimated from measurement errors and the ensemble simulations representing transport variability. On 28 May, our estimates indicate large contributions from ventilation associated with the southwestern part of the USCB. On 6 June, the more northern and central mines show emissions greater than the prior, which is zero for some of the mines (e.g., Centrum Witczak/Staszic and Julian II/Rozbark Barbara). On 7 June, the number of contributing shafts is less than for the other days, and the estimates point at large emissions from ventilation of Ziemowit, Janina, and Chwałowicze V.

Discussion and conclusion
We estimated CH 4 emissions from coal mine ventilation facilities in the USCB during the CoMet campaign in May-June 2018. To this end, we deployed four EM27/SUN spectrometers in a network configuration enclosing the USCB. Combining pairs of upwind-downwind XCH 4 observations  with trajectory simulations by the Lagrangian particle dispersion model FLEXPART, we attributed the observed XCH 4 enhancements to the emission rates of the facilities. The trajectory calculations were driven by wind fields simulated by WRF under assimilation of wind measurements by three wind lidars distributed throughout the region. The trajectory calculations also enabled us to take into account the travel time of air masses between the upwind and downwind stations. For emission estimation, we used a Phillips-Tikhonov regularization approach together with the L-curve criterion for selecting a well-suited regularization parameter. The regularization scheme comes with the emissions averaging kernel diagnostics that allow for quantifying which facilities or groups of facilities the measurements are sensitive to. For estimating errors, we constructed an ensemble of seven FLEX-PART runs, each with slightly perturbed atmospheric parameters (wind speed, wind direction, PBL height). Upscaling our emission estimates to annual emissions, we find higher emission rates than listed by E-PRTR. Other studies (Luther et al., 2019;Kostinek et al., 2021;Fiehn et al., 2020) have shown better agreement between their instantaneous estimates and the E-PRTR inventory, while, in our study, only two out of six cases are compatible with the E-PRTR inventory within the error range (see Table 2). The other four cases suggest 1.7 to 2.4 times higher emissions than reported by the E-PRTR. The comparisons with E-PRTR are only of an illustrative nature. Depending on the under-ground mining activities, the emissions from the mining process are highly variable in time, and thus our measurements are certainly not representative of the annual total reported by E-PRTR. In order to constrain the latter, a permanent observatory network would need to be operated with reasonably dense sampling throughout the year.
Our emissions estimates for the totals of all contributing facilities show errors in the range between 23 % and 36 %. Measurement errors (0.62 to 4.46 kt a −1 ) and backgroundvariability-induced errors (0.83 to 8.5 kt a −1 ) are negligibly small compared to errors induced by uncertainties of the wind fields (27 to 143 kt a −1 ). As hinted at by the emissions averaging kernels, estimates for individual shafts are correlated, and these correlations were taken into account when calculating the aggregated total emissions. For all cases, these background concentrations were taken from the measurements conducted at the eastern station The Glade since the trajectory calculations showed that it was more representative for the background than the northern station Za Figure 13. Shaft-wise emission estimates (lower sub-panels) for the three case studies at the western station Raciborz and the corresponding emissions averaging kernels (upper sub-panels). Colors of the shaft-wise emission estimates resemble colors of the emissions averaging kernel. Error bars contain atmospheric variability and observational uncertainty (measurement noise and background variability). Grey bars illustrate the a priori emission estimates.
Miastem. While for 28 May and 6 June, the two stations show similar concentrations (Fig. 7), The Glade (E) records higher XCH 4 than Za Miastem (N) on 7 June. Thus, the background concentration field has some spatial variability, which might be connected to remote sources such as the Krakow urban region affecting The Glade (E) but not Za Miastem (N) . In the future, we will aim at estimating errors due to spatial background variability by running a model that includes all known sources in a larger area around the target region or by deploying more spectrometers.
Our inverse method was set up such that the parameters to be estimated included all the individual CH 4 ventilation shafts for which the trajectory calculations indicated contributions to the measurements. The problem requires regularization since CH 4 plumes might mask each other (e.g., shafts located behind another shaft along the trajectory), since the detected contributions might have occurred only for a short period, and since, due to errors in the wind fields, the trajectories might indicate contributions when there were actually none or vice versa (e.g., for trajectories barely hitting the downwind station). The Phillips-Tikhonov regularization with L-curve criterion provides useful tools to diagnose the robustness and information content of the problem. The Lcurve shows distortions for episodes when the FLEXPART trajectories cannot represent the true link between emissions and downwind concentrations. Thus, inspecting the L-curve allows for excluding episodes affected by these simulation errors. The emissions averaging kernels are useful to identify the shafts or groups of shafts for which the emissions can be reliably estimated from the observations. For our configurations, we found that we can resolve emissions of individual shafts for the Pustelnik (S) downwind station where only a few comparatively close facilities contribute to the XCH 4 enhancements. For the Raciborz (W) station with a much larger catchment area, we can typically resolve groups of shafts along the main wind direction, but disentangling individual shafts is not possible when they are close to each other or behind each other along the air mass trajectory. Nonetheless, we argue it is useful to set up the problem in terms of individual shafts instead of a priori aggregating shaft clusters since inspecting the emissions averaging kernels provides a tool to check which shafts can be resolved and since the aggregated total emissions and their errors can be calculated a posteriori.
CH 4 emissions from sources other than coal mining are not considered in our FLEXPART simulations; i.e., the pairwise XCH 4 enhancements are assumed to be caused only by the ventilation shafts. Kostinek et al. (2021) have inspected EDGAR v4.3.2 for CH 4 emissions from sources that are not related to coal mining activities such as landfills, wastewater treatment, and livestock in the USCB. They found that these other sources amount to roughly 14 % of the USCB total CH 4 emissions and that the larger contributions stem form the northwest of the basin, which our measurements are not sensitive to. The CoMet inventory  lists annual emission estimates for 4 of 16 landfills relevant for our case studies, which in total amount to 0.97 kt a −1 . We expect similar emission estimates for the other 12 landfills reported with undefined emissions. Thus, we assume that contributions from sources other than coal mining are small compared to our error bars. We did not consider chemical reactions such as oxidation of methane with OH, since the timescales relevant for this study are too short for oxidation to have a significant influence on the emission estimations.
Overall, our study shows that deploying sun-viewing spectrometers in an ad hoc network configuration around the USCB allows for estimating CH 4 emissions from coal mine ventilation facilities with some resolution for individual facilities and groups of them, depending on the deployment location. Given that the errors are dominated by uncertainties in the wind fields driving the trajectory calculations, it is essential to validate model winds by local wind data or, as in our case, to assimilate local wind lidar measurements. In the view of developing such networks toward a monitoring capacity, priority should be put on making the networks permanent for better temporal representativeness of the observations and on making the networks denser in order to gain sensitivity to more shafts and to better quantify spatial variability in the background concentrations. In addition to the monitoring aspect, permanent networks could also validate satellite missions, and the different ground-based and spaceborne emission estimation approaches could be consolidated. Figure A1. Calibrated XCH 4 measurements from 26 May. The instruments for The Glade (E) and Pustelnik (S) measured side-byside at the station Pustelnik. The other two instruments measured at their respective campaign locations. The calibration of The Glade (E) towards the other instruments is based on these measurements. Data availability. WRF v4.0 can be downloaded from https: //www.mmm.ucar.edu/weather-research-and-forecasting-model (University  Corporation  for  Atmospheric  Research,  2020).
Author contributions. AL and AB wrote the paper. AL, JK, RK, LS, SD, MS, AF, AD, and DD operated the EM27/SUN spectrometers in the field during the campaign and collected and shared the data. NW provided the wind lidar data. AB, FH, MF, JC, and FD supported preparations for the measurement campaign, contributed to the spectral retrievals, and assisted with data postprocessing. JN and JS provided detailed information about coal mining and CH 4 ventilation and functioned as local advisers. CK and JK provided the WRF and FLEXPART framework. SV contributed to the discussion. AB and AR developed the research question.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement.
This article is part of the special issue "CoMet: a mission to improve our understanding and to better quantify the carbon dioxide and methane cycles". It is not associated with a conference.