Systematic detection of local CH₄ anomalies by combining satellite measurements with high-resolution forecasts

In this study, we present a novel monitoring methodology that combines satellite retrievals and forecasts to detect local CH₄ concentration anomalies worldwide. These anomalies are caused by rapidly changing anthropogenic emissions that significantly contribute to the CH₄ atmospheric budget and by biases in the satellite retrieval data. The method uses high-resolution (7 km × 7 km) retrievals of total column CH₄ from the TROPOspheric Monitoring Instrument (TROPOMI) on board the Sentinel 5 Precursor satellite. Observations are combined with high-resolution CH₄ forecasts (∼ 9 km) produced by the Copernicus Atmosphere Monitoring Service (CAMS) to provide departures (observations minus forecasts) at close to the satellite's native resolution at appropriate time. Investigating these 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 synoptic-scale and meso-alpha-scale biases in both forecasts and satellite observations. We then apply a simple classification scheme to the filtered departures to detect anomalies and plumes that are missing (e.g. pipeline or facility leaks), underreported or overreported (e.g. depleted drilling fields) in the CAMS emissions. The classification method also shows some limitations to detect emission anomalies only due to local satellite retrieval biases linked to albedo and scattering issues.

1 Introduction

Atmospheric methane (CH₄) is the second most important anthropogenic greenhouse gas after carbon dioxide; it contributes significantly to changes in radiative forcing and climate change. CH₄ 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 CO₂ per unit mass (Hartmann et al., 2013). There are numerous natural and anthropogenic CH₄ sources, which vary in location and areal extent. The related anthropogenic emissions, such as oil and gas production, coal mining and biomass burning, tend to be geographically localised, e.g. over a plant facility, a pipeline or a field of extraction. However, methane emissions relating to biological fluxes such as livestock, landfills and rice fields can be either geographically localised over narrow areas or more widespread; for example, there are extensive patterns of microbial respiration in wetlands across the globe (Saunois et al., 2016). Atmospheric methane concentrations have more than doubled since preindustrial times because of an imbalance between methane sources and sinks (Hartmann et al., 2013) caused by increases in oil and gas production, rice crop cultivation, livestock farming and landfills. Methane has a relatively short atmospheric lifetime (with respect to climate scales) of around 9 years, meaning that targeted emission reductions could be an effective way to limit the rate of warming over the coming decades (Shoemaker et al., 2013).

Greenhouse gas emission inventories are generated through the aggregation and extrapolation of regional and national data. These data are reported individually by countries using the guidelines provided by the United Nations Framework on Climate Change (UNFCC) and the Intergovernmental Panel for Climate Change (IPCC). The reporting follows a bottom-up approach that utilises activity data and emission factors from individual emissions sectors. The process of officially reporting and processing these data to build these bottom-up inventories can cause a significant lag, so information can be out of date for certain sectors once publicly released. This can become an issue in the context of rapidly changing emissions from large point sources, for example in the oil and gas sectors (Alvarez et al., 2018). In the case of atmospheric composition modelling, emissions inventories are used as input surface fluxes to simulate atmospheric concentrations. Within the Copernicus Atmosphere Monitoring Service (CAMS), these simulations are used to provide routine real-time forecasts of greenhouse gas concentrations. The CAMS greenhouse gas forecasting system integrates satellite observations (Massart et al., 2014, 2016) to generate initial conditions for high-resolution forecasts at about 10 km (Agustí-Panareda et al., 2019). The lack of up-to-date emissions inventories impacts and likely degrades simulated CH₄ concentrations in areas where the local contribution from anthropogenic emissions is significant.

Many studies have demonstrated the rapidly changing and event-based nature of CH₄ anthropogenic emissions, especially when identifying the locations of “super-emitter” point sources. Conley et al. (2016) used aircraft measurements to characterise a blowout of a well connected to the Aliso Canyon gas storage facility in California from October 2015 to February 2016. Pandey et al. (2019) showcased the detection of high methane emissions 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 CH₄ source associated with a gas compression station using a combination of satellite instruments over Central Asia (western Turkmenistan). These types of suddenly occurring CH₄ emissions cannot be or are not reported/detected in time to be included in the bottom-up inventories, but they are seen from space. Other studies have shown the ability of satellite measurements to detect CH₄ emissions from extensive drilling and fracking areas. Kort et al. (2014) identified a large methane anomaly over the Four Corners region of the USA; more recently, de Gouw et al. (2020), Zhang et al. (2020) and Schneising et al. (2020) reported the satellite detection of large and extended CH₄ enhancements from various US oil- and gas-producing regions, such as the Permian Basin. While these satellite-based studies focused on specific events and locations, none of them systematically detected such anomalies at a global scale; nor did they provide a method to do so.

