Systematic detection of local CH 4 emissions anomalies combining satellite measurements and high-resolution forecasts

In this study we present a novel monitoring methodology to detect local CH4 concentration anomalies worldwide that are related to rapidly changing anthropogenic emissions that significantly contribute to the CH4 atmospheric budget. The method uses high resolution (7 km x 7 km) retrievals of total column CH4 from the Tropospheric Monitoring Instrument (TROPOMI) onboard the Sentinel 5 Precursor satellite. Observations are combined with high resolution CH4 forecasts (~9km) produced by the Copernicus Atmosphere Monitoring Service (CAMS) to provide departures (observations minus forecasts) 15 close to the native satellite resolution at appropriate time. Investigating the departures is an effective way to link satellite measurements and emission inventory data in a quantitative manner. We perform filtering on the departures to remove the large-scale biases on both forecasts and satellite observations. We then use a simple classification on the filtered departures to detect anomalies and plumes coming from CAMS emissions that are missing (e.g. pipeline or facility leaks), under-reported or over-reported (e.g. depleted drilling fields). Additionally, the classification helps to detect local satellite retrieval errors due 20 to land surface albedo issues.


Introduction
Atmospheric methane (CH4) is the second most important anthropogenic greenhouse gas after carbon dioxide and contributes significantly to changes in radiative forcing and climate change. CH4 is estimated to account for at least a quarter of the present-day warming (Myhre et al., 2013) and has a near-term global warming potential that is 84 times larger than CO2 25 per unit mass (IPCC 2013). There are numerous natural and anthropogenic CH4 sources, which vary in location and areal extent. The anthropogenic emissions related such as oil and gas production and coal mining and biomass burning tends to be geographically localised, e.g. over a plant facility, a pipeline or a field of extraction. Methane emissions however related to biological fluxes such as livestock, landfills and rice fields which can also be either geographically localised over narrow areas https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. generate initial conditions for high-resolution forecasts at about 10km (Agustí-Panareda et al., 2019). The lack of up-to-date emission inventories will impact and likely degrade simulated CH4 concentrations in areas where the local contribution of anthropogenic emissions is significant.
Many studies have shown the rapidly changing and event-based nature of CH4 anthropogenic emissions, especially in the case of identifying the location of 'super-emitter' point source locations. Conley et al. (2016) used aircraft measurements to 50 characterise a blowout of a well connected to the Aliso Canyon gas storage facility in California from October 2015 to February 2016. Pandley et al. (2019) showcased detection of large methane emission from a gas well blowout in Ohio during February to March 2018 using satellite measurements. More recently, Varon et al. (2019) detected an anomalously large CH4 source using a combination of satellite instruments over Central Asia (western Turkmenistan) associated with a gas compression station. Those types of suddenly occurring CH4 emissions cannot be or are not reported/detected in time to be included in the 55 bottom-up inventories but are seen from space. Other studies showed the capability of satellite measurements to detect CH4 emissions related to extensive drilling and fracking areas. Kort et al. (2014) identified a large methane anomaly over the Four Corners region of the USA and more recently de Gouw et al. (2020) showed satellite detection of large and extended enhancements in the San Juan, Uintah and Permian basin in the USA. While these satellite-based studies focused on specific events and locations, none of them systematically detected such anomalies at global scale, nor did they provide a method to 60 do so.
Systematic detection of large point sources of anthropogenic CH4 emissions using a combination of satellite observations and modelling could enable rapid action to reduce emissions from the oil and gas sectors. Two recent developments allow for https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. systematic detection of unreported CH4 atmospheric anomalies linked to small scale and point sources emissions. Firstly, newly available high resolution (7 km x 7 km) satellite observations from the TROPOspheric Monitoring Instrument 65 (TROPOMI, Veefkind et al., 2012) on board the Sentinel-5p platform. Secondly, improved real-time forecasting at highresolution (~9 km) provided by CAMS (Agustí-Panareda et al., 2019). In this paper we present a novel methodology to routinely compare the satellite observations with the model forecasts in order to perform a systematic detection of atmospheric CH4 anomalies related to emission changes from small scale and point sources emissions that are not reported or lack timely update. The paper is organised as follows: Section 2 describes the setup that includes the TROPOMI observations, the 70 forecasting and monitoring configurations, Section 3 presents the detection method, Section 4 discusses several case studies followed by conclusions where we discuss briefly the benefit of our approach with coarse resolution inverse modelling.

