Gabriel Calassou
email: gabriel.calassou@onera.fr

Pierre-Yves Foucher
email: pierre-yves.foucher@onera.fr

Jean-François Léon
email: jean-francois.leon@aero.obs-mip.fr

remote sensing

Keywords: aerosol, plume, hyperspectral, multi-spectral, stack emissions

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 reflectance imagery derived from the Sentinel-2 Multi-Spectral Instrument (MSI). The proposed method detects the plume footprint and estimates the surface reflectance under the plume, the aerosol optical thickness (AOT), and the modal radius of the plume. Hyperspectral surface reflectances 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 reflectance 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 filter (CTMF). The OEM degrees of freedom are then analysed, in order to refine 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-fine 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 [START_REF]Health Risks of Particulate Matter from Long-Range Transboundary Air Pollution[END_REF][START_REF] Burnett | Global Estimates of Mortality Associated with Long-term Exposure to Outdoor Fine Particulate Matter[END_REF]. 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 [START_REF]National Emissions Reported to the Convention on Long-Range Transboundary Air Pollution (LRTAP Convention)-European Environment Agency[END_REF], 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 [START_REF]EMEP/EEA Air Pollutant Emission Inventory Guidebook 2019: Technical Guidance to Prepare National Emission Inventories[END_REF][START_REF] Andre | Rapport Ominea[END_REF]. However, the temporal variability of instantaneous emission fluxes remains a challenge to assess; in particular, in the case of unexpected atmospheric emissions due to industrial accidents or uncontrolled releases.

Ground-level sensor networks, including remote sensing-based techniques, have been commonly used for fine particle monitoring [START_REF] Amaral | An Overview of Particulate Matter Measurement Instruments[END_REF]. Since the signing of the Convention of Long-range Transboundary Air Pollution in 1979, several monitoring programs have been deployed in Europe, such as the European Monitoring and Evaluation Programme (EMEP). However, such ground-based networks are not always able to appropriately capture the spatial dispersion of the plume from the emission source. Airborne and satellite remote sensing imagery may, therefore, be considered as complementary to ground-based systems.

The dilution of a stack plume in the ambient air depends on the wind field, the atmospheric stability [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], the effluent buoyancy, and momentum [START_REF] Gordon | A Comparison of Plume Rise Algorithms to Stack Plume Measurements in the Athabasca Oil Sands[END_REF]. With dilution, the optical signature of the plume reduces in amplitude. Industrial plumes from an industrial stack generally have a rather limited visible spatial extent, only about a few hundred meters. Detecting a plume in a remote sensing image remains a challenge, and it is even more difficult to quantify the aerosol properties inside the plume. Hyperspectral visible/short-wave infrared (VIS-SWIR) imagery with a high spatial resolution (of about 1-10 m) is well-suited for the plume detection task and the study of aerosol properties [START_REF] Alakian | Retrieval of Microphysical and Optical Properties in Aerosol Plumes with Hyperspectral Imagery: L-APOM Method[END_REF][START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF][START_REF] Foucher | Aerosol Plume Characterization from Multitemporal Hyperspectral Analysis[END_REF]. Moreover, recent and planned high spectral and spatial resolution space missions have offered new opportunities for the study of stack emissions on a global scale.