The systematic detection of large point sources of anthropogenic CH₄ emissions using a combination of satellite observations and modelling could prompt rapid action to reduce emissions from the oil and gas sectors. Two recent developments allow for the systematic detection of unreported CH₄ atmospheric anomalies linked to small-scale and point-source emissions: the newly available high-resolution (7 km × 7 km) satellite observations from the TROPOspheric Monitoring Instrument (TROPOMI; Veefkind et al., 2012) on board the Sentinel-5P platform and the improved real-time forecasting at high resolution (∼ 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 systematically detect atmospheric CH₄ anomalies relating to emission changes from small-scale and point sources that are not reported or are not updated in a timely manner. The paper is organised as follows: Sect. 2 describes the setup, including the TROPOMI observations and the forecasting and monitoring configurations; Sect. 3 presents the detection method; and Sect. 4 discusses several case studies that showcase the capabilities but also the limitations of the detection method. This is followed by our conclusions, where we briefly discuss the benefit that our approach can bring combined with coarse-resolution inverse modelling.

Figure 1 Global average TROPOMI XCH₄ column-averaged dry-air mixing ratios for the full year 2019, July 2019 and 1 July 2019 (a–c respectively).

2 Setup

2.1 TROPOMI CH₄ observations

The TROPOMI (Veefkind et al., 2012) instrument was launched on 13 October 2017 on board the Sentinel-5 Precursor satellite, a low Earth orbiter with a Sun-synchronous orbit that overpasses at 13:30 local solar time. 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, sulfur dioxide, formaldehyde, aerosol and methane (Hu et al., 2018). Current CH₄ observations, which are available for the inner two-thirds of the swath and only over land, are vertically integrated columns sensitive to the troposphere (surface to 200 hPa). With a swath that is around 1750 km (normally 2600 km) wide from the along-track position and a ground pixel size of 7 km × 7 km, TROPOMI CH₄ data can provide near-global daily coverage at high horizontal resolution over land, but they are limited by cloud cover and retrieval quality. In this study, we use the bias-corrected version of the product and 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-Readme-File, last access: January 2021). That document states that an overall bias of −0.3 % was found in a comparison against independent data, which is well within the mission requirements of ≤ 1.5 % (24 ppb). The scatter of the data around this bias also complies with the mission requirements of ≤ 1.0 % (18 ppb). Figure 1 illustrates the CH₄ satellite observation coverage that TROPOMI provides per year, per month and per day.

Figure 2 Regional zooms of TROPOMI XCH₄ columns for the full year 2019 in the Middle East (a) and central/western North America (b).

The measurements show clear geographical variation in the CH₄ column-averaged dry-air mixing ratios (XCH₄), which is driven to some extent by atmospheric transport but mostly by the spatial and temporal variability of surface fluxes and emission variations. Figure 2 shows the 2019 annual averages for the Middle East and the western USA. In these regions, spatial variability results in XCH₄ enhancements of up to 50 ppb over emission hotspots. The average concentrations in these regions – approximately 1825 ppb over the USA and approximately 1875 ppb over the Middle East – are also significantly different. The strong local enhancements are an indication of strong local surface fluxes and emissions of CH₄ from oil and gas activities, mining, agriculture or wetlands. XCH₄ retrievals can also be prone to some systematic residual errors, especially those relating 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 shortwave-infrared spectral band during winter. Even though they are generally reduced in the bias-corrected product, such retrieval biases need further investigation (see Sect. 4.3). Nevertheless, the TROPOMI data are sufficiently accurate to show local enhancements linked (but not limited) to oil and gas production. We show in Sect. 3 how to isolate these small-scale signals of interest and remove the contribution of synoptic-scale (more than 2000 km) and meso-alpha-scale (between 2000 and 200 km) biases. In the rest of the paper, we define “large scale” as the combination of synoptic scale and meso-alpha scale.

2.2 CAMS high-resolution CH₄ 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. As part of the CAMS greenhouse gas services, the IFS provides 5 d forecasts for CO₂ and CH₄ (Agustí-Panareda et al., 2019) along with other species that are relevant to air quality (Flemming et al., 2015).

Figure 3 Flow chart of the CAMS greenhouse gas analysis and forecast system.

