Industrial Plume Properties Retrieved by Optimal Estimation Using Combined Hyperspectral and Sentinel-2 Data

: Stack emissions from the industrial sector are a subject of concern for air quality. However, the characterization of the stack emission plume properties from in situ observations remains a challenging task. This paper focuses on the characterization of the aerosol properties of a steel plant stack plume through the use of hyperspectral (HS) airborne remote sensing imagery. We propose a new method, based on the combination of HS airborne acquisition and surface reﬂectance imagery derived from the Sentinel-2 Multi-Spectral Instrument (MSI). The proposed method detects the plume footprint and estimates the surface reﬂectance under the plume, the aerosol optical thickness (AOT), and the modal radius of the plume. Hyperspectral surface reﬂectances are estimated using the coupled non-negative matrix factorization (CNMF) method combining HS and MSI data. The CNMF reduces the error associated with estimating the surface reﬂectance below the plume, particularly for heterogeneous classes. The AOT and modal radius are retrieved using an optimal estimation method (OEM), based on the forward model and allowing for uncertainties in the observations and in the model parameters. The a priori state vector is provided by a sequential method using the root mean square error (RMSE) metric, which outperforms the previously used cluster tuned matched ﬁlter (CTMF). The OEM degrees of freedom are then analysed, in order to reﬁne the mask plume and to enhance the quality of the retrieval. The retrieved mean radii of aerosol particles in the plume is 0.125 µ m, with an uncertainty of 0.05 µ m. These results are close to the ultra-ﬁne mode (modal radius around 0.1 µ m) observed from in situ measurements within metallurgical plant plumes from previous studies. The retrieved AOT values vary between 0.07 (near the source point) and 0.01, with uncertainties of 0.005 for the darkest surfaces and above 0.010 for the brightest surfaces.


Introduction
Human activities are responsible for significant emissions of particulate matter (PM). Atmospheric PM (also known as aerosols) is responsible for air quality degradation and the exposure of populations to ambient fine aerosols (PM 2.5 ), thus posing a global health concern [1,2]. The industrial sector was the second largest source of primary coarse aerosol (PM 10 ) emissions and the fourth largest source of PM 2.5 in Europe between 2013 and 2018. According to the emissions inventory of the national emissions reported at the Convention on Long-range Transboundary Air Pollution [3], PM 10 and PM 2.5 from the industrial sector represented 15% and 6.5% of the total emissions in Europe over the period 2003-2018, respectively. Emission inventories can be generated from on-site measurements, modelling, and reporting by the industry sector. Industrial emissions depend on manufacturing processes, fuel, emission factors, and technological abatement. Industrial emission inventories are updated each year according to specific guidelines [4,5]. However, the temporal variability of instantaneous emission fluxes remains a challenge a single hyperspectral image or a combination of the hyperspectral and multi-spectral Sentinel-2 observations. The plume aerosol optical properties are first estimated using a sequential method and, then, an iterative optimal estimation method, taking into account the observations and forward model uncertainties.

Theory
A radiance signal measured by an airborne or satellite sensor on a flat, homogeneous, and Lambertian surface can be expressed [41] as: where L = at-sensor radiance in W·m −2 ·sr −1 ·µm −1 L atm = atmospheric radiance, without interactions with the ground E d = direct part of the solar irradiance E s = scattering part of the solar irradiance T d = direct part of the atmospheric transmittance T s = scattering part of the atmospheric transmittance S = atmospheric spherical albedo ρ = surface reflectance of the studied pixel ρ e = surface reflectance from the environment.
In the presence of a homogeneous semi-transparent plume between the sensor and the soil, L atm , E s , E d , T s , and T d are modified. The total solar irradiance is given by E tot = E d + E s , and the total atmospheric transmittance is given by T tot = T d + T s . The pixel environment is assumed to be homogeneous, which means that ρ e is equivalent to ρ. The surface reflectance, ρ, is estimated by applying homogeneous atmospheric correction to the radiance image. The variational radiance signal due to the plume aerosols ∆L p is modelled by: where L p is the radiance signal measured for a single pixel in the same conditions as L, but with the presence of a plume, ∆L atm is the atmospheric radiance variation due to the plume aerosol, and ∆(E tot T tot ) is the variation of the product of the total irradiance with the total atmospheric transmittance. Philippets et al. [10] and Alakian et al. [42] have considered the plume as an infinite horizontal homogeneous layer. The plume is assumed to be here a finite horizontal homogeneous layer; this means that the direct part of the solar irradiance seen by a pixel under the plume may or may not pass through the plume. The variation of the scattering part of the solar irradiance seen by a pixel under the plume is a fraction of that seen in the case of an infinite plume layer. The total variation of the solar irradiance due to the plume and seen by a pixel can be expressed as follows: where α is equal to 0 or 1, β is a scalar between 0 and 1, ∆E d is the direct part of the total irradiance variation, and ∆E s is the scattering part of the total irradiance variation in the case of an infinite plume layer. In the same way, the surface reflectance variation, ∆ρ p , due to the plume aerosol is: where ρ p is the surface reflectance impacted by the plume.

