Joint Use of in-Scene Background Radiance Estimation and Optimal Estimation Methods for Quantifying Methane Emissions Using PRISMA Hyperspectral Satellite Data: Application to the Korpezhe Industrial Site

Methane (CH4) is one of the most contributing anthropogenic greenhouse gases (GHGs) in terms of global warming. Industry is one of the largest anthropogenic sources of methane, which are currently only roughly estimated. New satellite hyperspectral imagers, such as PRISMA, open up daily temporal monitoring of industrial methane sources at a spatial resolution of 30 m. Here, we developed the Characterization of Effluents Leakages in Industrial Environment (CELINE) code to inverse images of the Korpezhe industrial site. In this code, the in-Scene Background Radiance (ISBR) method was combined with a standard Optimal Estimation (OE) approach. The ISBR-OE method avoids the use of a complete and time-consuming radiative transfer model. The ISBR-OEM developed here overcomes the underestimation issues of the linear method (LM) used in the literature for high concentration plumes and controls a posteriori uncertainty. For the Korpezhe site, using the ISBR-OEM instead of the LM -retrieved CH4 concentration map led to a bias correction on CH4 mass from 4 to 16% depending on the source strength. The most important CH4 source has an estimated flow rate ranging from 0.36 ± 0.3 kg·s−1 to 4 ± 1.76 kg·s−1 on nine dates. These local and variable sources contribute to the CH4 budget and can better constrain climate change models.