TROPOMI CH4 observations
The TROPOMI (Veefkind et al., 2012) instrument was launched 13 October 2017 onboard the Sentinel-5 Precursor 75 satellite, a low earth orbiter with a Sun-synchronous orbit that overpasses at 13:30 local solar time. Currently operational since the end of April 2018, the instrument is an imaging spectrometer with a wide spectral range: ultraviolet, visible, near-infrared and shortwave infrared. This allows TROPOMI to measure a variety of atmospheric chemical species such as: ozone, nitrogen dioxide, carbon monoxide, sulphur dioxide, formaldehyde, aerosol and methane (Hu et al., 2018). Current CH4 observations, which are available for the inner two thirds of the swath and only over land, are vertically integrated columns sensitive to the 80 troposphere (surface to 200hPa). With a swath of around 1,750km (normally 2,600 km) wide from the along track position and a ground pixel size of 7 km x 7 km, TROPOMI CH4 data can provide near global daily coverage at high horizontal resolution over land but is limited by cloud cover and retrieval quality. In this study, we use the bias corrected version of the product and we apply the most stringent quality flagging possible, selecting only pixels that have the qa_value = 1.0 (see Product Readme Methane V01.03.02, https://sentinel.esa.int/documents/247904/3541451/Sentinel-5P-Methane-Product-85 Readme-File). Figure 1 illustrates the CH4 satellite observation coverage that TROPOMI provides over a year, a month and a day.
The measurements show clear geographical variation of the CH4 column-averaged dry-air mixing ratios (XCH4) that are driven by the atmospheric transport but most importantly by the spatial and temporal variability of the surface fluxes and emissions variations. Figure 2 shows the 2019 annual average zoomed over the Middle East region and the western USA 90 regions. Over these regions, spatial variability results in XCH4 enhancements of up to 50 ppb over emission hotspots. Differences in the average concentrations from region to region are also significant, from approximately 1825 ppb over the USA to 1875 ppb over the Middle East. The strong local enhancements are an indication of strong local surface fluxes and emissions of CH4 from oil and gas activities, mining, agriculture or wetlands. XCH4 retrievals can also be prone to some systematic residual errors especially related to surface albedo (Hasekamp et al., 2019). De Gouw et al. (2020) for instance mentioned the possibility of retrieval biases due to low surface albedo in the short wave infrared spectral bands in the winter. Such retrieval biases, even though mostly reduced by the bias corrected product, need further investigation and are outside the scope of this paper. Nevertheless, the TROPOMI data are sufficiently accurate to show local enhancements linked (but not limited) to oil and gas production. We show in section 3 how to isolate these small-scale signals of interest and how to remove the contribution of large-scale biases. 100

CAMS high-resolution CH4 forecasting suite
In this study we use the ECMWF Integrated Forecasting System (IFS), which is used in different configurations for the operational Numerical Weather Prediction (NWP) system as well as for the Copernicus Atmosphere Monitoring Service (CAMS) atmospheric composition analyses and forecasts (e.g. Flemming et al, 2015). As part of the CAMS greenhouse gases services, the IFS is used to provide 5 days CO2 and CH4 forecasts (Agustí-Panareda et al., 2019) jointly with other species 105 relevant for air-quality (Flemming et al., 2015).
The IFS model cycle used in this paper is CY45R1 and is run routinely with a TCo1279 horizontal resolution which is a cubic octahedral reduced Gaussian grid at approximately 9km (Holm et al., 2016) with 137 vertical levels from the surface to 0.01hPa. Details about the transport and meteorological configuration can be found in Agustí-Panareda et al. (2019). The CAMS greenhouse gases (GHG) operational suite is composed of an analysis and forecasts at medium and high resolution 110 (see Fig. 3). The analysis is based on the IFS 4D-Var assimilation system which was adapted to assimilate retrieved columnsaveraged mole fractions of CO2 and CH4 together with all the operational meteorological observations (Engelen et al., 2009, Massart et al., 2014. The analyses are produced every 12hours (00:00UTC and 12:00UTC). A 4-day forecast is then issued daily after the 00:00UTC analysis on a TCo399, a cubic octahedral grid corresponding to approximately 25 km x 25 km with the same 137 model level configuration. Two satellite observation streams are currently assimilated, the Infrared 115 Atmospheric Sounding Interferometer (IASI) for CH4 on the MetOp satellites and the Thermal And Near-infrared Sensor for carbon Observations (TANSO) on the GOSAT satellite for both CO2 and CH4 (see Massart et al. (2014) for further details).
The processing and acquisition of the level 2 data in 2019 provided the satellite XCH4 data 4 days behind real time. The highresolution forecast is then coupled to the analysis experiment by merging the 4-day lower resolution forecast from the CO2 and CH4 analysis with the previous 1-day high resolution forecast ( Fig. 3) in order to preserve the fine-scale features of the 120 high-resolution forecast. Additionally, the high-resolution forecast coming from the operational NWP runs is used to reset the initial meteorological conditions in order to ensure the best possible accuracy of the transport. In this paper we will focus on using the CH4 forecasts at high-resolution coming from the setup described above. The high-resolution forecasts are run on a TCo1279 L137 grid of approximately 9 km x 9 km for a 5-day period and are initialized approximately 4 hours behind realtime every day from 00:00UTC. 125 Both high resolution forecasts and analysis use prescribed CH4 surface fluxes. The anthropogenic emissions including fossil fuel emissions, agriculture and landfill/waste emissions are from the annual EDGARv4.2FT2010 data set (Olivier and https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.

Monitoring suite
To monitor and compare the TROPOMI XCH4 retrievals with the IFS CH4 9 km forecasts we re-use a part of the IFS 145 assimilation system in a so-called monitoring mode. The system recomputes a high-resolution trajectory at 9 km initialised from the forecasts over a 12-hour monitoring window to calculate so-called first guess departures (difference between the observation and the model forecast) with the observations at the appropriate time. At each observation location the departure can be written as follows, where is the departure, the observation, the observation operator, the model integration or trajectory and ! the initial CH4 condition at the beginning of the monitoring window. If we inject the retrieval equation (Rodgers, 2000) the departure becomes, where " is the true CH4 concentration state (which is never exactly known), is the averaging kernel matrix which represents the sensitivity of the retrieval on the vertical profile with respect to the true state, the identity matrix, # the apriori information used in the retrieval and the retrieval error term. The equation then simplifies to, which is the difference between the true state and the forecast smoothed by the averaging kernel function plus the retrieval 160 error term. Those departure values are thus strongly dependent on the averaging kernel function shape. For the TROPOMI XCH4 retrievals the mean averaging kernel function shows a homogenous sensitivity to the entire troposphere up to 200hPa where the sensitivity decreases in the stratosphere. The averaging kernel function is not very variable between pixels or between different regions of the globe (not shown). Figure 6 shows the departures over various time scales (yearly, monthly and daily) for the global domain. Overall the departures (observation minus forecast) show a global positive bias of around 25 165 ppb (meaning observation values are above the model values) which could be attributed to model biases (Ramonet et al., 2019) and/or observation biases (Langerock et al., 2019). Ramonet et al. (2019) compared the CAMS CH4 forecasts with independent total column data. Results showed that the forecasts continuously underestimate the CH4 total columns by 5-20 ppb. Langerock et al. (2019) showed that the averaged total column bias for the TROPOMI CH4 retrievals bias is -0.32% (i.e. around -5ppb) but with respect to ground-based measurements. of large-scale error seen in the departures in not fully understood yet and is beyond the scope of this paper, although an understanding of these biases is crucial to further improve the quality of the CAMS CH4 forecasts and TROPOMI retrievals.
At finer scales, structures are seen on the yearly average comparison and become more evident on the monthly timescales.
Local differences are even stronger on a daily basis but recognising fine scale structures is challenging due to the lack of daily coverage. For those reasons a spatial filtering and temporal averaging of the departures is performed to extract and use the 180 small-scale features seen in the departures.

Filtering the signal
To remove the large-scale features seen in the departures we have implemented a high pass Gaussian filtering. The filter uses a convolution of a 2D Gaussian kernel on a given averaged and binned departure field. In this study we use a 0.1° latitude-185 longitude binning. Due to ocean, cloud cover and quality control flagging a number of bins of the departure will show missing values that will jeopardize the convolution. This problem is solved technically by creating two auxiliary matrices that have missing values replaced by 0. The two auxiliary matrices are then defined as where $ and is the average departure in the given bin, the threshold of minimum number of observations in a given bin.
In this study, we have chosen = 2 in order to avoid smoothing with very isolated pixels that can be faulty but also keep as much data as possible. Replacing the missing values by zeros in introduces an error after convolving (inducing low values 195 due to smoothing out with zeros) in the filtered departures %& . This can be compensated by applying the same Gaussian filter on a matrix representing the selected bins for filtering (where number of counts are above ) and using the ratio of the two filtered matrices to compensate for the missing value errors. Then a high pass filtering on a given observation space field (here departures) can be formulated as follows, where ( ) is a 2D Gaussian kernel function with a length scale. The same filtering is also applied on the observation values and the first guess values ( ' ). Figure 7 shows the effect of the filtering on the observation-space data using a 30-day window and a length-scale of 2°. Firstly, we can see that the large-scale features in the departures such as the overall bias and regional variations are removed. Secondly, the departures, observations and first guess distributions are put towards gaussianity, centred around zero and displaying more a symmetrical shape and tails. This then makes the processing and the 205 classification of the data much easier (see section 3.2).
To decide on the appropriate window length and Gaussian kernel length scale we have conducted sensitivity tests with different length scales (s= [0.5,1.0,2.0,5.0] degrees) and a window length of 10, 30 and 90 days. Figure 8 shows the resulting 12 possible sensitivity tests. If the Gaussian kernel is extended, i.e. 5 degrees, large-scale structures tend to be kept in the signal. If the kernel is narrow, i.e. 0.5°, only noise seems to be kept with a few very local minima and maxima. In the 210 latter case too much of the signal contained in the departures seems to be lost. If the time window is short, e.g. 10 days, coverage is low and isolated data points that are spotting possible outliers might be filtered out towards 0. Conversely, if the time window is long, i.e. 90 days, the sharp spatial structures that correspond to more recent or sporadic emission events are smoothed in the time averaging effect and not detected as outliers. For those reasons, we found that a time window of 30 days and a Gaussian kernel of 2.0° provides the most reasonable results. 215

Outlier classification
The final step is an outlier detection of the filtered departures. We choose to retain the values beyond the ±3σ range (ruling out 99.7% of the values). We found that a narrower range retains too many values and starts to fail isolating the important anomalies and conversely a wider range might fail to capture useful information. Further refinements to the current methodology could be done to find an optimal range for outlier detection using more advanced statistical methodologies. In 220 the present study we found that ±3σ provides suitable results. In addition to the outlier detection we perform a classification https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
given the relative values and sign of the filtered observations and first guess values. This allows, us to define the following four categories: • high observations (red in Fig. 9): where positive filtered observations are higher than filtered first-guess. This class is representative of high XCH4 values detected by TROPOMI that are not seen as high or at all in the forecasts. These are 225 likely originating from emissions that are not reported or under-estimated in the inventories. However, high observation categorisation may also be caused by poor quality observations category (see section 4.3).
• high forecasts (green in Fig. 9): where positive filtered first-guess are higher than filtered observations. This class is representative of high CH4 values in the forecasts but not seen as strong or at all in the TROPOMI XCH4 retrievals. High forecasts categorized data points are likely originating from emissions that are over-estimated or no longer being produced 230 or even mis-located in the emission inventory.
• low observations (blue in Fig. 9): where negative filtered observations are lower than filtered first-guess. This class is representative of locally low XCH4 values detected by TROPOMI but are not seen to be as low or at all in the forecasts.
Poor-quality observations influenced by low surface albedo likely fall in that category (see section 4.3).
• low forecasts (gold in Fig. 9): where positive filtered first guess are lower than filtered observations. This class is 235 representative of low XCH4 values in the forecasts but not seen as low or at all in the TROPOMI XCH4 retrievals. This category has generally much fewer data points. Orography could be a reason for data points to fall in that category, i.e. model surface height value that are higher than the observation value. Further developments of the method will likely use orography to improve the filtering.
In the maps, shades of the colours indicate the intensity of the offset, i.e. how far from perfectly matching observation 240 versus forecasts the filtered departure is. The size of the points indicates the number of samples. A larger dot indicates more data points within the 30-day window to compute the statistics hence is more robust. In the next section we will focus on the under-reported or missed plumes (red) category and over-reported or under-reported plumes (green) category to showcase the usefulness of the method.

Under-estimation of local sources in the forecasts
South Western USA and Mexico: In figure 9, the method detects under-predicted local CH4 concentrations (in red) in the forecast system in three areas. This occurs in the Permian Basin region, located around the Texas-New Mexico border, where multiple oil drilling sites are currently operating. Those enhancements have been documented by de Gouw et al. (2020) and Zhang et al, (2020) showing the reliability of the presented method. Two other regions with a smaller bias and extent can 250 be identified around the southern tip of Nevada and northern Baja California close to the US-Mexican border. To our https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. knowledge those two cases have not been investigated or documented yet. The latter case could be due to local albedo properties that could create local biases in the retrievals (see section 4.3), as we did not identify a facility responsible for those enhancements. This needs further investigation.
Western Turkmenistan: To confirm the ability of this methodology for the detection of large point-source emitters we 255 also showcase very strong detection of anomalous concentrations over the western Turkmenistan. Our system detects strong enhancement during all of 2019 (Fig.10). The filtered departures can be very large (above 50 ppb) with a high number count in the bins (large size of the dots). As mentioned earlier in this paper, anomalously large CH4 sources from oil and gas production in this location have been documented and detected using TROPOMI combined with private sector satellite data by Varon et al. (2019). 260

Over estimation of local sources in the forecasts
Western Russia: Our detection system shows two local point sources, that show large first-guess values that are not seen by TROPOMI XCH4 (green dots in Fig. 11). The features do not show large sampling (small dots) in time but do exhibit the shape of plumes, with strong departures near the point sources. One is very close to Moscow and corresponds to the Domodedovo airport surroundings. The other source detection is near the Volga river with a location matching small drilling 265 fields seen in visible satellite images. In these two locations the detection method suggests that emission inventories are overestimating local sources, which in reality are now producing reduced emissions or are no longer active emitters (at least during the period of monitoring).
Los Angeles: Similar features can occur in the area of Los Angeles. Figure 9 shows significant over-prediction of CH4 (green dots) specifically over San Bernardino and Palmdale. Both towns have industrial facilities and regional airports. Strong 270 high forecasted values are also detected in those locations on other dates (not shown) without being seen in the TROPOMI XCH4 measurements.
Such cases in very different locations show the capability of the method to detect not only missing or underreported point sources but also overreported cases. This can only be achieved with combining numerical models forecasts and satellite measurements at close-matching high horizontal resolution (9km and 7km respectively). It is also important to mention that 275 the method presented here is subject to uncertainties due to both model transport errors and representation error, although the error associated to emission generally dominate.

Local retrieval issues
The retrieval can be affected by albedo surface issues (see section 2). The filtering is not able to remove features with geographical extent smaller than the size of the Gaussian kernel (see section 3.1). Figure 12 presents a typical example of such 280 an issue where the same pattern is seen repeatedly in the outlier detection. This persistent pattern in shape and intensity corresponds to a land surface feature that produces consistently higher TROPOMI XCH4 values than its surroundings. The inverse can be true as well displaying local areas with consistently lower local XCH4 retrievals (not shown). Thus, great care https://doi.org/10.5194/acp-2020-550 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
should be taken when diagnosing such filtered departures. Features with a consistent shape and intensity are retrieval error artefacts, as atmospheric plumes would show more variability and not a consistent shape over months as illustrated in Fig. 12. 285 Further improvements of the method could be implemented with pattern shape recognition to discard persistent shapes over time. Satellite retrieval providers could be notified about of such biases, in order to improve the quality of the satellite product.

Conclusions
In this paper we have shown the potential of systematic detection of anthropogenic CH4 point and local source emissions relative to known emission inventory data using the TROPOMI satellite measurements in combination with high resolution 290 CH4 forecasts. While many studies have shown detailed analysis of a few case studies using TROPOMI observations, this is the first time that a systematic way to detect strong anthropogenic local emitters of CH4 and to compare results with emission inventories is presented. The method presented here does not only allow for the detection of unreported or missing sources but also targets over-reported sources in the inventories. The method also has the potential of detecting systematic local retrieval errors which can help to improve the satellite product. 295 Our method is novel by combining information from multiple sources (emission inventory, modelled surface fluxes, and observations) in a data assimilation framework to detect and analyse observed anomalies. We have used global emission inventories and fluxes that were the best possible global estimates we had available at the time when running our system. Using different emission inventories from research specific activities that are more specific to local regions, for instance, could provide different answers. In that way our methodology could provide an efficient way to validate improvements in sector-300 specific emission inventories. For example, using revised CH4 inventories such as presented by Maasakers et al. (2016) over the USA or more recently by Scrapelli et al., 2020 globally could lead to different detection patterns. Bottom up inventories will always lag in time and therefore cannot track rapid emission changes such as pipeline and gas facility blowouts. Satellite measurements have a clear added value for timely detection in the case of large emissions.
Combining satellite measurements, forecasts and emission inventories partially using a data assimilation system paves the 305 way to estimate the emissions themselves. Inverse modelling studies to estimate CH4 emissions have been done with SCIAMACHY and GOSAT CH4 satellite data generally performed at rather low resolution and focus specific study sites (e.g. Jacob et al., 2016). To our knowledge no published studies showed global inversions using TROPOMI data updating emissions close to the 10 km scale globally. Inverse modelling is computationally expensive and in the case of running operations beyond 10km scales to close-match satellite observations is a challenge that needs to be overcome over the next decade. Efforts are 310 underway to implement a sector-specific inverse high-resolution modelling monitoring system as part of the CAMS service evolution at ECMWF and the future Copernicus CO2 service at global and regional scales (e.g. Barré et al., 2019, Bousserez et al., 2019, Pinty et al., 2019, Janssens-Maenhout et al., 2020. Approaches combining global and regional modelling could be adopted to perform inversion at fine scales but at the cost of missing fine-scale detection outside the regional domains. Large and local CH4 emissions events could to occur in very remote areas, which are typically not considered in regional 315 modelling setups (e.g. West Turkmenistan, Varon et al.,2019). Systematic detection will then require setting up many regional subdomains leading again to computational burden for a single monitoring entity. We have demonstrated that monitoring of satellite XCH4 departures at high resolution at global scale using already existing parts of a forecasting chain remains an affordable solution to develop a much needed capability: tracking rapidly changing CH4 sources across the world and support the urgently needed effort on developing climate policies for reducing anthropogenic CH4 emissions.