Definition of the Optimal Estimation Method
The retrieval method is based on the optimal estimation method defined by Rodgers [28]. The OEM problem can be straightforwardly represented by Equation (5). The measure y is explained through the forward model function, F, associated with the state vector x and a random noise . In this study, the state vector x corresponds to atmospheric and surface properties, while the forward model is a radiative transfer model.
In the linear case, the forward model can be expressed as F = Kx, where K is the Jacobian matrix containing the partial derivatives corresponding to the sensitivities of the direct model to the state parameters. In this form, the problem is ill-posed, as there are more unknowns than observations. To reduce the number of unknowns, prior knowledge of the state vector is used. The optimization method consists of minimizing the cost function to a global minimum. According to Rodgers [28], the cost function is expressed by Equation (6): where the first term of the equation represents the difference between the state vector x and the a priori state vector, given the a priori variance-covariance matrix S a ; while the second term represents the error between the forward model F(x) and the observations y, given the variance-covariance matrix of the observations S . The variance-covariance matrices S and S a represent the uncertainties of the observation and the a priori state vector, respectively. According to the formalism of Rodgers, S can be decomposed into: where S y is the variance-covariance matrix representing the uncertainties due to the sensor and S b represents the effects of unknowns; that is, parameters that have an impact on retrieval uncertainties but were not retrieved during the OEM processing. K b is the Jacobian matrix associated with the unknown parameters. Rodgers [28] theoretically defined an estimator of the state vector of the retrievals, as follows: The Jacobian matrix is also defined as the pseudo-inverse of the gain matrix, G, which represents the sensitivity of the estimated state vectorx to the observations y. The gain matrix G is computed by the following expression: The elements of G correspond to the partial derivatives of the estimated state vector x, with respect to the observation y. The gain matrix gives access to the averaging kernel matrix A = GK. The diagonal elements of A represent the degrees of freedom (DOF) associated with the state parameters. The trace of A gives the total DOF of the problem (i.e., the number of independent parameters). Finally, the rows of A give the sensitivity of the retrieved parameters to the true state. In the case of a perfect estimate of x, the matrix A would be the identity matrix. The true error corresponding to a random measurement noise is defined by: and the posterior distribution of the estimated parameters has a covariance matrix given by: The uncertainty given byŜ depends on the resolving error caused by the lack of resolution in the inverse process and measurement errors due to the measurement noise. The measurement error and resolving error covariance matrices are S m and S s , respectively, which are defined by:Ŝ

Forward Model Definition and State Vectors
The forward radiance model,ŷ, is defined by the following equation: where L (see Equation (1)) is the estimated radiance without the plume and ∆L p (see Equation (2)) is the radiance differential due to the aerosol plume. The forward radiance modelŷ is estimated using the radiative transfer model with the following input parameters: Atmospheric background properties (water vapor, background aerosol properties), the finite plume geometry parameters (i.e., the direct and scattering parts of the irradiance differential, defined as α and β in Section 2.1), the aerosol plume refractive index R index , the plume modal radius r, the log-normal size distribution standard deviation σ, and the plume aerosol optical thickness (AOT) at 550 nm τ 550 . We assume that ∆L p is linear, with AOT in the range of [0, 0.5] [10].
where γ(τ 550 ) represents the AOT ratio between a reference AOT τ 550 re f at 550 nm and the observed AOT τ 550 at 550 nm. The atmosphere modelling was performed with the COMANCHE software [43], a frontend of MODTRAN model version 5.2 [44]. The plume was modelled by a homogeneous layer with a defined vertical extent, a height above the ground level, and an aerosol optical thickness τ 550 re f at 550 nm, defined by the user. The aerosol plume optical properties were simulated using Mie theory, considering a mono-modal size distribution for a given modal radius r, with an associated standard deviation σ set to 1.5. The simulations were performed for 3 different refractive indices R index , corresponding to sulphate, brown carbon, and soot particles ( Table 1).
The plume reference AOT was set to 0.1 at 550 nm (optically thin plume). The plume had a vertical extent of 100 m and was located 10 m above ground level. The vertical extent of the plume was within the common range of industrial stack plume extension [8]. Moreover, the plume reference AOT and vertical extent values have no impact on the retrieved plume properties, as long as the AOT value is less than 0.5 [10].
The state vector x is associated with the OEM retrieved parameters, with a prior distribution x a , and with a prior variance-covariance matrix S a . The state vector x includes: (i) The surface reflectance vector ρ, (ii) the plume AOT at 550 nm τ 550 , and (iii) the aerosol plume modal radius r. The prior distribution of the retrieved parameters and their variancecovariance matrix are given by: Pixel-by-pixel estimation of the prior value is defined in the next section. The covariance matrix, S ρ , of the surface reflectance ρ a was computed class-by-class. The covariance matrix S τ and S r were estimated from the standard deviation of the prior values τ 550 a and r a , respectively, and from the forward model sensitivity analysis (see Section 2.5).
The state vector b corresponds to the forward model inputs that are not retrieved by the OEM. The associated variance-covariance matrix, S b , represents the variance of the error of those unknown parameters. The state vector b includes the atmospheric water vapor content, background aerosol Ångström coefficient, and background atmospheric visibility. Visibility is related to ground surface aerosol extinction at 550 nm by the Koschmeider equation, and was used to scale the MODTRAN aerosol extinction profile.