Introduction
Climate change has become a major preoccupation for society, leading nations to develop common environmental commitments, such as the Paris Agreement (2015). Methane (CH 4 ) is one of the greenhouse gases (GHGs) chiefly under scrutiny due to its key role in several activity domains. Firstly, CH 4 has significant explosive properties that may result in major infrastructure damage and security issues [1]. CH 4 also has an impact on public health; inhaling high concentrations can lead to respiratory complications. CH 4 also contributes to the formation of tropospheric ozone (O 3 ), which may have a much higher impact in terms of respiratory and cardiovascular deaths [1][2][3]. Regarding the environment, CH 4 is the second most important anthropogenic GHG in terms of impact on climate and influence on air quality [4]. There is a high concentration of carbon dioxide (CO 2 ) in the atmosphere, yet CH 4 has a global warming potential approximately 28 times higher on a 100-year scale than that of CO 2 [5][6][7][8][9]. Calculations of CH 4 impact on climate change must consider indirect effects. The oxidation of CH 4 in the presence of nitrogen oxides and sunlight lead to the formation of ozone (O 3 ) [10,11]. O 3 adds approximately 38% to CH 4 force [2]. Considering all these impacts, the broad social cost of methane resulting from one ton of emission is 50-100 times greater than that of CO 2 [12]. Fortunately, CH 4 has an atmospheric lifetime of about a decade, while the CO 2 lifetime is approximately a century [10,13,14], meaning that reducing emissions will have a quicker impact for mitigating climate change [15][16][17].
Nevertheless, the CH 4 concentration in the atmosphere estimated by ice cores has changed dramatically, from approximately 650 parts per billion (ppb) in the year 1000 to approximately 1700 ppb shortly before the year 2000 [18]. The highest increase in the National Oceanic and Atmospheric Administration's 37-year record [2] was observed in 2020, with an atmospheric concentration of 1890 ppb. The tripling of CH 4 concentration in a thousand years, as well as the increase in anthropogenic CH 4 compared to natural CH 4 , can be explained by human activity. Natural CH 4 sources are wetlands (swamps), termites and oceans, which are relatively time-stable [19]. Anthropogenic sources are distributed in the sectors that followed the industrial revolution: energy (fossil fuel production, etc.), industry, waste treatment (15 to 18% of the anthropogenic balance [20]) and agriculture (rice farming, ruminant farming, etc.) [21]. CH 4 sinks are divided into three categories [22]; the main cause of methane destruction (90% of the sinks [23]) is the oxidation of CH 4 by the hydroxyl (OH) radical in the troposphere, while absorption by soils and oxidation in the stratosphere represent the rest of the CH 4 sinks.
Improved detection and quantification of CH 4 emissions would provide a better constraint for inventories of large-scale CH 4 releases into the atmosphere. Due to CH 4 s explosive properties and the limited accessibility of the source (in the case of major gas leaks), remote sensing is a powerful tool, already validated up to quantification [24]. Satellite sounders, such as TROPOMI (TROPOspheric Monitoring Instrument) [25] on ESA's Sentinel-5 Precursor satellite, are suitable instruments to study CH 4 emissions on a global scale [9]. However, their kilometric spatial resolution is unsuitable for mapping anthropogenic CH 4 from industry. The dimensions of an industrial plume are typically less than a few kilometers. Detecting industrial plumes is an economic and environmental challenge. Firstly, permanent CH 4 leakage represents a loss of revenue for industries [26]. Secondly, although large emitters dominate the industrial methane emissions budget [27], small sources also contribute to increased CH 4 emissions. Identifying a leak is always a good step towards reducing GHG emissions. Other instruments, called imagers, have a smaller ground pixel size than sounders. The most common are airborne imagers that provide a ground pixel size of approximately a few meters. Airborne imagers, such as AVIRIS-NG (Airborne Visible/Infrared Imaging Spectrometer-Next Generation) [28], can detect and quantify industrial plumes [29]. These imagers can detect small-scale leakages on industrial infrastructure. Once located, dedicated teams can visit the site or a campaign can focus on certain emission points [30]. However, these airborne campaigns are expensive, difficult to perform, and cover only a small portion of the area of interest.
Satellite imagers are a compromise between satellite sounders and airborne imagers. Satellite imagers developed rapidly since 2020; however, it is important to note that two of them are pioneers. Hyperion [31], launched in the late 2000s, detects and quantifies industrial plumes, but the results are established with a super-emitter source [29] and the signal-to-noise ratio (SNR of 20 [32]) of this instrument has limited further studies [33] on weaker sources. Launched in 2016, the GHGSat-D (GreenHouse Gases Satellite) satellite demonstrated the feasibility of detecting CH 4 plumes with a spatial resolution of 50 m, but the data are not freely accessible [34]. The new imaging satellite, PRISMA (PRecursore IperSpettralle della Missione Applicativa) [35], launched in 2019, has a spatial resolution of 30 m and hyperspectral images are available. Recent work has demonstrated the feasibility of using this satellite to quantify methane plumes [36][37][38].
In the literature, two concentration map methods are widely used: the Matched Filtering (Linear Method (LM)) [23,29,[39][40][41][42][43][44][45][46] and the Optimal Estimation (OE) method [37,[47][48][49][50]. The first method linearizes the mathematical expression of the gas transmission, resulting in an underestimation of the retrieved concentration. In the second method, the standard OE approach [28] requires significant computation time due to the non-linear radiative transfer model computation and the number of parameters to retrieve. In this paper, we propose a new algorithm called ISBR-OE that builds on the LM and OE approaches. This paper focuses on the Korpezhe oil and gas field, located in a desert area of Turkmenistan. This site was chosen because several significant emission points are present and the background is homogeneous [36]. Based on these data, this paper first highlights the drawbacks of the LM method widely used in the literature, resulting in an underestimation of the retrieved concentration. Then, the second part details a new plume quantification method, called the ISBR-OE method, that builds on the LM and OE approaches in order to obtain results with uncertainty control and reasonable computation time. In order to understand the improvement in the new method, we compared the integrated path concentration estimated from the linear and the new methods using synthetic data and real PRISMA data on different CH 4 plumes. Furthermore, two state-of-the-art flow rate estimation methods and a new method derived from the literature were applied to PRISMA data. We have analyzed CH 4 emission uncertainties derived from integrated path concentration maps for the industrial area studied. All three flow rate estimation methods are based on wind speed. Some approaches try to circumvent using this information by performing large-eddy simulations (LES) to derive the flow rate as a function of the plume shape [51], but these methods are not covered here.

Hyperspectral Images
Launched in March 2019, the PRISMA satellite is equipped with hyperspectral and panchromatic sensors. The Agenzia Spaziale Italiana's (ASI) policy is to give access to the data to the largest number of users, and it is easy to register on http://prisma-i.it/ index.php/en/ (last accessed on 15 June 2021) to obtain data (PRISMA Data Policy Issue 1.0 Date 20 May 2020) which can be provided with various levels of pre-processing (L0, L1, L2B, L2C, L2D) [52]. The hyperspectral camera covers VNIR (visible/near-infrared) and SWIR (short-wave infrared) with respectively 92 and 157 channels [35]. As this paper focuses on CH 4 quantification, only the 113 channels between 1500 and 2500 nm are retained. Indeed, CH 4 has several vibrational absorption bands at wavelengths of approximately 1650 nm and 2350 nm (plus 3400 nm and 7700 nm in mid-wave infrared). The CH 4 monochromatic absorption of the Pacific Northwest National Laboratory (PNNL) database [53], in ppm −1 ·m −1 , is shown in blue in Figure 1. This quasi-monochromatic absorption is then convolved with the PRISMA spectral response function (red line). In this spectral range, the PRISMA satellite has a spectral sampling distance of approximately 9 nm and a spectral resolution (Full Width of Half Maximum) between 9 and 15 nm [54]. The hyperspectral camera's spatial resolution is approximately 30 m (ground sampling distance) and an image covers an area of 30 km × 30 km. SNR specifications of PRISMA are described in Labate [35] and listed in Table 1. Using the Homogeneous Regions Division And Spectral De-Correlation (HRDSDC) method, Cogliati [55] estimated the PRISMA pixel SNR at each wavelength in the spectral dimension (by using the adjacent wavelengths) and scene heterogeneity (by using neighboring pixels in the same homogenous region). The results are provided in Table 1 and comply with the instrument's specifications. Table 1. Signal-to-noise ratio (SNR) specifications for the PRISMA satellite according to spectral range; requirement specifications come from Labate [35] and the observation-based specifications are estimated in Cogliati [55].

Spectral Range (nm) Requirement [35] Estimation [55]
SNR 1500-1750 ≥200 200-400 1950-2350 ≥100 ≈100 The images analyzed here focus on an industrial area, with four identified emission points in the west of Turkmenistan, as presented in Figure 2. These sources belong to the Korpezhe oil and gas field (38.52 • N, 54.21 • E). One of these emission points (D panel in Figure 2) is near a compressor station and has been identified as a strong CH 4 emitter [56]. On this industrial site, we were able to recover eleven images with cloudless sky or with few clouds between 2020/04/19 (YYYY/MM/DD) and 2021/06/22. Table 2 lists the dates these different images were taken.

Wind Information
Using discrete (in time or in space) gas concentration measurements, wind speed and direction are required to estimate the gas flow rate (see Section 3.4). In this paper, European Centre for Medium-Range Weather Forecasts (ECMWF) forecasts are used. ECMWF wind data was calculated at 00:00 or 12:00. The forecasts were then used to estimate the 10 m wind components for each hour with a grid of 0.125 • × 0.125 • , i.e., 13 km × 13 km. The Korpezhe site is located between four ECMWF grid points. The wind direction is expressed in a meteorological data frame. Angles of 0 • , 90 • , 180 • and 270 • correspond to wind blowing from the north, east, south and west, respectively. Due to the grid and time of ECMWF data, spatial and temporal interpolations had to be applied. Spatially, the Korpezhe values were obtained by distance-weighted linear interpolation. Temporally, all the images of the analyzed data set were acquired at around 7:20 UTC. The 7:00 and 8:00 UTC values were used to calculate image wind speed by a weighted average. Table 2 lists the wind speed U 10 and wind direction for the various PRISMA images.

Methods
The determination of the total mass of a CH 4 plume and the estimation of the source flow rate require a plume segmentation step followed by a quantification step, which have been implemented in the CELINE code. These two steps are based on the information contained in the radiance spectrum of the hyperspectral image. In the CELINE code, once the pixels affected by CH 4 are identified, the observed radiance is again used to calculate the CH 4 concentration for each plume pixel.

Direct Model of Radiance as a Function of Plume Concentration
In VNIR-SWIR, the radiance measured by an imaging instrument without any gas plume in the line of sight can be expressed as the sum of three different radiances (in W·m −2 ·sr −1 ·µm −1) (scheme on Figure A1 in Appendix A). The first radiance, L atm , is the atmospheric radiance without any interaction with the ground. The second one, L di f f , describes the photons that have been scattered towards the sensor after an interaction with the ground. The last one, L dir , represents the upwelling radiance directly transmitted from the ground. L di f f and L dir are both dependant on the downward solar irradiance and the ground reflectance [57,58]. Therefore, the total measured radiance L no pl (in W·m −2 ·sr −1 ·µm −1 ) can be expressed as: where R is the ground reflectance, E o the downward irradiance at ground level (in W·m −2 ·µm −1 ), τ dir the direct transmission, τ di f f the diffusive contribution, and S the spherical albedo of the atmosphere. The total atmospheric transmission τ atm is the sum of τ dir and τ di f f and can be expressed as τ atm = (1 + η)τ dir with η = τ di f f /τ dir . Scattering effects can be neglected in SWIR. A Radiative Transfer Model (RTM) (CO-MANCHE [59]) was effectively used to quantify the impact of each term. A sensitivity study was conducted to verify that additional simplifications could be made in the case of the airborne sensor. Firstly, the (1 − RS) term can be approximated by 1 with an error of approximately 0.05% for clear-sky observation conditions. Secondly, the diffusive contribution τ di f f is negligible compared to the direct transmission τ dir (η 1), which leads to an atmospheric transmission τ atm = τ dir .
In a near-nadir viewing geometry and taking into consideration the presence of a CH 4 plume in the line of sight, the atmospheric transmission when a plume is present can be expressed as τ dir = τ dir ·τ CH4 with: where ρ is the integrated line-of-sight concentration (in ppm·m), A the absorption (in ppm −1 ·m −1 ) at the instrument wavelengths, and θ the solar zenith angle. In the following, we replace the term of "integrated line-of-sight concentration" by the term "concentration (in ppm·m)". Therefore, the radiance in the presence of plume is L pl (see scheme on Figure A1 in Appendix A): Finally, between 1500 and 2500 nm, the atmospheric radiance represents less than 2% of the total radiance in clear-sky conditions. In the following, we assume that τ atm is equivalent to τ dir . These reasonable hypotheses result in the following expressions: In specific cases where the atmospheric radiance L atm cannot be neglected, a single RTM calculation of L atm for the whole image can be run (or derived from the image outside the plume) and the whole image can be corrected as follow: where L is the measured radiance and L * is the corrected radiance. In this case, we assume:

Plume Segmentation
The plume segmentation method consists of calculating a probability score, called the Clutter Matched Filter (CMF) score [43,44,60], used for hyperspectral image analysis. This technique estimates the probability of the presence of a spectral signature b (target gas signature) in the radiance spectrum L for each pixel, considering the covariance of the image background. Assuming that the background spectrum follows a normal distribution with a mean background spectrum L bkg and a covariance matrix C, the CMF method tests the most likely hypothesis between H 0 where the target gas is absent (no-plume pixel), and H 1 where the gas is present, e.g., The CMF method estimates the value of the parameter α that is quantifying the strength of the gas signature b in the background signal. The closer it is to the b-signal gas signature in the image, the greater the value of α. An improvement of the CMF method, quickly introduced by Funk [43] and widely applied in the literature, uses image clustering as pre-processing [23,29,[39][40][41][42]. Once the image is divided into K clusters (or classes thanks to K-means algorithm), an average background spectrum and a covariance matrix are calculated for each class. This improvement is called the Cluster-Tuned Matched Filter or CTMF method. The expression of the CTMF score α j for pixel (x, y) among z j pixels of the class j is: where the optimal filter q j of the class j is expressed as [43]: and L mean,j is the mean radiance of the class j and C j the class covariance matrix [61], such that: The gas signature b must be homogeneous to L(x, y) − L mean,j to which the optimal filter q j is applied (Equation (9)). The observed radiance for plume pixels corresponds to the L pl of the direct model (Equation (5)). An important hypothesis is the approximation of the no-plume radiance of Equation (4) by L mean,j to determine the expression of b (right term of the following equation): The CTMF score results in a detection map of pixels impacted by CH 4 . A threshold value of the score and morphological transformations are applied to isolate the CH 4 plume in the hyperspectral image (see Discussion in Section 5.3).

Quantification Step: Concentration
Once the mask of the plume pixels has been established, the CELINE code starts the quantification of the CH 4 plume in terms of concentration.

Linear Method (LM)
One of the most widely used methods for quantifying the plume concentration is based on Equation (12). The non-linearity of the problem is the main difficulty of the inversion. This method, which we call the linear method, uses a linearized expression for the methane transmission. Indeed, in the first-order Taylor series approximation, we can express τ CH4 as: In this case, the Equation (12) becomes: A linear inversion gives the concentration directly as: The LM is often presented as a method derived from the CTMF score (a derivative of the Matched Filter) [27,33,46,62]. Indeed, it is possible to express the concentration as the CTMF score divided by the optimal filter q j applied to the theoretical signal b j . When the hypothesis of low plume concentrations is no longer valid, an underestimation of the concentrations is observed as discussed below.

Optimal Estimation (OE) Method
To retrieve plume concentrations, OE is also commonly used in the literature [47][48][49]. The aim of this method was initially to retrieve the atmospheric composition by minimizing a cost function [24]. The standard OE approach uses iterative RTM computations to find the smallest difference between observed and modeled radiances [48]. In this approach, the state vector x contains all parameters that contribute to the radiance collected by the instrument, in particular the atmospheric composition and the ground reflectance. The inversion requires a model F(x) to fit the observed radiance y as: where is the sum of instrument and model errors. As the forward model is non-linear, the state vector is obtained iteratively. The Gauss-Newton solution is presented in Rodgers [24] as: In this equation, x i is the state vector at the iteration i, x a the prior estimation of the state vector, and S a the prior error covariance matrix. S is the covariance matrix of the image's pixels and K i = ∂F ∂x | x=x i the Jacobian matrix calculated for each iteration. This iterative inversion is applied to each plume pixel. OE inversion provides accurate results and has been used, for example, to determine the potential of next-generation imagers, such as PRISMA, to detect and quantify CH 4 point sources [37]. In this latter reference, the authors use the Iterative Maximum a Posterior Differential Optical Absorption Spectroscopy (IMAP-DOAS) algorithm [63,64] and a polynomial decomposition of the surface reflectance to characterize the algorithm's response to the surface type (grass, bright, dark or urban). In this case, the forward model is expressed as: where I 0 is the incident top-of-atmosphere (TOA) solar intensity, A the geometric air mass factor, τ n,l the default optical depth from the US Standard Atmosphere (from MOD-TRAN [65]) for trace gas element n = [1; 3] of the state vector at vertical level l = [1; 72], s n the scaling factor to that default optical depth optimized retrieval, P k the k th Legendre polynomial, and a k the coefficients optimized in the retrieval. However, iterative RTM computations for each pixel can be time-consuming and require good knowledge of nonretrieved parameters. In this paper, we combine the standard OE approach with a method to reconstruct the background (see below).

in-Scene Background Radiance (ISBR)-Optimal Estimation Method (OEM)
The method used in this paper is based on OE and the ISBR method developed for this study which can be referred to as ISBR-OEM. Nevertheless, the ISBR method can be applied to other algorithms, such as LM, to result in ISBR-LM.
The ISBR method is based on three assumptions. First, the radiance observed by the sensor in the presence of the plume can be written as L obs = L bkg ·τ CH4 (Equations (4) and (5)). Second, the background radiance L bkg can be rebuilt from the spectra contained in the image itself without RTM. Finally, the reconstruction of background radiance reduces the number of parameters to be retrieved to one, i.e., the state vector x of the standard OE approach only comprises the plume concentration ρ. The inverse problem's (Equation (16)) solution thus becomes: The main advantage of ISBR-OEM is that, unlike in previous papers, an RTM is not required, avoiding a complex and time-consuming step. The forward model F ISBR (ρ) is written as: For an in-plume pixel, background radiance cannot be obtained directly because this would require an image taken just before the start of the emission and temporally very close to the time of the plume image studied. Therefore, the ISBR method's background estimation steps are: (i) image classification using CH 4 non-absorbing channels, (ii) for each in-plume pixel, a set of no-plume pixels spectrally similar (corresponding to the same class as the target in-plume pixel) in the area surrounding the plume are selected; (iii) root mean square error (RMSE) is used to identify the "selected" pixel; (iv) in-plume pixel background radiance for the CH 4 absorbing channels is assumed to be the corresponding radiance of the "selected" pixel.
The a priori concentration x a = ρ LM is derived from the linear method even if the a priori has to be independent of the observation in the standard OE approach. To limit the share of observation in a priori knowledge, we used a prior error of 100% of the a priori concentration in the prior error covariance S LM . The a priori concentration is also the initial state x 0 . The ISBR-OEM's solution is then expressed as: The covariance matrix S is calculated for each class without using the pixels identified as belonging to the plume. The a posteriori uncertainty of the retrieved concentration is [47,48]. This a posteriori covariance matrix is providing a robust confidence interval for the concentration (for each plume pixel), from which a reliable uncertainty on the flow rate is estimated.

Quantification Step: Flow Rate
Once the wind information associated with the image is available, an estimation of the source flow rate is calculated by the CELINE code and we compare three methods of flow rate estimation.

Integrated Mass Enhancement (IME) Method
The IME method is commonly used to calculate flow rates from a gas concentration map [27,34,36,51], and it does not use wind direction. The IME method relates the flow rate Q I ME to the plume's total mass IME (in g), the wind speed U (in m·s −1 ) and the characteristic size of the plume L (in m) as: where ΣC plume is the concentration of the plume-pixel (in g), and A plume the surface of the plume (in m 2 ). The characteristic size of the plume L is defined as the square root of the plume surface area [37]. F I ME is the IME factor (in g·m −1 ) similar to a mass per unit length of the following methods.

Cross-Sectional Flux (CSF) Method
The CSF method was developed to be efficient computationally and conceptually simple to estimate a source's ozone and aerosol flow rates [66]. This method is commonly used in the literature to estimate gas flow rates [7,[67][68][69][70]. The product of the mass per unit length m ul CSF (or line density) by the wind speed U results in a flow rate value q at each distance d of the source as followed: To obtain the source's flow rate Q (in g·s −1 ), an average of q(d) is calculated over a distance range from the source. In our study, the wind speed grid is larger than the image dimension. The wind speed U (in m·s −1 ) is assumed to be constant in the entire image. Averaging q(d) is the same as averaging m ul CSF (d) (in g·m −1 ). Then, Q is defined as: where < > D is the mean operator between the source point and the distance D.
Estimating the principal wind direction is a critical point of the CSF method. The mass per unit length m ul CSF (d) is effectively obtained by summing the mass present on a slice perpendicular to the wind direction.

Rings Decomposition of Mass (RDM) Method
The local wind direction in the case of industrial plumes is not easy to estimate. Turbulence effects or the presence of obstacles (buildings, etc.) create variations in wind direction over very short timescales. In order to avoid this variability, a method independent of wind direction was developed during this work and based on a flow rate retrieval algorithm developed by SPASCIA society (see authors affiliation). After selecting the source pixel, the RDM method calculates the sum of the mass inside successive concentric rings (with the source point as centre) to obtain the mass per unit length m ul RDM (d) as a func-tion of the distance to the source. As previously, after spatial averaging of m ul RDM (d), a multiplication by the wind speed leads to the RDM flow rate Q RDM as: The radius of the rings is an input parameter. Choosing a small step (30 m or 1 pixel between successive rings) brings this method closer to the CSF method, while a step size covering the entire plume results in a calculation comparable to the IME method. The step has an impact on the uncertainty associated with the measurement.

LM versus ISBR-OEM
LM is the most common plume concentration quantification method. However, transmission linearization results in an underestimation of the retrieved concentration, especially for high concentration values. In order to illustrate this effect, we used a reflectance image using a land surface reflectance database and standard atmospheric parameters to create a radiance image (using the formula in Section 3.1) with a synthetic plume of known concentration (see Figure A3 in Appendix C) added to an atmospheric signal. The background reference signal corresponds to a standard atmosphere of MODTRAN code (mid-latitude summer atmospheric profiles and concentrations) [65]. This image is then processed by LM and ISBR-OEM. The results are plotted in Figure 3; the true and retrieved ISBR-OEM concentrations are shown as a function of the retrieved LM concentration. This particular graph compares the real case (PRISMA) and the synthetic case. The reflectance used in Figure 3a is based on the reflectance spectrum of a brick material (ECOSTRESS JPL Spectral Library at https://speclib.jpl.nasa.gov/ [71,72], last accessed on 25 January 2021) with synthetic noise according to PRISMA SNR (SNR of 150 at 2300 nm), while Figure 3b is based on a real reflectance map from a SWIR image over an industrial site [73]. The synthetic concentration curves (in red) are above the first bisector, meaning an underestimation of the concentration obtained with the LM as expected [73]. The higher the synthetic concentration, the higher the underestimation. Concerning ISBR-OEM (blue dashed line), the synthetic concentration is perfectly retrieved in the simple reflectance map case, and retrieval is fairly similar to the simulation for the real reflectance map despite scene heterogeneity and the interclass variability. In the test cases, the ISBR-OEM-retrieved concentration is equal to the input synthetic concentration, and the gap between ISBR-OEM and LM increases with concentration. The test cases give confidence in the results for the real case; the gap between ISBR-OEM and LM concentrations is almost the same with a real PRISMA image (green dots in Figure 3b) corresponding to the 2020/07/21 plume. As before, LM concentrations are lower than ISBR-OEM concentrations with the PRISMA image. Figure 4 corresponds to a spatial comparison of the plume concentration maps for the 2020/07/21 PRISMA image using LM (Figure 4a) and ISBR-OEM (Figure 4b). The absolute difference between ISBR-OEM and LM concentrations is shown in Figure 4c. The difference between the two methods is less significant in the low concentration parts than in the more concentrated parts of the plume. Nevertheless, the spatial smoothness of the plume with LM seems more realistic than with ISBR-OEM. This is mainly because in case of non-convergence of the OE, the concentration values are set to 0. The highest concentration pixel has a mass difference of 30%, which can be significant in terms of the estimated total mass emitted by the source and thus the flow rate. LM can be used up to a maximum concentration of approximately 6000 ppm·m to be less sensitive to underestimation (and for a signal sufficiently higher than the SNR). The analysis in terms of total mass and flow bias is presented for two plumes in Section 4.3. As the concentration map of the real plume is unknown, comparison between real and retrieved concentrations must be replaced by a comparative study of the measured and estimated radiances. Figure 5a represents (in red) the instrument-measured radiance for the highest concentration pixels near the source (approximately 20,000 ppm·m for the 2020/07/21 plume) for the wavelengths zoomed from 2000 to 2500 nm. The other colors portray LM and ISBR-OEM (green and blue, respectively) results, while other lines represent different types of radiance: the dashed lines represent the estimated background radiances L bkg used in the method, and the solid lines represent the instrument-estimated radiances calculated using the retrieved concentration value. The L bkg LM , i.e., the mean radiance of the class, is close but not equal to the L bkg ISBR−OEM obtained by the ISBR method. When these spectra are used to estimate the radiance with the associated retrieved concentration (respectively 16,000 ppm·m for LM and 20,000 ppm·m for ISBR-OEM), the ISBR-OEM spectrum is closer to the measured radiance than the LM spectrum. The residuals of the two methods are shown in Figure 5b. The residual's associated norm is 0.32 W·m −2 ·sr −1 ·µm −1 for the ISBR-OEM case compared to 0.39 W·m −2 ·sr −1 ·µm −1 for the LM case. The ISBR-OEM residual's associated norm is closer to the expected noise value in this spectral range for PRISMA data. This noise is estimated to be around 0.25 in the CH 4 absorption bands. ISBR-OEM avoids any linearization of the gas transmission. Using an estimation of the background radiance that varies for each pixel reduces residuals. The CH 4 signature has been largely suppressed in the residuals. Furthermore, ISBR-OEM benefits verification values from the standard OE approach, such as the degree of freedom or the χ 2 value. Examples are shown in Figure 6 for the 2020/07/21 PRISMA image. In the case of a state vector with one parameter, the concentration in this study, the degree of freedom can be interpreted as the share of the retrieved concentration coming from the measurement and independent to a priori information. This parameter must tend towards unity to ensure that the inversion is not over constrained by the a priori. Figure 6a shows a degree of freedom of approximately 0.99 in the most concentrated area and higher than 0.5 in the rest of the plume. The χ 2 value corresponds to the normalized mean square deviation between the observed and the modelled radiances for each wavelength. A χ 2 value close to unity means a correct estimation of retrieval uncertainty. Figure 6b shows a few pixels between 0.5 and 1, which can be interpreted as a small overcorrection, while a few others near the plume head are insufficiently corrected (χ 2 > 1). Nevertheless, the results remain close to unity which gives confidence to inversion results. In addition, the standard OE approach estimates a posteriori uncertainty [47]. For each plume pixel, retrieved concentration uncertainty is obtained in ppm·m, which can provide an estimation of the accessible concentration range. The a posteriori uncertainty for the plume in Figure 6 has an average value of approximately 350 ppm·m.

Overview of the Concentration Map
Here, only ISBR-OEM results are provided and discussed. Of the eleven images in the dataset, nine show a plume corresponding to source A (Figure 2), while the plume corresponding to source D is only seen in two images. Source B and C plumes are seen only once.
Independently of the sources and dates, thirteen plumes were detected. Focusing on the plume head, the plume from source A on the 2020/07/21 image, used as an illustration in the previous sections, shows the highest concentration with a maximum pixel of approximately 20,000 ppm·m. The maximum concentration values for each detected plume are listed in Table 3. To summarize, three plumes are marked for source A. Sources B and D each have one plume that appears significant, while the source C plume is less intense. For illustration, Figure 7 shows plume concentrations maps on the four sources. Plumes not presented in Figures 4b and 7 are shown in Figure A2 of Appendix B. A mask delimiting the plume contours was applied. These masks including the whole plume are called "full" masks. It is possible to see on this figure, especially on the plume of source D, estimation errors due to the terrain.
Indeed, these sources are distributed in a desert with roads and bright sand banks. The dispersion of a gas by wind advection and turbulence leads to a progressive decrease of the concentration as a function of the distance from the source. However, we observe a small local disturbance in this variation when the plume is covering road pixels. This disturbance is larger over a bright sand bank corresponding to high CTMF values (false alarms), probably due to the intra-class variability or "confuser" pixels that lead to a poor reconstruction of the background signal (see Discussion Section 5.4). The case of the 2020/10/10 image corresponding to source D illustrates the consequence of the presence of bright sand banks. The rectangle and the ellipse on the D source plume correspond to overestimated concentration areas and are linked to the presence of bright sand banks. If not considered, these artefacts have an impact on the estimation of the average mass and consequently on the estimation of the flow rate.

Mass per Unit Length or Wind Normalized Flow Rate Analysis
This section presents results for the flow rate quantification methods presented in Section 3 in terms of mean mass per unit length < m ul > and F I ME . The mean mass per unit length corresponds to the flow rate normalized by the wind speed (Section 2.2). The results are detailed in the case of ISBR-OEM retrieval and then compared to the LM retrieval. To analyze these methods, the most stable plumes were focused on. A stable plume must be as close as possible to a Gaussian plume, which settles in a constant emission regime and for which the mass spatial distribution follows a Gaussian model [74,75]. Thus, this analysis is based on two plumes of source A from the 2020/07/21 and 2021/06/22 images, shown respectively in Figures 4 and 7, and respectively referred to as P1 and P2 below.
The CSF and RDM methods are first considered. Figure 8 shows the spatial evolution of the mass per unit length of P1 as a function of the distance d to the source (located at d = 0) for the CSF and RDM methods shown in a bold line. The mass per unit length m ul (d) is obtained by summing the mass present on a slice of the plume (perpendicular to the wind direction for CSF and in the ring for RDM). The areas in blue (CSF method) and with orange dots (RDM method) correspond to the sum of a posteriori uncertainties from the OE approach. These uncertainties do not include the loss mass due to the detection threshold, which increases with the distance to the source. The CSF and RDM methods provide equivalent results as shown in Figure 8, but the RDM method does not use the wind direction estimation. The mean of the mass per unit length is calculated from these curves, and in order to minimize the errors due to the ground artefacts presented in Section 4.2, two masks have been established. The "small" mask covers the immediate environment of the source with an extension of 1400 m (excluding the road pixels or bright sand banks) and the "full" mask includes most of the concentration (see Figure 7). These masks have a direct impact on the mean of the mass per unit length and the associated uncertainty. The < m ul > uncertainty is lower with the small mask than with the full mask. Indeed, the associated uncertainty varies from 21.8 (23.4) g·m −1 to 36.8 (39.9) g·m −1 for the CSF (RDM) method from the small to the full mask, respectively. For the IME method, the uncertainty varies from 181.5 to 455.2 g·m −1 with the small and the full mask, respectively. As shown in Table 4, F I ME values are more than the double of CSF and RDM values due to the fact that the effective wind values to be considered are higher for the last two methods than for IME (see Section 5). The other cases are listed in Tables A1 and A2 in Appendix D. Table 4. Estimation of the wind normalized source flow rate corresponding to plumes P1 and P2 using ISBR-OEM retrieved concentration and CSF, RDM or IME flow rate estimation methods.

Plume
Mask <m ul CSF > (g·m − There is a difference in the estimated mean of the mass per unit length between the full and small masks, but the difference is relatively minor compared to the difference in surface area between the two masks (see number of pixels in the table). This means that the mean of the mass per unit length is relatively stable regardless of the area considered. It is therefore preferable to use the small mask for the flow rate estimation because the errors due to certain background elements are minimized. Table 5 is the equivalent of Table 4 when using the concentration obtained with the LM for flow rate estimates. The relative differences between the flow rate estimates from ISBR-OEM and LM concentrations (for the three flow rate methods) are equal to the relative difference between the plume total mass from the ISBR-OEM and the LM presented in Table 6). Here, this bias is estimated to be 16% (for P1) and 7% (for P2) in the case of a concentrated mask, and drops to 7% (for P1) and 4% (for P2) with the "full" mask.

Source Wind Normalized Flow Rate Estimation
We have shown that the ISBR-OEM can correct biases on pixel-to-pixel concentrations compared to a linear inversion. Tables 4-6 show, in a real case, the impact of the correction of this bias on the mass per unit length. This is smaller (5%) if all detected pixels are used for flow rate estimation. This value is directly related to the bias on the total mass included in the study mask. However, we have also seen that the uncertainty on the total mass increases when using a larger mask. It increases from 12.6% and 9.3% to 24.3% and 28.3% for P1 and P2, respectively. There is therefore a strong interest in using a "small" mask to reduce the uncertainties on the mass per unit length. In this case, it becomes necessary to correct the bias that can become greater than the uncertainty.

Wind Information and Flow Rate Estimation
Wind information and effective wind estimation uncertainties remain critical points for flow rate estimation. In this paper, ECMWF forecasts are used to derive the wind direction and wind speed. Due to spatial and temporal sampling of these forecasts, an interpolation was performed on the Korpezhe site.
Regarding the direction, when estimating the flow rate using the CSF method, we found, depending on the case, a significant difference between the wind direction given by the forecast, and that given by the main direction of the plume as seen in the images. As such, the plume direction used in this method was determined from the images themselves.
The wind speed U is derived from the 10 m ECMWF wind components. The spatial and temporal averaging and a 50% random error (1-69% for some authors [46]) added to the forecast derived values resulted in low confidence in the value of U [36]. Other authors suppose a 1.5 m·s −1 standard deviation in the wind speed [7,76]. In addition, the value used to calculate the flow rate is an effective wind speed, and a function of the 10 m wind speed. In the literature, authors used fluid dynamics simulations to calibrate the relationship between the 10 m wind speed and the effective U involved in plume dispersion [7,27,34,37,68,76]. Large-eddy simulation methods were used to determine the relationship, which depends on the flow rate estimation method, instrument specifications, and measurement conditions. The various expressions for U, according to the 10 m wind speed U 10 , bring uncertainty to the wind speed value. For example, regarding GHGSat-D, Varon (2018) [68] proposes U e f f I ME = α logU 10 + 0.6 m·s −1 with α between 0.9 and 1.1 and U e f f CSF = 1.47 U 10 m·s −1 . For the same instrument, Varon (2019) [34] uses U e f f I ME = log(U 10 ) + 0.5 m·s −1 for the IME method and U e f f CSF = 1.5 U 10 m·s −1 for the CSF method. For PRISMA, LES expresses effective wind speed as U e f f I ME = 0.34 U 10 + 0.44 m·s −1 [36]. Considering U 10 = 3 m·s −1 , the effective wind speed can vary from 2.05 to 2.37 m·s −1 for the IME method and is approximately 7.4 m·s −1 for the CSF method. Figure 9 represents the flow rate estimation for source A and uncertainties for each method in nine images from 2020 to 2021 with the ISBR-OEM retrievals. Figure 9. Temporal evolution of source A's flow rate for the three methods: CSF (blue), RDM (orange) and IME (red). The points of the three colors represent the flow rate obtained with the mean mass per unit length and the mean effective wind (average of all the U e f f calculated with U 10 from Table 2). This mean effective wind speed can be interpreted as the wind value assuming no uncertainty. The areas in blue, in orange dotted and in red correspond to the uncertainty on the mean mass per unit length and the same mean effective wind speed. The arrow error bars cover maximum flow rate value uncertainty. They are calculated by combining the uncertainty on the mean mass per unit length and the extremes of U e f f derived from U 10 ± 1.5 m·s −1 .
The flow rate obtained with the three methods are of the same order of magnitude with lower values for the IME method (around 30% lower than CSF). Taking only into account the uncertainty on the mean mass value and considering a perfect estimation of the effective wind speed, the uncertainty is reasonable as shown by the blue, red and orange (dotted) areas. The associated uncertainties are about 122, 135 and 196 g·s −1 for the CSF, RDM and IME methods, respectively. However, effective wind speed uncertainty has a significant impact on these values. Adding this U uncertainty results in the error bar value in Figure 9 that becomes rather large (863, 903 and 691 g·s −1 with the same order). Nevertheless, considering all these uncertainties, an overlap appears between IME and CSF flow rate estimates. The uncertainties are consistent with each other. If only one value of the source flow rate were to be given, an average of the retrieved values should be calculated because, statistically, the true value of the flow rate must lie between the retrieved values (range that contains the uncertainties of both estimation methods).
The effective wind is the largest source of uncertainty. The use of better quality wind products associated with the measurement should be feasible and would reduce the uncertainty associated with the estimated source flow value [77]. In the event that such wind products are not available, it is also possible to use methods that circumvent the need for wind speed information by using LES, providing a more robust uncertainty [51].

Manual Intervention
The eye is a very powerful tool to identify an area of interest to focus on. In the algorithm set up to shift from a hyperspectral image to a map of plume concentration and then flow rate, some steps require manual checking. This is the case for the definition of the plume contour in the detection part and then the selection of an area of interest to translate the concentration map into a flow rate. Selecting a plume area is a common problem in the literature [38,78]. Monitoring methane emissions using satellite images has the advantage of a short revisit time and global coverage of the Earth. Manually adapting processing to each image is not optimal for high-frequency or large-scale monitoring. These manual operations can be removed by optimizing the algorithm, which should be possible.

Overestimation Area
In some images and for some areas, an overestimation of the methane concentration was identified over specific areas in the desert region under study. This may be related to the presence of transport infrastructure as well as to spatially homogeneous sand banks with higher reflectance. Overestimations like this should not occur due to the classification that correctly groups these pixels into specific classes. A possible explanation is that, even if these pixels are grouped into a class, the step for reconstructing the background radiance using similar pixels around the source has limitations. If a few pixels of the same class surround the plume, then the most similar pixels selected may not be spectrally ideal. Further investigation is required to determine the cause of this local overestimation. However, some authors have shown that the surface reflectance features in SWIR can lead to unexpected retrieved results [79]. For example, calcium carbonate has already been identified as a "confuser", meaning a source of errors in plume concentration [64].

Conclusions
This paper shows a quantitative analysis of CH 4 plumes from an industrial site in the Turkmenistan desert based on PRISMA hyperspectral images at a spatial resolution of 30 m.
In open literature, a linear method based on CTMF score is used, providing a CH 4 concentration for each pixel in the hyperspectral image. However, this method uses a CH 4 transmission linearization that results in an underestimation of gas concentration. The proposed ISBR-OEM is divided into two steps. The first is the detection for which the CTMF score and standard image processing are used to delineate the plume. Detecting the plume's boundaries is required to perform the next step of quantification on a reduced area of the image, therefore improving efficiency. The plume quantification step is divided into two parts: plume concentration mapping and source flow rate estimation. Instead of using a RTM in OE, a radiance expression is developed from an estimation of background radiance and CH 4 transmission. ISBR-OEM uses an a priori concentration map derived from the map obtained by LM. ISBR-OEM has two advantages. In addition to having a relatively fast computing time (no RTM), a posteriori uncertainty and verification values (such as the degree of freedom) can be studied. The a posteriori uncertainty map provides a precise and controlled framework of uncertainty on the retrieved concentration value. We detected thirteen plumes in the eleven PRISMA satellite images acquired over the Korpezhe site from 2020 to 2021. Four sources were identified, and the major source (Source A, Figure 2) is seen in nine images, with the stronger plume associated with a surface area of 432 km 2 and a total mass of approximately 3.2 Tg, and with the highest concentration for a pixel with 20,010 ppm·m. In this paper, three flow rate estimation methods were compared against each other with consistent results given the associated uncertainties. A sensitivity study was conducted on the impact of the mask applied to the concentration map and on the use of a threshold value for this map. For source A, the higher flow rate was estimated at 3.9 ± 1.7 kg·s −1 using the CSF method, 4.0 ± 1.8 kg·s −1 using the RDM method and 2.6 ± 1.3 kg·s −1 using the IME method. The lower flow rate for the nine plumes was approximately 0.45 ± 0.35 kg·s −1 using the CSF method, 0.47 ± 0.35 kg·s −1 using the RDM method, and approximately 0.36 ± 0.3 kg·s −1 using the IME method. It has been shown that the ISBR-OEM reduces flow rate bias from 4% (plume P1) to 16% (plume P2).
The largest source of uncertainty in flow rate remains the effective wind speed estimation, which is responsible for more than 75% of the total uncertainty in the total budget error. Without wind uncertainty, the higher flow rate for source A was estimated at 3.9 ± 0.2 kg·s −1 using the CSF method and 2.6 ± 0.4 kg·s −1 using the IME method.

Data Availability Statement:
The data presented in this study is available on request from the corresponding author. The original data is not publicly available under ONERA's internal policy.

Conflicts of Interest:
The authors declare no conflict of interest. Appendix A Figure A1. Scheme of the radiances according to the direct model. The term E TOT describes the total solar irradiance received at the target pixel. It takes into account the direct and diffusive irradiances and the irradiance resulting from the ground-atmospheric coupling and can be written as E TOT = E 0 /(1 − RS).