The IFS model cycle used in this paper is CY45R1. This is run routinely at TCo1279 horizontal resolution, which is a cubic octahedral reduced Gaussian grid at approximately 9 km resolution (Holm et al., 2016), 137 vertical levels from the surface to 0.01 hPa, and a time step of 450 s. Details about the transport and meteorological configuration can be found in Agustí-Panareda et al. (2019). The CAMS greenhouse gas (GHG) operational suite is composed of an analysis and forecasts at medium and high resolution (see Fig. 3). The analysis is based on the IFS 4D-Var assimilation system, which was adapted to assimilate retrieved column-averaged mole fractions of CO₂ and CH₄ together with all the operational meteorological observations (Engelen et al., 2009: Massart et al., 2014, 2016). The analyses are produced every 12 h (00:00 and 12:00 UTC). A 4 d forecast is then issued daily after the 00:00 UTC analysis on a TCo399 cubic octahedral grid corresponding to approximately 25 km × 25 km with the same 137-model-level configuration. Two satellite observation streams are currently assimilated: the Infrared Atmospheric Sounding Interferometer (IASI) for CH₄ on the MetOp satellites and the Thermal And Near-infrared Sensor for carbon Observations (TANSO) on the GOSAT satellite for both CO₂ and CH₄ (see Massart et al., 2014, for further details). In this configuration, only the concentrations are corrected by the assimilation; the emissions and surface fluxes remain unchanged. Processing and acquisition of the level 2 data in 2019 yielded the satellite XCH₄ data 4 d behind real time. The high-resolution forecast is then coupled to the analysis experiment by merging the 4 d lower-resolution forecast from the CO₂ and CH₄ analysis with the previous 1 d high-resolution forecast (Fig. 3) in order to preserve the fine-scale features of the high-resolution forecast. Additionally, the high-resolution forecast 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 CH₄ forecasts at high resolution from the setup described above. The high-resolution forecasts are run on a TCo1279 L137 grid of approximately 9 km × 9 km for a 5 d period and are initialised approximately 4 h behind real time every day from 00:00 UTC. The impact of the TANSO and IASI CH₄ retrieval assimilation is not strong close to the surface (see Fig. 8 in Massart et al., 2014). Analyses performed at a resolution of 25 km do not correct the emissions and hence mainly provide a correction to the forecast initial condition concentrations in the free troposphere and above. At lower altitudes, the emissions are the dominant influence on the 4 d forecast at 25 km, which is used to initialise the high-resolution forecast at 9 km that does not include data assimilation.

Figure 4 Examples of combined net fluxes (only positive values are shown due to the logarithmic scale) that constitute the surface boundary conditions of the IFS high-resolution CH₄ forecast. Global- and regional-scale examples for 15 January and 15 July 2019 at 12:00 UTC are presented.

Figure 5 Examples of the IFS high-resolution CH₄ forecast output, showing snapshots at global and regional scales of the total column mean molar fractions for 15 January and 15 July 2019 at 12:00 UTC. The lower panels show parts of Europe (left) and the Middle East (right).

Both high-resolution forecasts and analysis use prescribed CH₄ surface fluxes. The anthropogenic emissions data, including fossil fuel, agricultural and landfill/waste emissions, are from the EDGARv4.2FT2010 dataset (Olivier and Janssens-Maenhout, 2012) for 2010 with 0.1^∘ × 0.1^∘ and monthly resolution for the rice emissions (Matthews et al., 1991). Monthly mean wetland emissions come from a climatology (1990–2008) based on the LPJ-WHyMe model that is constrained by SCIAMACHY observations during the HYMN project (Spahni et al., 2011) and has a resolution of 1^∘ × 1^∘. The biomass burning emissions are from GFASv1.2 (Kaiser et al., 2012). Other sources and sinks include a monthly soil sink (Ridgwell et al., 1999), annual mean oceanic fluxes (Houweling et al., 1999; Lambert and Schmidt, 1993), and monthly mean fluxes from termites (Sanderson, 1996) and wild animals (Houweling et al., 1999). The chemical sink in the troposphere and stratosphere is represented by a climatological monthly mean chemical loss rate (Bergamaschi et al., 2009). This is based on OH fields optimised with methyl chloroform using the TM5 model (Krol et al., 2005) with prescribed concentrations of stratospheric radicals obtained using the 2D photochemical Max Planck Institute model. Figure 4 shows the geographical and seasonal structure of the surface fluxes. Large-scale and smoother structures are representative of wetland, soil and agricultural fluxes, whereas the finer-scale and sharper structures are representative of anthropogenic and fire emissions. Figure 5 shows the usefulness of the high-resolution forecasts at global and regional scales. Global seasonal cycles and synoptic-scale concentrations are represented, as well as concentrations at smaller scales, such as plumes from point-source emissions and orographic effects. Large point sources and associated plumes can be seen over Europe, for example over Madrid, Paris and Tours (western France). Inventory estimates suggest that the modelled hotspot region near Tours is probably the result of solid-waste landfill emissions. Other possibilities include emissions from the enteric fermentation and wastewater treatment sectors, all of which may be linked to a landfill site. Over the Middle East region zoom, sharp point sources are seen in Teheran and southern Iran, as well as over Pakistan (Karachi) and closer to the Himalayan region.