Forward Model Initialization: First Guess
This step aims to initialize b and x and their associated uncertainties. Atmospheric background properties were estimated using MSI surface reflectances and the MODTRAN aerosol database. The MODTRAN background aerosol models were of three types: "maritime", "rural", or "urban". The rural aerosol model was composed of a mixture of 70% water soluble substances and 30% dust-like aerosols. The urban aerosol model was a mixture of rural aerosol particles and soot-like aerosols produced by industrial activities. The proportions of rural aerosols and soot-like aerosols in the urban aerosol model were 80% and 20%, respectively. The maritime aerosol model was composed of sea-salt particles created by the evaporation of sea-spray droplets and of continental aerosols, which were almost identical to those of the rural model. The exception was that the largest particles in the rural aerosol model were removed. Each aerosol model was defined by the normalized extinction and absorption coefficients at 550 nm. Figure 1a,b presents the normalized extinction and absorption coefficients for each model.  The state vector b and the associated S b matrix were initialized at this stage. The water vapor uncertainty (see Section 4.1) was set to 10% of the initial concentration of the mid-latitude winter profile defined by MODTRAN. Channels above 920 nm that were strongly affected by water vapour were not considered. The visibility error was fixed at 15%, corresponding to an absolute error of 5 km (see Section 4.1). The error in the background aerosol model was defined using the Ångström coefficient Å [45]. The Ångström coefficients were 0.48, 1.32, and 1.15 for the "maritime", "rural", and "urban" aerosol models, respectively. The standard deviation of the error associated with the Ångström coefficients was set to 10%. Surface reflectances were estimated by compensating for the atmospheric effects in the same way for the whole image, including the plume area. An extra step was needed to estimate the "off plume" surface reflectance in the plume area. In the case of a hyperspectral image alone, the "off plume" surface reflectance was assumed to be the estimated average spectrum, class-by-class. The images were classified with a random forest algorithm. The classification was performed with six user-defined classes: "water", "sparse vegetation", "dense vegetation", "concrete soils", "dark soils", and "bright soils". The training data sets were polygons drawn by the user using the open source geographical information system QGIS [46].
In the case of a single hyperspectral image associated with a multi-spectral image, different surface reflectance estimation methods may be considered. Different fusion algorithms exist, such as MAP-SMM [47] or FUSE [48], which are based on Bayesian approaches, or methods like Hysure [49], ICCV'15 [50], or CNMF [51], which are based on unmixing analyses. The coupled non-negative matrix factorization (CNMF) algorithm of Yokoya et al. [51] does not depend on image classification and the prior unmixing is unsupervised. CNMF merges a hyperspectral image with a multi-spectral image to obtain a final image, where each pixel is computed as a linear combination of the end members of the hyperspectral image. Vertex components analysis (VCA) [52] is used to obtain an initial set of endmembers for the hyperspectral image. The extracted hyperspectral endmembers and their weights are refined by alternating unmixing of hyperspectral and multi-spectral images by non-negative matrix factorization (NMF) [53,54]. Then, CNMF calculates the abundance matrix of hyperspectral endmember spectra using the multi-spectral image. The ground reflectance, ρ a (see Equation (16)), is then estimated, by combining the abundance matrix and the hyperspectral endmembers. The associated covariance matrix, S ρ , is then computed class-by-class.
First guesses of the AOTs and mean radii were estimated with a sequential approach based on the hyperspectral image corrected for the atmosphere. Atmospheres containing different types of plumes were simulated using the forward model. Three types of aerosol were considered (see Table 1). The standard deviation of the log-normal size distribution was equal to 1.5 and the modal radius varied from 0.025 µm to 1 µm. The sequential approach was performed using the surface reflectance differential defined in Equation (4), where ρ is the prior reflectance estimated with the CNMF method and ρ p is the hyperspectral data corrected for the atmosphere. The reflectance differential was compared to the simulations ∆ρ p simu (τ 550 re f ), using the cluster-tuned matched filter (CTMF) or the root-mean-square error (RMSE) as a metric. The analysis of the metric scores led to a first pixel-by-pixel estimation of the optical thickness, aerosol type, and modal radius of plume particles.  [55] and is used to retrieve gas thermal signatures (CO 2 , SO 2 , N 2 O, and O 3 ) in the thermal infrared. It was later adapted by Thorpe et al. [56,57] and Dennison et al. [58], in order to detect other gases in the reflective domain. Most recently, Philippets et al. [10] used CTMF to detect and characterize the aerosol signature of industrial plumes. The CTMF model can be described as a combination of the surface reflectance ρ with the spectral reflectance signature of the plume aerosols ∆ρ p . Both vectors ρ and ∆ρ p have the same dimension (i.e., the number of channels m in the hyperspectral image). After the classification of a hyperspectral image, an optimal filter q j for a soil class j can be defined as: where C −1 j represents the inverse of the correlation matrix of the soil class j and b j represents the mean aerosol spectral signature equivalent to ∆ρ p simu (τ 550 re f ). From the optimal filter q j associated to a modal radius and an aerosol type, the CTMF score, f i,j , describes the correlation between the simulated and the measured aerosol signal, which can be expressed as: where γ(τ 550 ) is the optical thickness ratio between the true optical thickness τ 550 and the reference optical thickness at 550 nm τ 550 re f , which was fixed to 0.1 in the direct model. A CTMF score map was computed for each aerosol radius and type. As the aerosol type was assumed to be homogeneous inside the plume, the aerosol type associated to the most frequent best score inside the plume was chosen for the entire plume. Once the aerosol type was selected, a second application of the CTMF was performed. The radius map was computed by selecting the best CTMF score in the set of CTMF maps for each pixel. The AOT was then deduced, by comparing the estimated surface reflectance differential map ∆ρ p (τ 550 ) with the simulated variation of the surface reflectance ∆ρ p simu (τ 550 re f ) corresponding to the selected modal radius r pixel-by-pixel: The RMSE defined by Equation (20) was used as an alternative metric to estimate the aerosol type, radius, and AOT.
where n is the number of channels λ of the spectrum and r is the modal radius. ∆ρ p simu was computed, by Equation (4), for each pixel and for the different aerosol properties. Aerosol properties were set pixel-by-pixel to the minimum RMSE value.