Remote sensing of atmospheric aerosols has benefited from space missions dedicated to Earth observation (e.g., the NASA EOS program). First considered as a corrective parameter for estimating surface reflectances, aerosols have become an important focus of research, due to their important role in climate change. Kaufman and Sendra [START_REF] Kaufman | Algorithm for Automatic Atmospheric Corrections to Visible and Near-IR Satellite Imagery[END_REF] and Kaufman et al. [START_REF] Kaufman | Passive Remote Sensing of Tropospheric Aerosol and Atmospheric Correction for the Aerosol Effect[END_REF] have used the atmospheric correction approach for aerosol characterization over dark, densely vegetated pixels. This method is based on dedicated spectral channels or channel ratios of multi-spectral sensors and, since then, has been improved and successfully applied to multi-spectral images [START_REF] Levy | Global Evaluation of the Collection 5 MODIS Dark-target Aerosol Products Over Land[END_REF][START_REF] Remer | The Dark Target Algorithm for Observing the Global Aerosol System: Past, Present, and Future[END_REF]. Variants of this method use channels with wavelengths below 500 nm over bright soils [START_REF] Hsu | Aerosol Properties Over Bright-Reflecting Source Regions[END_REF][START_REF] Hsu | Deep Blue Retrievals of Asian Aerosol Properties during ACE-asia[END_REF][START_REF] Sayer | Stability, and Consistency of MODIS Collection 6.1 and VIIRS Version 1 Deep Blue Aerosol Data Over Land[END_REF].