2.3 Monitoring suite

To monitor and compare the TROPOMI XCH₄ retrievals with the IFS CH₄ 9 km forecasts, we reuse part of the IFS assimilation system in a so-called monitoring mode. The system recomputes a high-resolution trajectory at 9 km (which is a model integration), initialised from the forecasts in a 12 h monitoring window, to calculate so-called first-guess departures (a departure is the difference between an observation and the corresponding model forecast) using the observations at the appropriate model time step. At each observation location, the departure can be written

$$\begin{array}{}\text{(1)}& d=y-HM\left(x_i\right),\end{array}$$

where d is the departure, y is the observation, H is the observation operator, M is the model integration or trajectory, and x_i is the initial CH₄ condition at the beginning of the monitoring window. If we inject the retrieval equation (Rodgers, 2000), the departure becomes

$$\begin{array}{}\text{(2)}& d=\mathbf{A}{x}_{\mathrm{t}}+(\mathbf{I}-\mathbf{A})x_{\mathrm{a}}+\mathit{ϵ}-\mathbf{A}M\left(x_i\right)-(\mathbf{I}-\mathbf{A})x_{\mathrm{a}},\end{array}$$

where x_t is the true CH₄ concentration state (which is never known exactly), A is the averaging kernel matrix (which represents the sensitivity of the retrieval on the vertical profile to the true state), I is the identity matrix, x_a is the information known a priori that is used in the retrieval, and ϵ is the retrieval error term. The equation then simplifies to

$$\begin{array}{}\text{(3)}& d=\mathbf{A}(x_{\mathrm{t}}-M(x_i\left)\right)+\mathit{ϵ},\end{array}$$

which is the difference between the true state and the forecast smoothed by the averaging kernel function plus the retrieval error term. Those departure values are thus strongly dependent on the averaging kernel function shape. For the TROPOMI XCH₄ retrievals, the mean averaging kernel function shows a homogeneous sensitivity to the entire troposphere, slightly decreasing in the stratosphere (see Fig. 2 in Hu et al., 2016). The averaging kernel function does not vary markedly between pixels or regions of the globe (not shown). Figure 6 shows the departures over various timescales (yearly, monthly and daily) for the global domain. Overall, the departures (observation minus forecast) show a global positive bias of around 25 ppb (meaning that the observation values are higher than the model values), which can be attributed to model biases (Ramonet et al., 2019) and/or observation biases (Langerock et al., 2019). Ramonet et al. (2019) compared the CAMS CH₄ forecasts with independent total column data. Results showed that the forecasts continuously underestimated the CH₄ total columns by 5–20 ppb. Langerock et al. (2019) showed that the average total column bias for the TROPOMI CH₄ retrievals was −0.32 % (i.e. around −5 ppb) with respect to ground-based measurements.

Figure 6 Departure values computed with the observations displayed in Fig. 1 for the full year 2019, July 2019 and 1 July 2019.

Regional-scale error structures are evident from the observation–model comparison. For example, boreal regions show a band of negative values that are potentially attributable to systematic errors caused by surface albedo values during winter (see Sect. 2.1) in the TROPOMI retrieval algorithm. Alternatively, they could be caused by CH₄ biases at the tropopause and lower stratosphere levels in the IFS model. Also, a possible time lag in the wetland emissions, which are calculated offline and provide boundary conditions in the IFS forecasting chain (see Sect. 2.2), could cause such errors. The origin of this type of large-scale error in the departures is not yet fully understood and is beyond the scope of this paper, although an understanding of these biases is crucial to further improving the quality of CAMS CH₄ forecasts and TROPOMI retrievals. At finer scales, structures are seen on the yearly average comparison that become more evident at the monthly timescale. Local differences are even stronger at the daily timescale, but recognising fine-scale structures is challenging due to the lack of daily coverage. For these reasons, spatial filtering and temporal averaging of the departures are performed to extract and use the small-scale features seen in the departures.