Sensitivity Analysis
Retrieval uncertainties, in terms of modal radius and AOT, were due to: (i) The signal-to-noise ratio (SNR); (ii) the ground reflectance error (retrieved by the OEM); and (iii) parameters not retrieved by the OEM (forward model assumptions). The OEM nonretrieved parameters were the atmospheric background properties (b state vector), the finite plume geometry parameters (i.e., the direct and scattering parts of the irradiance differential, defined as α and β in Section 2.1), the aerosol plume refractive index R index , and the aerosol plume size distribution standard deviation σ.
The instrument noise S y is a diagonal matrix with diagonal elements σ 2 l , representing the square of the noise equivalent delta radiance (NEDL). The NEDL σ l was computed by σ l = √ a 1 + a 2 L, where a 1 represents the residual noise and a 2 represents the photonic noise. The coefficients a 1 and a 2 were given by the noise model of the HYSPEX sensor. Figure 2 represents the modelled SNR of the HYSPEX sensor for the mean radiance of water pixels. The reference simulation state corresponded to seawater reflectance, a background aerosol model set to the MODTRAN rural mode with a visibility of 15 km. The scattering aerosol plume modal radius was set to 0.15 µm and the AOT to 0.1, respectively. The coefficients α and β were set to 0.0 and 0.3, respectively. The reference state was closed to the observation case. From the reference state, we estimated the retrieval uncertainties due to a nominal variation (or uncertainty) of the forward model inputs or assumptions. Estimated retrieval uncertainties corresponding to the HYSPEX noise were around 0.0024 for AOT and 0.001 µm for the modal radius. Table 2 shows the modal radius and AOT retrieval uncertainties due to forward model input uncertainties. Parameter uncertainties were defined consistently with the data and the initialisation step (Section 2.4). The surface reflectance uncertainty was set to 10% and uncertainties corresponding to the water vapor, the background aerosol Ångström coefficient, and the visibility used in the forward model were set to 10%, 10%, and 15%, respectively. An uncertainty of 10% was set for the standard deviation of the plume aerosol size distribution. Then, an uncertainty of 10% was set for the value of the real part of its refractive index. The finite plume geometry parameter uncertainties were set to 0.2 for β and to 1 for α.
We observed that most of the uncertainties due to the model assumption were higher than the measurement noise for both the AOT and modal radius. The maximum AOT error of 0.012 (i.e., a relative error of 12%) was due to a reflectance offset error: An error of 0.01 corresponds to more than 20% of water mean reflectance in the VIS spectral domain. As expected, the maximum modal radius error of 0.045 (i.e., a relative error of 30%) was due to σ uncertainty, as the retrieved modal radius is closely linked with the input size distribution standard deviation. The least sensitive parameters (errors with the same range of magnitude as noise error) were water vapor, β for AOT, and β and visibility for modal radius. In the case of other parameter uncertainties, retrieval error due to model uncertainties were quite similar, with a mean value around 0.005 (10%) for AOT and around 0.015 (10%) for modal radius. The forward model assumption uncertainties defined in this analysis led to acceptable uncertainties, in terms of the modal radius and AOT. Table 2. AOT and modal radius uncertainties corresponding to the uncertainties on the background aerosol model with the visibility vis and Ångström coefficient, the water vapor, the plume geometry α and β, the real part of the refractive index R index , the standard deviation of the plume aerosol size distribution, the surface reflectance relative error δρ, and the surface reflectance offset.

