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 2Regional 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 3Flow chart of the CAMS greenhouse gas analysis and forecast system.

Download

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 4Examples 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 5Examples 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

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

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

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 6Departure 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.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

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:

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 7Example 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 8High-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 9Examples 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 10Outlier detection and classification over the southwestern USA. Dates indicate the end date of the 30 d time window.

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

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 12Outlier detection and classification over western Russia. Dates indicate the end date of the 30 d time window.

Figure 13Outlier 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 15Outlier detection and classification over Siberia. Dates indicate the end date of the 30 d time window.

Figure 16Albedo 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 17Outlier detection and classification over Australia. Dates indicate the end date of the 30 d time window.

Figure 18Aerosol 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.