Inverse model analysis applied to aerosol satellite remote sensing is an ill-posed problem, for which statistical frameworks, such as optimal inversion [START_REF] Thomas | The GRAPE Aerosol Retrieval Algorithm[END_REF][START_REF] Dubovik | Statistically Optimized Inversion Algorithm for Enhanced Retrieval of Aerosol Properties from Spectral Multi-Angle Polarimetric Satellite Observations[END_REF][START_REF] Wurl | Optimal Estimation Retrieval of Aerosol Microphysical Properties from SAGE~II Satellite Observations in the Volcanically Unperturbed Lower Stratosphere[END_REF][START_REF] Zheng | Optimal Estimation Retrieval of Aerosol Fine-Mode Fraction from Ground-Based Sky Light Measurements[END_REF][START_REF] Govaerts | Joint retrieval of surface reflectance and aerosol optical depth from MSG/SEVIRI observations with an optimal estimation approach: 1[END_REF], have been proposed. Thompson et al. [24,[START_REF] Thompson | Optimal Estimation of Spectral Surface Reflectance in Challenging Atmospheres[END_REF] and Hou et al. [START_REF] Hou | An Algorithm for Hyperspectral Remote Sensing of Aerosols: 1. Development of Theoretical Framework[END_REF][START_REF] Hou | An Algorithm for Hyperspectral Remote Sensing of Aerosols: 2. Information Content Analysis for Aerosol Parameters and Principal Components of Surface Spectra[END_REF] have presented a method based on the formalism of Rodgers [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding[END_REF], in order to simultaneously retrieve the atmospheric properties and the surface reflectance from VIS-SWIR images. Several studies have used the formalism of Rodgers to retrieve gas, aerosol, or surface properties from space missions such as SAGE II [START_REF] Wurl | Optimal Estimation Retrieval of Aerosol Microphysical Properties from SAGE~II Satellite Observations in the Volcanically Unperturbed Lower Stratosphere[END_REF], IASI [START_REF] Kangah | Measurement of Nitrous Oxide (N 2 O) from Space[END_REF][START_REF] Hurtmans | FORLI Radiative Transfer and Retrieval Code for IASI[END_REF][START_REF] Barret | The Detection of Post-monsoon Tropospheric Ozone Variability Over South Asia Using IASI Data[END_REF], and GOSAT [START_REF] Yoshida | Retrieval Algorithm for CO 2 and CH 4 Column Abundances from Short-wavelength Infrared Spectral Observations by the Greenhouse Gases Observing Satellite[END_REF]. Sequential algorithms, such as MAJA [START_REF] Hagolle | Algorithm Theoretical Basis Document[END_REF] or GRASP [START_REF] Dubovik | GRASP: A Versatile Algorithm for Characterizing the Atmosphere[END_REF], allow for the retrieval of atmospheric properties with a multi-temporal comparison and a combination of several images. Sequential atmospheric correction algorithms have been used on hyperspectral data [START_REF] Qu | The High Accuracy Atmospheric Correction for Hyperspectral Data (hatch) Model[END_REF][START_REF] Adler-Golden | Atmospheric Correction for Shortwave Spectral Imagery Based on MODTRAN4[END_REF][START_REF] Staenz | ISDAS-A System for Processing/analyzing Hyperspectral Data[END_REF][START_REF] Richter | Atmospheric Correction of DAIS Hyperspectral Image Data[END_REF][START_REF] Richter | Correction of Satellite Imagery Over Mountainous Terrain[END_REF][START_REF] Richter | Geo-atmospheric Processing of Airborne Imaging Spectrometry Data. Part 2: Atmospheric/topographic Correction[END_REF].

The retrieval of the aerosol optical properties of an industrial stack plume requires both surface and atmospheric correction. The plume properties can be estimated by using the difference between a reference image without the plume and the actual image with the plume. However, while atmospheric correction is usually homogeneous over the scene, the surface properties can be highly heterogeneous. Therefore, an adequate method is needed to infer the surface reflectance under the plume. The reference signal can be estimated from pixels out of the plume, through using a mean spectrum considering the class of surface reflectance [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF] or by using multi-temporal data [START_REF] Foucher | Aerosol Plume Characterization from Multitemporal Hyperspectral Analysis[END_REF]. The method using the mean spectrum by class is appropriate for rather homogeneous classes. The reference can also be estimated from ancillary radiometric observations coming from other sensors, such as multi-spectral images, as is presented in this study.

In this paper, a method to improve the estimation of plume aerosol properties from hyperspectral image is described. The proposed framework is applied to a plume observed in February 2016 over a steel plant in France. The method is based on multi-temporal and multi-resolution observations, as well as an optimal estimation method (OEM). We first describe how the surface reflectance under the plume is determined, using either 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.

Methods

Theory

A radiance signal measured by an airborne or satellite sensor on a flat, homogeneous, and Lambertian surface can be expressed [START_REF] Chandrasekhar | Radiative Transfer[END_REF] as:

L = L atm + ρ (E d + E s )T d π(1 -ρ e S) + ρ e (E d + E s )T s π(1 -ρ e S) , (1)
where 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:

L = at-sensor radiance in W•m -2 •sr -1 •µm -1 L atm =
∆L p = L p -L = ∆L atm + ρ ∆(E tot T tot ) π(1 -ρS) , (2)
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. [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF] and Alakian et al. [START_REF] Alakian | Remote Sensing of Aerosol Plumes: A Semianalytical Model[END_REF] 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:

∆E tot = α∆E d + β∆E s , ( 3
)
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:

∆ρ p = ρ p -ρ = L p -L atm E tot T tot + πS(L p -L atm ) - L -L atm E tot T tot + πS(L -L atm ) , (4)
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 [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding[END_REF]. The OEM problem can be straightforwardly represented by Equation [START_REF] Andre | Rapport Ominea[END_REF]. 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.

y = F(x) + . ( 5
)
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 [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding[END_REF], the cost function is expressed by Equation ( 6):

χ 2 = 1 2 (x -x a ) T S -1 a (x -x a ) + 1 2 (y -F(x)) T S -1 (y -F(x)), (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:

S = S y + K b S b K T b , (7)
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 [START_REF] Rodgers | Inverse Methods for Atmospheric Sounding[END_REF] theoretically defined an estimator of the state vector of the retrievals, as follows:

x = x a + (K T S K + S -1 a ) -1 K T S -1 (y -Kx a ), (8)
The Jacobian matrix is also defined as the pseudo-inverse of the gain matrix, G, which represents the sensitivity of the estimated state vector x to the observations y. The gain matrix G is computed by the following expression:

G = (K T S -1 K + S -1 a ) -1 K T S -1 . ( 9
)
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:

xx = (GK -I)(xx a ) + G , [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF] and the posterior distribution of the estimated parameters has a covariance matrix given by:

Ŝ = (K T S -1 K + S -1 a ) -1 . ( 11
)
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:

Ŝ = GS G T + (A -I)S a (A -I) T (12)
= S m + S s . (13)

Forward Model Definition and State Vectors

The forward radiance model, ŷ, is defined by the following equation:

ŷ = L + ∆L p (τ 550 ), (14)
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] [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF].

∆L p (τ 550 ) = γ(τ 550 )∆L p (τ 550 re f ), [START_REF] Remer | The Dark Target Algorithm for Observing the Global Aerosol System: Past, Present, and Future[END_REF] 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 [START_REF] Poutier | Two reciprocal atmospheric codes for hyperspectral remote sensing[END_REF], a frontend of MODTRAN model version 5.2 [START_REF] Berk | MODTRAN® 6: A Major Upgrade of the MODTRAN® Radiative Transfer Code[END_REF]. 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 [START_REF] Gordon | A Comparison of Plume Rise Algorithms to Stack Plume Measurements in the Athabasca Oil Sands[END_REF]. 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 [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF].

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:

x a =   ρ a τ 550 a r a   S a =   S ρ 0 0 0 S τ 0 0 0 S r   . (16)
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 Å [START_REF] Ångström | On the Atmospheric Transmission of Sun Radiation and on Dust in the Air[END_REF]. 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 [START_REF]QGIS Geographic Information System[END_REF].

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 [START_REF] Eismann | Resolution enhancement of hyperspectral imagery using coincident panchromatic imagery and a stochastic mixing model[END_REF] or FUSE [START_REF] Lanaras | Hyperspectral Super-resolution by Coupled Spectral Unmixing[END_REF], which are based on Bayesian approaches, or methods like Hysure [START_REF] Simoes | A Convex Formulation for Hyperspectral Image Superresolution Via Subspace-based Regularization[END_REF], ICCV'15 [START_REF] Wei | Fast Fusion of Multi-band Images Based on Solving a Sylvester Equation[END_REF], or CNMF [START_REF] Yokoya | Coupled Nonnegative Matrix Factorization Unmixing for Hyperspectral and Multispectral Data Fusion[END_REF], which are based on unmixing analyses. The coupled non-negative matrix factorization (CNMF) algorithm of Yokoya et al. [START_REF] Yokoya | Coupled Nonnegative Matrix Factorization Unmixing for Hyperspectral and Multispectral Data Fusion[END_REF] 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) [START_REF] Nascimento | Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[END_REF] 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) [START_REF] Liu | An Approach Based on Constrained Nonnegative Matrix Factorization to Unmix Hyperspectral Data[END_REF][START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF]. 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. CTMF was developed by Funk et al. [START_REF] Funk | Clustering to Improve Matched Filter Detection of Weak Gas Plumes in Hyperspectral Thermal Imagery[END_REF] 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. [START_REF] Thorpe | Point Source Emissions Mapping Using the Airborne Visible/infrared Imaging Spectrometer (AVIRIS)[END_REF][START_REF] Thorpe | Mapping Methane Concentrations from a Controlled Release Experiment Using the Next Generation Airborne Visible/infrared Imaging Spectrometer (AVIRIS-NG)[END_REF] and Dennison et al. [START_REF] Dennison | High Spatial Resolution Mapping of Elevated Atmospheric Carbon Dioxide Using Airborne Imaging Spectroscopy: Radiative Transfer Modeling and Power Plant Plume Detection[END_REF], in order to detect other gases in the reflective domain. Most recently, Philippets et al. [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF] 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:

q j = C -1 j b j b T j C -1 j b j , (17)
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:

f i,j = γ(τ 550 )q T j ∆ρ p , ( 18
)
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:

τ 550 = γ(τ 550 )τ 550 re f = ||∆ρ p (τ 550 )|| ||∆ρ p simu (τ 550 re f )|| × τ 550 re f . (19)
The RMSE defined by Equation ( 20) was used as an alternative metric to estimate the aerosol type, radius, and AOT.

RMSE = 1 n λ n ∑ λ=λ 1 ∆ρ p (λ, τ 550 ) -γ(τ 550 ).∆ρ p simu (r, λ, τ 550 re f ) 2 , ( 20
)
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.

vis

Data 3.1. 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) [START_REF]Géoservices|Pourquoi #ChoisirGéoportail[END_REF]. 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 [START_REF] Brigot | Adaptation and Evaluation of an Optical Flow Method Applied to Coregistration of Forest Remote Sensing Images[END_REF][START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF].

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 [START_REF] Hagolle | Algorithm Theoretical Basis Document[END_REF]. 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. [START_REF] Hagolle | A Multi-temporal and Multi-spectral Method to Estimate Aerosol Optical Thickness Over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENµs and Sentinel-2 Images[END_REF]. 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.

Results

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).

Image Classification

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) [START_REF] Kruse | The spectral image processing system (SIPS)-interactive visualization and analysis of imaging spectrometer data[END_REF] 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. [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF]. 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 [START_REF] Leoni | Source Impact Determination Using Airborne and Ground Measurements of Industrial Plumes[END_REF][START_REF] Deschamps | Industrial Plume Characterization Using Hyperspectral Imagery[END_REF]. 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. [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF]. However, these values seemed to overestimate in situ measurement over metallurgic plants [START_REF] Leoni | Source Impact Determination Using Airborne and Ground Measurements of Industrial Plumes[END_REF][START_REF] Marris | Fast changes in chemical composition and size distribution of fine particles during the near-field transport of industrial plumes[END_REF]. 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 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.

Discussion

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 [START_REF] Philippets | Anthropogenic Aerosol Emissions Mapping and Characterization by Imaging Spectroscopy-Application to a Metallurgical Industry and a Petrochemical Complex[END_REF][START_REF] Deschamps | Industrial Plume Characterization Using Hyperspectral Imagery[END_REF]. In situ observations over two metallurgic sites in Europe [START_REF] Leoni | Source Impact Determination Using Airborne and Ground Measurements of Industrial Plumes[END_REF][START_REF] Marris | Fast changes in chemical composition and size distribution of fine particles during the near-field transport of industrial plumes[END_REF] 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 [START_REF] Leoni | Source Impact Determination Using Airborne and Ground Measurements of Industrial Plumes[END_REF] and contain metal-bearing nanoparticles [START_REF] Marris | Fast changes in chemical composition and size distribution of fine particles during the near-field transport of industrial plumes[END_REF]. The contribution of fine particles (PM 1 ) can be up to 60% for metallurgy, although the PM 1 fraction depends on the industrial process [START_REF] Ehrlich | PM2.5 and PM1.0-Emissions from industrial plants-Results from measurement programmes in Germany[END_REF]. A similar result has been obtained for a refinery plume using hyperspectral airborne imagery [START_REF] Foucher | Aerosol Plume Characterization from Multitemporal Hyperspectral Analysis[END_REF], 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 [START_REF] Ehrlich | PM2.5 and PM1.0-Emissions from industrial plants-Results from measurement programmes in Germany[END_REF]. 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 [START_REF] Pignatti | The PRISMA hyperspectral mission: Science activities and opportunities for agriculture and land monitoring[END_REF], 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 [START_REF] Leoni | Source Impact Determination Using Airborne and Ground Measurements of Industrial Plumes[END_REF][START_REF] Marris | Fast changes in chemical composition and size distribution of fine particles during the near-field transport of industrial plumes[END_REF].

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.

  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.

Figure 1 .

 1 Figure 1. Normalized (a) extinction and (b) absorption coefficients of MODTRAN background aerosol model.

Figure 2 .

 2 Figure 2. Water pixel: (a) Mean radiance spectrum; and (b) associated signal-to-noise ratio.

Figure 3 .

 3 Figure 3. HYSPEX hyperspectral image (opaque) superposed on the Sentinel-2 image (transparent). Points A and B, respectively, represent areas of interest for the spectra used in Figure 7.

Figure 4 .

 4 Figure 4. Snapshots of Surface reflectances obtained from radiances acquired by HYSPEX camera; (A) surface classification using HYSPEX image alone; and (B) surface classification using a combination of HYSPEX and MSI surface reflectances.

Figure 5 .

 5 Figure 5. Plume mask estimated from (A) HYSPEX image and (B) HYSPEX + MSI. The green dot is the location of the source point.

Figure 6 .

 6 Figure 6. Retrieved plume aerosol parameters: (a) model, (b) AOT, and (c) modal radius, using the CTMF compared with the RMSE.

Figure 7 .

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

Figure 8 .

 8 Figure 8. OEM-retrieved: (a) AOT and modal radius; and (b) radius DOF map.

Figure 9 .

 9 Figure 9. Mean standard deviation of posteriori uncertainties on surface reflectance as a function of wavelength and by surface category.

Table 1 .

 1 Complex aerosol refractive index at 550 nm used in the plume simulations.

	Type	Real Part	Imag. Part
	Brown carbon	1.55	1.2 × 10 -2
	Sulphate	1.52	5 × 10 -4
	Soot	1.83	0.74

Table 3 .

 3 Root-mean-square error (RMSE) and spectral angle mapper (SAM) for surface reflectances outside the plume for HYSPEX image alone and the combination of HYSPEX and MSI images.

	Classes	HYSPEX Alone RMSE SAM	HYSPEX + MSI RMSE SAM
	Water	0.0054	6.68	0.0057	6.10
	Sparse vegetation	0.0426	14.06	0.0304	9.13
	Dense vegetation	0.0414	7.55	0.0393	7.79
	Concrete soils	0.0711	10.94	0.0473	7.60
	Dark soils	0.0337	19.51	0.0267	14.28
	Bright soils	0.1139	9.93	0.1040	7.53

Table 4 .

 4 Mean AOT and radius by class of soil, retrieved with the CTMF and RMSE methods.

	Classes	CTMF	AOT	RMSE	Radius (µm) CTMF RMSE
	Water	0.012		0.011	0.488	0.112
	Dense vegetation	0.028		0.027	0.296	0.268
	Sparse vegetation	0.029		0.026	0.346	0.357
	Concrete soils	0.046		0.039	0.354	0.285
	Dark soils	0.023		0.021	0.217	0.182
	Bright soils	0.095		0.072	0.510	0.286

Table 5 .

 5 Mean radius (µm) by class of soil, retrieved with the OEM methods and from the initial guess inside the final plume mask.

	Classes	OEM	Radius (µm)	Initial Guess
	Water	0.118		0.137
	Dense vegetation	0.151		0.197
	Sparse vegetation	0.231		0.285
	Concrete soils	0.184		0.341
	Dark soils	0.158		0.208
	Bright soils	0.263		0.423

Table 6 .

 6 Mean standard deviations of the posterior uncertainties σ in AOT and modal radius by surface categories. The uncertainty from the measurement σ m and from the resolving error component σ s are associated with this standard deviation.

	Classes	σ	AOT σ m	σ s	σ	Radius (µm) σ m	σ s
	Water	0.003	0.003	0.001	0.049	0.038	0.031
	Dense vegetation	0.004	0.003	0.003	0.064	0.037	0.052
	Sparse vegetation	0.005	0.003	0.004	0.051	0.030	0.042
	Concrete soils	0.007	0.003	0.006	0.048	0.019	0.044
	Dark soils	0.005	0.002	0.004	0.053	0.030	0.041
	Bright soils	0.013	0.006	0.011	0.057	0.022	0.054

Acknowledgments:

The authors would like to thank Yokoya for sharing their code on CNMF. Sentinel-2/MSI data used in this paper are Copernicus data process by CNES for the THEIA thematic center.

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.

Conflicts of Interest:

The authors declare no conflict of interest.