Airborne Hyperspectral Data
A plume emitted by the stack of the sinter plant at the ArcelorMittal site was imaged by the HYSPEX hyperspectral camera aboard the SAFIRE ATR-42 research aircraft on the 17th of February at 11:00 LT. The ArcelorMittal site is located in the Fos-sur-Mer district area and is the second largest steel plant in France, which produces steel, coils, and tubes. According to the IREP database (Registre français des rejets et des transferts de polluants), the Fos-sur-Mer ArcelorMittal industrial site was, respectively, the first and the second emitter of PM 10 (1230 t/year) and CO 2 (7,460,000 t/year) for France in 2018. The HYSPEX camera has a spectral range between 0.41 and 2.5 µm with 160 spectral channels in the visible/near-infrared (VNIR) range (410-996 nm) and 262 in the short-wave infrared (SWIR) range (970-2500 nm). The HYSPEX spatial resolution is 0.5 m × 1 m in the VNIR and 2 m × 2 m in the SWIR.
The georeferencing of the hyperspectral VNIR and SWIR images in UTM zone 31N was performed with QGIS, using a 50 cm spatial resolution ortho-rectified picture of the scene given by the Institut Géographique National (IGN) [59]. The georeferencing was performed with wedge points positioned manually on the IGN ortho-image. The VNIR image resolution was set to the SWIR resolution by using a spatial matching algorithm based on pixel matching through an optical flow calculation called GEFOLKI [60,61].

Sentinel-2 Products
Sentinel-2 (S2) was first launched in June 2015 (first platform), and a second platform was launched in March 2017. Sentinel-2 is equipped with the multispectral instrument (MSI), which acquires Earth atmosphere-reflected radiances in 13 spectral bands ranging from 442 nm to 2202 nm. The spatial resolution of MSI ranges from 10 m to 60 m, depending on the spectral channel. We selected an MSI image acquired on the 5th of February, 2016, being the closest S2 acquisition date to our airborne acquisition. The MSI image was processed using the MAJA algorithm [33]. The MAJA algorithm performs a masking procedure to detect and remove clouds and their shadows from the Level 1C image. The surface reflectance was derived following the algorithm proposed by Hagolle et al. [62]. Hereinafter, we used the level 2A surface reflectance product delivered by the processing chain of the THEAI center [63]. We used the channels B2 to B12, except for B9 and B10. The spatial resolution was set to 10 m for the selected bands. A first atmospheric correction was performed on the HYSPEX image, in order to apply a spatial matching algorithm (GEFOLKI) between MSI and HYSPEX surface reflectances to produce geolocalized images at a spatial resolution of 10 m. The transformation matrix, computed by GEFOLKI spatial matching, was then applied to the HYSPEX radiance image. Figure 3 shows a superposition of an RGB composite image of the HYSPEX image on S2 over the studied area. The stack and associated plume are located in the center of the image. The plume moves south-east over a water channel, crosses a ground arm with wharves, and then moves over the sea.
The scene is complex, with varying surface reflectances due to the presence of industries, water bodies, vegetation, and communication networks.

Atmospheric Correction
The optimal estimation procedure was applied to the off-plume pixels in the HYSPEX image, in order to estimate the atmospheric radiative terms. MSI surface reflectances were used to constrain the retrievals. However, MSI reflectances in the SWIR for water and dark vegetation were overcorrected, consequently leading to an overestimation of the visibility. SWIR channels had an influence on the OEM within the plume, as the uncertainty on the surface reflectance was very large. SWIR channels were, therefore, only used for the surface classification (see below), and not for plume retrieval. The most likely background aerosol models were the "maritime" and "rural" models, with a visibility of 18 ± 5 km and 15 ± 4 km, respectively. We selected the rural aerosol model for atmospheric correction. As mentioned in the method section, an error term on the aerosol model selection was introduced through use of the Ångstrom coefficient. The algorithm did not converge when using the urban background model, due to its much higher absorption coefficient (see Figure 1b). Figure 4 shows the surface reflectance classification results using the random forest procedure applied to the HYSPEX and MSI images. One can notice that a large part of the plume remained visible over water when using a single hyperspectral datum, while it disappeared when using the combination of HYSPEX and MSI. Moreover, the hyperspectral classification was more heterogeneous and misclassified "bright soils" or "sparse vegetation". As expected, the footprint of the plume was missing in the MSI classification. The classification of the MSI image was more homogeneous than that of HYSPEX; however, it over-classified spectra in the "concrete soils" class. This overclassification mainly concerns spectra that should have been classified as "bright soils".