3 Detection method

3.1 Filtering the signal

To remove the large-scale features seen in the departures, we have implemented 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 0.1^∘ latitude–longitude binning. Due to the ocean, cloud cover and quality control flagging, a number of departure bins will have missing values, which will jeopardise the convolution. This problem is solved technically by creating two auxiliary matrices in which missing values are replaced by 0. The two auxiliary matrices are then defined as

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathbf{D}=\left\{\begin{array}{l}{d}_{\mathrm{m}}\text{if}n>N\\ \mathrm{0}\text{otherwise}\end{array}\right.\text{(5)}& {\displaystyle }& {\displaystyle }\mathbf{C}=\left\{\begin{array}{l}\mathrm{1}\text{if}n>N\\ \mathrm{0}\text{otherwise},\end{array}\right.\end{array}$$

where d_m (the subscript m refers to “mean”) is the average departure in the given bin and N is the minimum number of observations that must be included in a given bin. In this study, we use N=2 to avoid smoothing with very isolated (and potentially faulty) pixels while retaining as much data as possible. Replacing the missing values in D with zeros introduces an error in the filtered departures d_hp (where the subscript hp refers to “high pass”) after convolution (lower values are induced by the smoothing resulting from the insertion of zeros). This can be corrected for by applying the same Gaussian filter to a matrix C representing the bins selected for filtering (where the number of counts is at least N) and using the ratio of the two filtered matrices to compensate for the missing-value errors. A high-pass filter on a given observation space field (departures here) can then be formulated as follows:

$$\begin{array}{}\text{(6)}& d_{\text{hp}}=d_{\mathrm{m}}-{\displaystyle \frac{\mathbf{G}\left(\mathit{\sigma }\right)*\mathbf{D}}{\mathbf{G}\left(\mathit{\sigma }\right)*\mathbf{C}}},\end{array}$$

where G(σ) is a 2D Gaussian kernel function with a σ length scale. The same filtering is also applied to the observation values y and the first-guess values HM(x_b), as these will be used for classification in Sect. 3.2. Figure 7 shows the effect of the filtering on the observation-space data using a 30 d 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 departure, observation and first-guess distributions are made more Gaussian, such that the distributions is centred at zero and is more symmetrical, and display tails. This makes the processing and classification of the data much easier (see Sect. 3.2).

Figure 7 Example of the effect of high-pass filtering on observation-space data using a 30 d window ending on 1 July 2019 with a 2^∘ Gaussian kernel length scale: (a) the unfiltered departures, (b) the filtered departures, (c) histograms comparing unfiltered (red) vs. filtered (blue) departures, and (d) 2D distributions of the unfiltered (red) and filtered (blue) data points in the observation–forecast space. Note that the unfiltered data distribution is centred at the mean in this plot to facilitate comparisons with the filtered distribution.

Figure 8 High-pass filtering sensitivity tests of the departures normalised to the instrument precision. Window length (columns) is 10, 30 or 90 d, and kernel size (rows) is 0.5, 1.0, 2.0 or 5.0^∘. The selected filter parameter values are those for the plot outlined in bold.

To decide on the appropriate window length and Gaussian kernel length scale, we conducted sensitivity tests with different length scales (σ = [0.5, 1.0, 2.0, 5.0]^∘) and window lengths of 10, 30 and 90 d. Figure 8 shows the resulting filtered departures normalised to the instrument precision for the 12 possible sensitivity tests. For the tests with Gaussian kernel sizes of 0.5 and 1.0^∘, the filtered signal is generally weaker than the measurement precision, so very few or no detections of local anomalies will be made. Conversely, if the kernel is large, the signal is stronger relative to the instrument precision, but there is a risk of picking up larger patterns than the targeted features, i.e. features that are directly related to local emissions in the atmospheric distribution of CH₄. For these reasons, we found that a kernel of 2.0^∘ performed best. If the time window is short, e.g. 10 d, decreased coverage could hinder the correct detection of outliers, especially for isolated data points. Isolated data points that indicate possible outliers could be filtered out towards 0, as the convolution does not have neighbouring points to use within the kernel range. The shorter the window and the narrower the kernel, the more likely this is. Conversely, using a long time window, i.e. 90 d, maximises the chances of achieving good observation coverage, allowing the most effective use of the convolution filter, although this reduces the ability to provide information on temporal variability. Also, the sharp spatial structures that correspond to more recent or sporadic emission events are smoothed by the time averaging, decreasing the filtered departure over instrument precision ratio. For these reasons, we found that a time window of 30 d provides the most reasonable results.

Figure 9 Examples of outlier classification. (a, b) Global distributions in the observation–first-guess space for two different 30 d windows (end dates: 1 July and 1 September 2019). The four data classes are shown in different colours; the colour intensity reflects the number of outliers. (c, d) Locations of outlier classes around the globe during the two 30 d windows (end dates: 1 July and 1 September 2019, respectively). Darker dots show larger departures. Larger dots indicate that more occurrences were detected in the bin and time window.

3.2 Outlier classification

The final step is outlier detection of the filtered departures. We choose to retain the filtered departures with absolute values that are larger than the CH₄ measurement precision of TROPOMI. If the absolute value of a filtered departure is lower than the measurement precision, it is considered noise and ruled out. The current methodology could be refined to find a better outlier detection method based on more advanced statistical techniques. In the present study, we found that the measurement precision of the satellite product provides suitable results. In addition to outlier detection, we classify the departures into four categories based on the relative values and signs of the filtered observations and forecast values:

• High observations (red in Fig. 9), where the filtered observation values are higher than the absolute filtered forecast values. This class represents XCH₄ values that are high according to TROPOMI but are not seen as high or not seen at all in the forecasts. These likely originate from emissions that are not reported or are underestimated in the inventories. However, high observations may also be caused by poor-quality observations due to albedo and scattering issues (see Sect. 4.3).

• High forecasts (green in Fig. 9), where the filtered forecast values are higher than the absolute filtered observation values. This class represents CH₄ values that are high in the forecasts but are not seen as strong or not seen at all in the TROPOMI XCH₄ retrievals. High forecasts likely originate from emissions that are overestimated, no longer being produced, or even mislocated in the emissions inventory.

• Low observations (blue in Fig. 9), where the filtered observation values are lower than the absolute filtered forecast values. This class represents XCH₄ values that are locally low according to TROPOMI but are not seen to be as low or are not seen at all in the forecasts. Poor-quality observations that are influenced by low surface albedo likely fall in this category (see Sect. 4.3).

• Low forecasts (gold in Fig. 9), where the filtered forecast values are lower than the absolute filtered observation values. This class represents XCH₄ values that are low in the forecasts but are not seen to be as low or are not seen at all in the TROPOMI XCH₄ retrievals. This category generally includes relatively few data points, which are very sparsely distributed. Orography may explain the data points in this category (i.e. model surface heights that are higher than the corresponding observation values). Further improvements to our method will likely involve the use of orography to improve the filtering.

In the maps shown in Fig. 9, colour intensity indicates the magnitude of the offset, i.e. how far the observation value is from the forecast value. Point size indicates the number of samples included: a larger dot indicates that more data points within the 30 d window were used to compute the statistics, making them more robust. Figure 9 gives an overview of outliers detected around the world, which are many and varied. In the next section, we focus on specific case studies from the underreported/missed plumes (red) category and the overreported plumes (green) category to demonstrate both the usefulness of the method and its current limitations.

Figure 10 Outlier detection and classification over the southwestern USA. Dates indicate the end date of the 30 d time window.

Figure 11 Outlier detection and classification over Turkmenistan. Dates indicate the end date of the 30 d time window.

4 Case studies

4.1 Underestimation of emissions from local sources in the forecasts

Southwestern USA and Mexico. Figure 10 shows that the method detects underpredicted local CH₄ concentrations (in red) in the forecast system for three areas of the southwestern USA and Mexico. One area is in the Permian Basin, around the Texas–New Mexico border, where multiple oil drilling sites are currently operating. These enhancements have been documented by de Gouw et al. (2020) and Zhang et al. (2020), showing the reliability of the presented method. The two other regions have smaller biases and extents, and they occur around the southern tip of Nevada (by Lake Mead) and northern Baja California (close to the USA–Mexico border). To our knowledge, these two cases are yet to be investigated or documented, so they need further investigation. We could not identify any facilities that are responsible for those enhancements; neither could we find any matching land surface albedo features that could create local biases in the retrievals (see Sect. 4.3) from visible satellite imagery. The northern Baja California enhancements are observed to be correlated with scattering parameters such as the aerosol optical thickness (AOT, not shown), which should be investigated further to improve the method. A case study of the influence of the AOT on the detection method is described in Sect. 4.3.

Western Turkmenistan. To confirm the ability of this methodology to detect large point-source emitters, we also showcase the very strong detection of anomalous concentrations over western Turkmenistan. Strong enhancement occurred throughout 2019 (Fig. 11), although the enhancement changed in intensity and shape. The filtered departures can be very large (above 50 ppb), with high counts in the bins (i.e. large dots). As mentioned above, anomalously large CH₄ sources linked to oil and gas production in this location have been documented and detected by Varon et al. (2019) using TROPOMI in combination with private-sector satellite data.

Figure 12 Outlier detection and classification over western Russia. Dates indicate the end date of the 30 d time window.

Figure 13 Outlier detection and classification over Libya, Egypt and Niger. Dates indicate the end date of the 30 d time window.

4.2 Overestimation of emissions from local sources in the forecasts

Western Russia. Our detection system shows two local point sources with large forecast values that are not seen by TROPOMI XCH₄ (green dots in Fig. 12). The features are depicted as small dots (i.e. relatively few samples) that form a plume shape, with strong departures near the point sources. One is very close to Moscow and corresponds to the surroundings of Domodedovo Airport. The other source is near the Volga River; its location correlates with small drilling fields seen in visible satellite images. The detection method suggests that emission inventories are overestimating emissions from these sources, which were actually releasing lower emissions or were no longer active emitters during the monitoring period.

Los Angeles. Similar features can occur in the area of Los Angeles. Figure 10 shows significant overprediction of CH₄ (green dots) over San Bernardino and Palmdale. Both towns have industrial facilities and regional airports. Stronger detection is seen in the 1 September and 1 October 2019 windows than in the 1 July and 1 August 2019 windows. Differences in intensity may be attributed to changes in emission levels from month to month, as well as seasonal atmospheric transport changes due to different meteorological situations in different windows. For example, if the overall wind speed increases near the source, less accumulation of CH₄ would be seen, leading to smaller departures and less detection.

Such cases in very different locations show the ability of the method to detect not only missing or underreported point sources but also overreported cases. This can only be achieved by combining numerical model forecasts and satellite measurements that have closely matched high horizontal resolutions (9 and 7 km respectively). It is also important to mention that the method presented here is subject to uncertainties due to model transport error and representation error, although the error associated with emission levels generally dominates. Further work is needed to account for atmospheric transport and, more generally, weather variability in the detection method. Techniques such as those described in Barré et al. (2020) show great potential for application in this context.

Figure 14 (a, b) TROPOMI albedo in the near-infrared (NIR) band, (c, d) TROPOMI albedo in the shortwave-infrared (SWIR) band, and (e, f) TROPOMI XCH₄ columns. Maps display averages for the same windows as shown in Fig. 13.

4.3 Local retrieval issues

The TROPOMI XCH₄ retrieval can be affected by albedo issues (see Sect. 2), and the filtering is not able to remove features with geographical extents that are smaller than the size of the Gaussian kernel (see Sect. 3.1). Figure 13 shows patterns in the outlier detections where the filtered observations are lower than the filtered forecasts (low observation category; plotted in blue). Similar patterns are seen for the two months considered. In Fig. 14, the TROPOMI albedo in the near-infrared (NIR) and shortwave-infrared (SWIR) bands is displayed along with the XCH₄ column retrieval. There are large albedo variations in the two spectral bands in this area that affect the XCH₄ columns. The patterns for the albedo and the XCH₄ columns clearly match. The ability of the filtering algorithm to remove this effect depends on the structure of the pattern and how narrow or small it is. This is an issue with the outlier detection method in its current form.

Figure 15 Outlier detection and classification over Siberia. Dates indicate the end date of the 30 d time window.

Figure 16 Albedo in the shortwave-infrared (SWIR) band over Siberia, as provided by TROPOMI. Maps display averages for the same windows as shown in Fig. 13.

Another example is provided in Fig. 15, where the same pattern is seen repeatedly in the outlier detection over the Siberian region. In this case, the filtered observation values are higher than the filtered forecast values (high observation category; plotted in red). The pattern detected matches the SWIR albedo pattern shown in Fig. 16. The high-albedo areas in the SWIR albedo map are also the areas with relatively high TROPOMI XCH₄ compared to their surroundings, and the pattern is narrow enough to be missed by the filter. Thus, great care should be taken when diagnosing such filtered departures. Features with a consistent distinctive shape and intensity are potential retrieval error artefacts, as atmospheric methane signals generally show considerable variability as a result of meteorological dynamics, not a consistently distinctive shape over a period of months (8 and 10 weeks in Figs. 13 and 15 respectively). Further development of our detection method should include the incorporation of albedo information to try to account for such systematic local biases.

Figure 17 Outlier detection and classification over Australia. Dates indicate the end date of the 30 d time window.

Figure 18 Aerosol optical thickness (AOT) in the near-infrared (NIR) band over Australia, as provided by TROPOMI. Maps display averages for the same windows as shown in Fig. 15.

However, in certain cases, identifying such biases can be more complex than just correcting outlier detections using albedo information. In Fig. 17, the anomaly presents a consistent shape for the four months shown when using 30 d windows, but there is no feature in the corresponding NIR and SWIR albedo maps that is clearly correlated with the anomaly. Instead, investigations showed that this anomaly is associated with scattering parameters such as the aerosol optical thickness (AOT); it correlates well with the AOT in the NIR and SWIR bands. As an example, Fig. 18 shows the presence of a feature in the map of the AOT in the NIR band that correlates to high-XCH₄ observations with respect to the forecasts. A similar feature is also seen in the map of the AOT in the SWIR band (not shown). This shows that more work is needed to identify systematic local biases, as they are not all caused by albedo variations. Additionally, there is still a potential impact of temporally and spatially variable small-scale albedo features. For example, nonpersistent biases can arise with snow coverage. More research will be needed to improve automatic local bias detection for less persistent biases.

Currently, the method presented here helps to identify retrieval biases, but it does not do this systematically. Further improvements to the method could include the joint use of albedo and scattering parameters to perform additional filtering correction and flagging. The availability of such an automatic retrieval bias detection tool would help to improve the quality of the retrieval product.

5 Conclusions

In this paper, we have shown the potential for the systematic detection of anthropogenic CH₄ point and local source emissions relative to known emissions inventory data using TROPOMI satellite measurements in combination with high-resolution CH₄ forecasts. While many studies have included detailed analyses of a few case studies using TROPOMI observations, this is the first step towards providing a systematic way to detect strong local anthropogenic emitters of CH₄ and to compare the results with emissions inventories. The method presented here not only shows the potential for the detection of unreported or missing sources, but it also targets overreported sources in the inventories. It also has the potential to identify systematic local retrieval errors, which could help to improve the satellite product. However, the current method has some limitations, as it requires additional correlative analyses with albedo and scattering parameters to account for local biases in the retrieval values. Complications can arise if land-surface-induced biases and detected local emission features are collocated. The presence of multiple types of collocated emissions (i.e. anthropogenic plus wetland emissions) will also increase the complexity, as emission patterns could show albedo correlations in this case. We have, however, demonstrated the potential of the methodology by focusing on several case studies, and further work is required to provide a global assessment based on several years from this dataset.

Our method is novel in that it combines information from multiple sources (emissions inventories, modelled surface fluxes and observations) in a data assimilation framework in order to detect and analyse observed anomalies. We used the global emissions inventories and fluxes that were the best possible global estimates available to us when we were running our system. Using other emissions inventories from research activities that are more specific to local regions could provide different answers. In this way, our methodology could provide an efficient way to validate improvements in sector-specific emission inventories. For example, using revised CH₄ inventories such as those presented by Maasakers et al. (2016) for the USA or more recently by Scrapelli et al. (2020) for the entire globe 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 clear added value given that they can provide timely detection of large emissions.

Combining satellite measurements, forecasts and emissions inventories using a data assimilation system paves the way to estimating the emissions themselves. While inverse modelling studies to estimate CH₄ emissions have been performed with SCIAMACHY and GOSAT CH₄ satellite data, they have generally been performed at rather low resolution and have focused on specific study sites (e.g. Jacob et al., 2016). To our knowledge, there are no published studies that report global inversions based on TROPOMI data and update emissions at close to the 10 km scale globally. Inverse modelling is computationally expensive, and modelling at scales of more than 10 km to closely match satellite observations is a challenge that must be overcome in the next decade. Efforts are 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 CO₂ 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 regional domains. Large and local CH₄ emission events could occur in very remote areas that are typically not considered in regional modelling setups (e.g. western Turkmenistan; Varon et al., 2019). Systematic detection will then require the establishment of many regional subdomains, again increasing the computational burden on a single monitoring entity. We have demonstrated that monitoring satellite XCH₄ departures at high resolution at the global scale using pre-existing parts of a forecasting chain is an affordable solution to a highly desirable aim: to track rapidly changing CH₄ sources around the world, thus supporting efforts to urgently develop climate policies aimed at reducing anthropogenic CH₄ emissions.