Surface Reflectance
When considering the HYSPEX image alone, the spectra were averaged to provide a mean value for each class of soil. For the combination of the MSI and HYSPEX images, the surface reflectances were estimated using the CNMF method. In each case, the rootmean-square error (RMSE) and the spectral angle mapper (SAM) [64] were computed class-by-class (Table 3). The data were selected outside the plume area. For all the defined classes, the CNMF method reduced the heterogeneity of the classes and, thus, reduced the RMSE and the SAM. The scores of the combination or the single image were slightly similar for the water category, thus indicating that the averaged spectra can provide an estimate of the surface reflectance.

Plume Segmentation Using CTMF
The CTMF was used to identify the plume structure in the image. A mask was applied to the CTMF score corresponding to an aerosol modal radius of 0.2 µm. The 0.2 µm modal radius was chosen, in order to be of the same order of magnitude as the modal radius retrieved by Philippets et al. [10]. The objective was to have a first estimation of the plume footprint, in order to reduce the number of spectra to be inverted. Two masks were built: they corresponded, respectively, to the 5% and 30% of pixels with the best CTMF score. The thresholds were empirical. The pixels were then aggregated by region for each mask. Regions in the mask with the lowest threshold that did not intersect with the regions in the mask with the highest threshold were removed. The resulting mask was smoothed out with a 3 × 3 pixels median filter. Figure 5 shows the plume mask footprint obtained from the HYSPEX image alone and from the combination of HYSPEX and MSI. The identification of the plume structure was improved by using the combined HYSPEX/MSI product. Indeed, the number of false alarms was significantly reduced upstream of the source point and the plume was better defined downstream of the source point, compared to the plume segmentation from the HYSPEX image alone. The plume footprint was also visible over vegetation and over part of the artificial soils.

Plume Properties' Initialisation
In this stage, we applied the CTMF and RMSE sequential approaches to initialise the plume AOT, modal aerosol radius, and aerosol type pixel-by-pixel. Figure 6 shows maps associated with each metric, and Table 4 corresponds to the mean retrieved values class-by-class. The aerosol type most frequently obtained by the CTMF and the RMSE was a "scattering aerosol", as is expected from a sinter plant emitter [65,66]. Over the whole image, 70% of the pixels were classified as scattering aerosols for the RMSE method and 45% for the CMTF method. The retrieved modal radii differed between CTMF and RMSE. Over water bodies, the RMSE retrieved modal radii varying from 0.1 to 0.25 µm, being close to that found in previous studies over the same area Philippets et al. [10]. However, these values seemed to overestimate in situ measurement over metallurgic plants [65,67]. The CTMF provided higher aerosol radii (higher than 0.5 µm) over water bodies, where the AOTs of the plume were below 0.02, but provided similar radii as the RMSE method for larger AOTs in the densest part of the plume. Over artificial soils, the mean retrieved radius was about 0.4 µm for the CTMF and 0.25 µm for RMSE; additionally, 33% and 38% of the retrieved radii were higher than 0.5 µm for the RMSE and CMTF methods, respectively. As mentioned in Section 4.3, the surface reflectance correction was less consistent over artificial and bright soils. For the surface reflectances of the artificial and bright soils, both methods tended to overestimate the radius and AOT, compared to that for water surfaces, where the error of the estimated surface reflectance was low.
The retrieved mean AOTs were about 0.02 for both the CTMF and RMSE. The mean difference between AOTs retrieved was lower than 0.01 for all surface reflectance classes, except for the bright soil class (Table 4). For the bright soil class, the mean difference was 0.23.   Figure 7 represents the difference between spectra acquired inside and outside the plume (point B and point A in Figure 3, respectively). The point B spectrum corresponds to a mean of 3 × 3 pixels, having a mean AOT of 0.06. The spectra A and B were localised over the water surface and their differential was compared to the differential simulated using aerosol properties from the CTMF or RMSE methods. Mean radii retrieved using CTMF and RMSE were 0.13 µm and 0.16 µm, respectively. The normalized RMSE of the differential were 9.2% for the CTMF and 6.5% for the RMSE methods, and the corresponding SAM were 5.1 • and 3.6 • , respectively. Example of difference between spectra acquired inside and outside the plume over water surface (red), compared to the best differential model from the CTMF (orange) and RMSE (green) approaches.
In general, the differential calculated with the radius obtained by the RMSE method better corresponded to the observations for all the soil classes, although there were still high radius values for artificial soils. The CTMF formulation is well-adapted for gases with a signature defined by local absorption bands; however, due to the correlations between the spectral signatures of aerosols and their influence on the slope of the entire spectrum, the CTMF results led to a high error rate on the type of aerosol. The RMSE, therefore, seems to best provide a first estimate of the modal radius and the AOT. In addition, the RMSE provides an estimate of the α and β parameters. Minimization of RMSE in the densest part (i.e., at the center) of the plume corresponded to values of 0 for α and 0.3 for β.

Results of the Optimal Estimation Method
OEM (see framework in Section 2.2) accounted for the errors associated with the prior parameters, the observations, and the atmospheric correction. The OEM was initialized with the AOT and modal radius retrieved from RMSE sequential approach and with the surface reflectances under the plume estimated by the CNMF. AOT, modal radius, and surface reflectance maps, with their associated posteriori uncertainties and their DOF, were then retrieved. Figure 8a illustrates these results. Low DOF values indicate a retrieval being dependent on a priori information, meaning that, due to the large initial error (reflectance and/or aerosol properties), the residual was too high when compared to a priori uncertainties. Only pixels with a retrieved radius with DOF > 0.5 were used to post-process the retrieved maps of AOT and modal radius (see Figure 8). Pixels for which the OEM failed to converge were included in the rejection mask. This rejection mask was associated with pixels where prior information were not consistent with the measurement and the model. It corresponded, in particular, to pixels over artificial soils near the source point, where the prior radius from RMSE initialization seemed to be overestimated (values larger than 0.5 µm). Table 5 shows the mean radius retrieved by the OEM class-by-class. The mean retrieved radius in the entire plume was 0.125 µm, with a standard deviation of 0.05 µm. In comparison, inside the final plume mask, the initial mean radius was around 0.15 µm, with a standard deviation of 0.1 µm. This decrease in the mean radius led to a slight decrease in the AOT-in particular, for bright surfaces. Near the source, the mean retrieved radii were between 0.1 and 0.2 µm and the AOT was above 0.05. Downwind of the source point over water bodies, the mean retrieved radius was 0.12 µm, with a standard deviation of about 0.04 µm, while the retrieved AOT varied from 0.01 to 0.07. Over sparse vegetation soils, the DOF were lower than over water bodies, with values varying from 0.2 to 0.7; furthermore, the mean retrieved radius was around 0.15 µm, with a standard deviation of about 0.07 µm. Table 5 shows that the variation of the mean modal radius from one class to another was reduced from the initial guess. OEM-retrieved values were more homogeneous spatially than in the a priori. These results show that, over the whole plume, the size distribution was fairly stable-corresponding to an ultra-fine mode.  The residual errors associated with the retrieval for different surface categories are presented in Table 6. Figure 9 illustrates the mean standard deviations by class of soil of posterior uncertainties on the retrieved surface reflectances. The main part of the error came from the resolving error σ s (see Table 6), except for the water surface, where the error was driven by the measurement noise. Radius uncertainties from Table 6 were not really dependent on the class, with a mean value around 0.05; this value was coherent with the standard deviation estimated from the OEM radius map. However, radius prior uncertainties remained larger for bright surfaces and vegetation surfaces.
AOT posterior uncertainties depended on the type of surface, with values varying from 0.003 for water to 0.013 for bright soils. This means that the OEM surface reflectance retrieval did not compensate completely for the initial error in the surface reflectance estimation computed with the CNMF. As shown in Figure 9, the reflectance posterior uncertainty for bright surfaces, vegetation, and concrete soil were higher than for water reflectance spectra. For water surfaces, the mean posterior reflectance uncertainties were less than 0.01, which led to less uncertainty in the retrieved AOT. For bright surfaces, the mean reflectance uncertainties increased up to 0.045, leading to larger uncertainties ( Table 6). Finally, the OEM-retrieved plume parameters were less sensitive to ground heterogeneity than the initial sequential approach. The estimated mean modal radius was about 0.125 µm, with a standard deviation of about 0.05 µm. Moreover, the OEM helped to enhance the quality of the retrieval, by rejecting pixels where the prior information was not consistent with model uncertainties.

Retrieved Modal Radius
The size distribution retrieved in the plume corresponded to an ultra-fine log-normal mode, having a modal radius of 0.125 µm and a standard deviation of 1.5. The retrieved size distribution was in agreement with previous work based on airborne hyperspectral data on the same plant, which gave a modal radius around 0.2 µm [10,66]. In situ observations over two metallurgic sites in Europe [65,67] have shown that the particle size distribution inside the plume and near the source corresponds to an ultra-fine mode with a modal radius that is (most of the time) less than 0.1 µm, with a standard deviation around 1.5. Aerosols in the ultra-fine fraction are formed by secondary aerosols in the Aitken mode, which coagulate downwind the plume [65] and contain metal-bearing nanoparticles [67]. The contribution of fine particles (PM 1 ) can be up to 60% for metallurgy, although the PM 1 fraction depends on the industrial process [68]. A similar result has been obtained for a refinery plume using hyperspectral airborne imagery [11], which showed a retrieved particle size distribution corresponding to the ultra-fine mode (modal radius around 100 nm). An absorbing fraction of around 10% has been found in the refinery plume, while only a scattering aerosol was retrieved in the metallurgy plume.
Although the hyperspectral retrievals indicate the presence of ultra-fine particles, it is expected to have a contribution of larger particles (PM 10 ) inside the plume [68]. However, the sensitivity of VNIR radiances to the aerosol coarse mode with a modal radius higher than 1 µm is weak. Increasing the standard deviation of the size distribution (e.g., from 1.5 to 2.0) lowered the retrieved radius (by 0.05 µm, respectively) in areas with DOF > 0.5, indicating a competitive impact of the modal radius and the standard deviation on the retrieval. This confirms that, in the VNIR spectral domain, the observed plume signature has a very weak contribution of coarse aerosols. The addition of a coarse mode in the aerosol forward model, as associated with an extension of the retrieval spectral range to the SWIR domain, may help to infer the contribution of the different aerosol modes.

Ground Reflectance Estimation and OEM
The surface reflectance estimation accuracy in the initialization step has a major impact on the retrieval uncertainties. Particularly in the case of bright surfaces, any bias or offset on the estimated reflectances must be avoided. Using a single hyperspectral image to estimate the ground reflectance can lead to large errors, as was seen in Sections 4.2 and 4.3. The use of MSI data in the classification allows for a removal of the radiative impact of the plume (see Figure 4) from the background image and reduction of the intra-class surface reflectance variability, compared to a single hyperspectral acquisition (see Table 3). Moreover, adding the ground reflectance as a retrieved parameter in the OEM process decreases the reflectance uncertainties, compared to the initial stage. Consequently, the reflectance uncertainty can be reduced by a factor of 2 (see Figure 9) for bright and concrete soils. Considering a further application to hyperspectral space missions, such as PRISMA [69], the reflectance error estimation could be reduced, as: (i) the spatial heterogeneity will decrease, due to the lower spatial resolution (30 m instead 10 m); (ii) the signal-to-noise ratio is similar to HYSPEX noise (see Figure 2b); and (iii) the fusion algorithms will not have to deal with a change of spectral resolution. Additionally, we notice that the OEM also helps to reject the pixels which have an estimated ground reflectance error too high and/or with a false positive plume detection, thanks to the DOF and the algorithm convergence.

Model Assumptions and Uncertainties
From sensitivity and posteriori uncertainty analyses, we saw that the modal radius uncertainty of 0.05 µm was quite independent of the reflectance value, whereas the AOT error increased with the reflectance value (from 0.003 to 0.01). Moreover, the non-retrieved b state vector used in the results section did not include all the forward model assumption uncertainties or errors described in the method section. In particular, the aerosol plume refractive index uncertainty and size distribution standard deviation are important sources of error for the AOT (0.008) and modal radius (0.045 µm), respectively. This means that the a posteriori uncertainties estimated in the results section were underestimated, particularly for low reflectance values. We assumed that the real uncertainties were around 0.01 for AOT and around 0.06 µm for the modal radius. Some forward model assumptions are still missing in the total budget error, such as the environmental effects and the uncertainty in the imaginary part of the refractive index.

Conclusions
Stack plume detection and characterization over a metallurgic plant were performed with the use of airborne hyperspectral data. The proposed method relies on the combination of hyperspectral acquisition with additional Sentinel-2/MSI surface reflectances. The spatial resolution of the product was 10 m, matching the satellite data. The MSI data provided an off-plume background image of the scene and, so, reduced the error in both plume segmentation and the estimation of aerosol optical properties.
A first guess of the plume aerosol optical properties, in terms of modal radius and AOT, was performed using a sequential method. The previously proposed CTMF method was outperformed by the RMSE metric in the sequential approach. The OEM improved the consistency of the retrieval and provided an estimation of the error balance, by simultaneously retrieving the plume properties and the surface reflectances under the plume. OEM provides the opportunity to analyse uncertainties through an analytical equation system. It can provide information (gain matrix) and criteria (degree of freedom) on the retrieval sensitivity of plume properties above several ground types, with respect to the observation noise. The observed plume was optically thin, with AOT lower than 0.1 and a detection threshold around 0.01. The retrieved size distribution mode corresponded to a fine (i.e., ultra-fine) mode with a modal radius of 0.125 µm and a log-normal distribution standard deviation of 1.5. All of the detected pixels had a DOF higher than 0.5 and the average a posteriori uncertainties of the estimates were 0.01 for the AOT and 0.06 µm for the modal radius. These results were close to the ultra-fine mode (modal radius around 0.1 µm) observed from in situ measurements within the plumes of metallurgical plants [65,67].
Knowledge of PM atmospheric concentrations from industrial stack plumes and associated emission fluxes from satellite observations could allow for a better assessment of the impact of industrial emissions on air quality at a high spatial resolution. The proposed method to retrieve aerosol optical properties and evaluate the associated uncertainties serves as a first step toward this objective. We also showed the potential and feasibility of combining hyperspectral data with other satellite data with different spectral and spatial resolutions, but higher temporal frequency. The proposed framework can be applied to current hyperspectral satellite missions, in combination with the available Sentinel-2 products.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the ONERA internal policy.