Atmospheric correction over the ocean for hyperspectral radiometers using multi-angle polarimetric retrievals

We developed a fast and accurate polynomial based atmospheric correction (POLYAC) algorithm for hyperspectral radiometric measurements, which parameterizes the atmospheric path radiances using aerosol properties retrieved from co-located multi-wavelength multi-angle polarimeter (MAP) measurements. This algorithm has been applied to co-located spectrometer for planetary exploration (SPEX) airborne and research scanning polarimeter (RSP) measurements, where SPEX airborne was used as a proxy of hyperspectral radiometers, and RSP as the MAP. The hyperspectral remote sensing reflectance obtained from POLYAC is accurate when compared to Aerosol Robotic Network (AERONET), and Visible Infrared Imaging Radiometer Suite (VIIRS) ocean color products. POLYAC provides a robust alternative atmospheric correction algorithm for hyperspectral or multi-spectral radiometric measurements for scenes involving coastal oceans and/or absorbing aerosols, where traditional atmospheric correction algorithms are less reliable. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Oceans are a primary component of the Earth's climate and its ecosystems due to their immense interaction with other Earth components such as the atmosphere and biosphere through cycling, transporting and reserving energy and matter [1]. In order to monitor and evaluate the changes in the Earth system, it is crucial to perform continuous synoptic scale ocean observations, which is only possible with satellite remote sensing. Satellite ocean color remote sensing has permitted the identification and quantification of algal blooms, primary production and ecological variations [2] which further facilitates the understanding on global biogeochemical cycles [3] such as carbon and nitrogen cycles, ocean ecological response to climate change [4] and coastal water quality [5,6].
Ocean color sensors measure the total radiance emanating from the coupled ocean and atmosphere system. Atmospheric correction is performed to remove atmospheric path radiance and ocean surface contributions [7]. This enables estimation of the spectral water leaving signal, which is then used to derive ocean inherent and apparent optical properties that are crucial for ocean and climate studies. The validity of these properties depends on the accuracy of the atmospheric correction scheme that is used to retrieve water leaving signal. A byproduct of employed co-located RSP and SPEX airborne measurements from the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign in 2017 [41,42] to demonstrate the concept by considering SPEX airborne as a hyperspectral radiometer. The Multi-Angular Polarimetric Ocean coLor (MAPOL) joint atmosphere-ocean retrieval algorithm is used to retrieve aerosol properties including optical depth, size distribution, and refractive index from the RSP measurements [28,32,43]. The MAPOL algorithm generates aerosol properties from the discrete bands of RSP, which are then used to parameterize aerosol hyperspectral path radiance in the visible with the POLYAC scheme rooted in the polynomial based algorithm applied to MERIS (POLYMER) [44]. This procedure is very efficient because we avoid using the radiative transfer model to perform hyperspectral simulations for channels covered by hyperspectral radiometers.
We have also used MAPOL to retrieve aerosol and ocean color properties directly from the SPEX airborne polarimetric measurements in selected bands. The retrieval results are compared with the POLYAC algorithm to check the consistency of POLYAC algorithm against direct polarimetric retrievals. Aerosol products from the Aerosol Robotic Network (AERONET) [45,46] and High Spectral Resolution Lidar (HSRL) [47,48] and ocean color products from AERONET-OC [49] and VIIRS are used to validate our retrieval data products.
This paper is organized as follows. Section 2 overviews the ACEPOL field campaign and other datasets used in the study; Section 3 describes the methodology including the MAPOL and POLYAC algorithms; Section 4 presents and discusses the RSP and SPEX airborne retrieval results and the performance of POLYAC; and finally, Section 5 summarizes the conclusions.

Measurements
In this study, we have used polarimeter measurements acquired by RSP and SPEX airborne from the ACEPOL field campaign, which was held from October 19th to November 9th, 2017 (see [41] for download information). To validate our retrieval results, the HSRL-2 lidar measurements from ACEPOL, aerosol data products (level 1.5, version 3.0) from the AERONET USC_SeaPRISM and Monterey sites, and ocean color products from the VIIRS instrument on board the Suomi National Polar-orbiting Partnership (Suomi-NPP) satellite were used.

RSP
RSP is an along track scanner, which includes 152 viewing angles within ±60 • and an IFOV of 14 mrad. It has 9 spectral channels, with central wavelengths of each band located at 410, 470, 550, 670, 865, 960, 1590, 1880 and 2250 nm [50,51]. RSP is equipped with an in flight calibration system [52]. RSP-1 and RSP-2 are the two versions of the RSP instruments, which have different measurement uncertainty characterizations [53]. The RSP measurements over ocean have been used for aerosol and ocean color retrieval in multiple studies [28,29,31,32,[54][55][56] with promising performances. In the ACEPOL campaign the RSP-2 measurements were acquired with a relative radiometric characterization uncertainty of approximately 3% and polarimetric characterization uncertainty of about 0.2%. The instrument noise model for RSP is provided in [53]. In this study, the geolocated level 1B (L1B) measurements of RSP-2 from ACEPOL campaign were used. The RSP channels with wavelengths longer than 865 nm are excluded in the retrieval, as the SPEX airborne instrument does not have spectral sensitivity for longer wavelengths. Viewing angles affected by clouds were discarded. The isotropic Cox-Munk model was used to derive the reflectance and transmittance of the rough ocean surface [57], though sun glint is not included in the retrieval in order to reduce the retrieval uncertainty [28,32].

SPEX airborne
SPEX airborne is a hyperspectral imager developed under the lead of SRON Netherlands Institute for Space Research [52,58]. It covers the spectral range from 380 nm to 800, but the best performance is for 450-750 nm. SPEX airborne has a 10-25 nm resolution for degree of linear polarization (DOLP, the ratio of linearly polarized to total radiance) and 2 nm resolution for radiance. It performs multi-angular measurements at 9 fixed viewing angles: ±56 • , ±42 • , ±28 • , ±14 • and 0 • . It has a 4% accuracy in spectral radiance and the DOLP accuracy is 0.5% [59,60]. The SPEX airborne and RSP measurements were compared at 410, 470, 550, and 670 nm, which showed a systematic difference in radiance and a random difference in DOLP. The root mean square (RMS) difference for ocean scenes in DOLP at 410, 470, 555, and 670 nm were 0.008,0.006, 0.004, and 0.008, respectively, and the RMS of the percentage difference in radiance at 410, 470, 555 and 670 nm were 4%, 3%, 2%, and 3%, respectively [59]. For this reason, SPEX airborne polarimetric retrieval studies [26,27] have disregarded the wavelengths shorter than 450 nm due to the large discrepancy with the RSP measurements. They have also excluded wavelengths longer than 750 nm due to the channel cross contamination from the grating order overlap.
The spectral behavior of total radiance measurement is crucial for retrieving water leaving radiance accurately [61]. Therefore, the exclusion of wavelengths smaller than 450 nm or inclusion of defective measurements in the polarimetric retrieval can result in erroneous or highly uncertain retrievals. Despite previous studies discarding the shortest wavelengths, in this study we have used 5 discrete wavelength bands: 398, 410, 470, 555, and 670 nm to perform the SPEX airborne polarimetric retrieval after application of additional corrections (described in detail below). Viewing angles affected by clouds and sun glint are excluded in the retrieval. Basic information regarding the measurements for all the retrieval cases are summarized in Table 1 (see below) whereas the locations and scanning directions of each case is summarized in Fig. 1. The comparison between RSP and co-located SPEX airborne measurements [59] shows systematic and random errors in radiance and DOLP values respectively. Figure 2 shows the comparison of radiance measurements for the selected cases in this study, which shows small yet non-negligible differences, especially at 410 nm. The radiometric calibration difference between RSP and SPEX airborne has led to large difference in water leaving retrieval in SPEX airborne MAPOL retrieval (Sec. 4.2). Thus, we have rescaled/corrected the SPEX airborne radiometric measurements using a hyperspectral error function. The hyperspectral error function is based on the average of the relative differences between co-located SPEX L1C ( L1B measurements from different viewports aggregated on a spatial grid) and geolocated RSP L1B radiometric measurements in 410, 470, 555 and 670 nm wavelength bands for the nadir viewing angle. The hyperspectral error function was obtained using a second order polynomial fit with the errors corresponding to the above four wavelength bands (Fig. 2). In the SPEX airborne polarimetric retrieval, only the radiometric measurements in 410, 470, 555, and 670 nm wavelength bands are rescaled with respect to the error function to conform the co-located RSP measurements. The SPEX airborne DOLP data were not rescaled because the errors in DOLP measurements are random. The spectral correction was not applied to the 398 nm channel because the potential large errors associated with extrapolation (shortest wavelength of RSP is 410 nm). Our sensitivity study shown in Section 4.1 and 4.2 demonstrates that the inclusion of the defective wavelengths after spectral correction can improve the polarimetric retrieval performance significantly at these discrete wavelengths. We compared the retrieval of remote sensing reflectance with and without the SPEX airborne measurements at 398 nm and found that its inclusion improves the retrieval performance slightly. As such, we chose to include the 398 nm channel, as it increases the information content of the aerosol properties, which potentially outweighs its calibration error.
In the POLYAC retrieval, the retrieved hyperspectral water leaving signal below 450 nm still shows some calibration features due to the imperfect correction. Therefore, we excluded the signals below 450 nm in the POLYAC results. This won't be an issue for applying the POLYAC to the PACE mission as the calibration of the spectrometer will be much improved.
The SPEX airborne and RSP measurements that are coincident in temporal and spatial domains were used to implement the POLYAC scheme. Co-located pixels were selected based on geolocation overlapping within a maximum distance of 300 m. The instruments mounted in opposite airplane wings may have orientation shift due to wing flex, which was not considered in co-location. The different ground pixel sizes in different angles are also ignored. The nadir ground pixel of RSP in the ACEPOL campaign is around 280 m whereas the spatial sampling of the SPEX airborne L1C data is 1 km x 1km (across x along track).

VIIRS, AERONET, and HSRL-2
VIIRS is a single along track view angle multi-spectral imager. The VIIRS data were processed by the standard NASA atmospheric correction algorithm [7,62]. The Level 2 VIIRS data provides ocean color products such as remote sensing reflectances (R rs ) at 5 visible wavelengths (410, 443, 486, 551, and 671 nm) and [chla] (OCI algorithm; [63]). The USC_SeaPRISM AERONET site has aerosol optical depth (AOD) data available for 10/23/2017 and 10/25/2017 and Monterey AERONET site has AOD data available for 11/07/2017. The USC_SeaPRISM has normalized water leaving radiance ([L w ] ex N ) data available for 10/23/2017 and 10/25/2017, which can be easily rescaled to Rrs. HSRL-2 provides AOD at 532 nm for all the cases. During the ACEPOL campaign HSRL-2 experienced an interference that appeared to be related to atmospheric turbulence. This impacted the ability to use the 532 nm and 355 nm molecular channels to derive aerosol extinction with the HSRL technique. In such cases, the extinction profiles were estimated by multiplying the aerosol backscatter profiles with an assumed aerosol lidar ratio. For low AOD cases, such as cases investigated in this study, the AOD product derived using an assumed lidar ratio is expected to have lower uncertainties. RSP and SPEX airborne measurements were co-located with the HSRL-2 and VIIRS data within a maximum distance of 500 m. For both the VIIRS and AERONET products, measurements within 1 hour from the RSP scanning time were selected.

MAPOL algorithm
MAPOL is a joint retrieval algorithm that utilizes RSP measurements to retrieve both aerosol and ocean properties. It has been validated with both synthetic and in-situ measurements [28,32,43]. We have further extended the algorithm to employ the SPEX airborne measurements in selected wavelength bands. The algorithm is based on a Levenberg -Marquardt non-linear least squares optimization [64], which minimizes the following cost function; where ρ t = πL t /µ 0 F 0 is the total reflectance; ρ p = π √︂ Q 2 t + U 2 t /µ 0 F 0 is the polarized reflectance; and σ t and σ p are the total uncertainty of total and polarized reflectance measurements, respectively. The superscript m denotes the measurement; f denotes the forward model simulation; x is the state vector of the retrieval; i is the measurement index corresponding to a particular angle or wavelength of the polarimeter; and N is the total number of measurements used in the retrieval. L t , Q t , and U t are the first three Stokes parameters measured at sensor level; µ 0 is the cosine of the solar zenith angle, and F 0 is the extraterrestrial solar irradiance corrected for the Sun-Earth distance. Total uncertainty is a combination of uncertainties due to measurement noise [53], spatial averaging, and forward modeling.
The forward model in MAPOL [43] is the vector radiative transfer model for coupled oceanatmosphere systems based on the successive order for scattering method [65]. We assume the coupled ocean and atmosphere system as a three-layer system: an upper molecular layer, a middle layer of mixed molecules and aerosols, and a lower ocean layer topped by a rough ocean surface. The aerosol size distribution is represented by six aerosol sub-modes; three fine mode and three coarse mode aerosols, each with a log-normal distribution. The complex refractive index spectra of aerosols is represented by two principle components obtained from a principle component analysis (PCA) performed on aerosol refractive index measurements of fine and coarse mode aerosols [66,67]. MAPOL retrieves the refractive index spectra of aerosols by using two PCA components for both the fine and coarse modes, and real and imaginary parts, which includes a total of 8 parameters (2 (fine and coarse) modes × 2 PCA × 2 parts (real and imaginary) of refractive index).
MAPOL includes two ocean bio-optical models to treat open and coastal waters [11]) separately. These models represent the inherent optical properties (IOPs) of phytoplankton, CDOM and sediments/non-algal particles (NAP). The bio-optical model of open ocean waters solely depends on [chla] whereas the generalized coastal model is characterized by seven bio-optical parameters [43].
In this study, the RSP and SPEX airborne polarimetric retrievals were performed with the [chla] based open water bio-optical model, as the available test cases are more likely to associate with clear waters. The retrieval parameter space of this study therefore includes a total of 16 parameters: 8 parameters to represent PCA based aerosol refractive index spectra, 6 parameters to represent aerosol volume distribution, 2 more parameters on ocean: wind speed and [chla]. In order to ensure that the global minimum is achieved in the retrieval, MAPOL was run several times with different initial retrieval parameter values (Sec. 4.). (For further information regarding the algorithm, readers are referred to [43]) MAPOL is computationally demanding as it needs to iteratively run the radiative transfer forward model for coupled atmosphere and ocean systems. Typically, it takes several hours for a single core CPU to process one pixel retrieval acquired by RSP, which is too slow for processing satellite measurements operationally. We are developing a new version of MAPOL, which replaces the radiative transfer forward model with a neural network algorithm. This new version, called FastMAPOL [68], can process one RSP pixel within a few seconds in CPU, and with retrieval quality equivalent to MAPOL. Currently, FastMAPOL uses [chla] to parameterize ocean water optical properties. To extend FastMAPOL to cover coastal water, we will integrate a recently developed ocean reflectance model for coastal waters based on neural network [69]. FastMAPOL is the subject of a different manuscript in preparation, which we will not discuss in detail in this paper.

Hyperspectral atmospheric correction
Heritage retrieval algorithms for spectrometers with single viewing angle employ observations in the NIR to determine aerosol properties, which are extrapolated to visible wavelengths for atmospheric correction purposes. For pixels contaminated by sun glint, the sun glint signal is too large that the ocean signal becomes too small to be reliably extracted from the total measurements. To improve this issue, Steinmetz et al. proposed the POLYMER atmospheric correction algorithm, which is based on spectral optimization [44]. The POLYMER algorithm utilizes the entire spectral range of the measurements and is applicable over the whole sun glint pattern. It includes a water reflectance model targeting coastal waters and a polynomial (Eq. (8)) to model the residual sun glint and spectral reflectance from the atmosphere. This approach has been validated with MEdium Resolution Imaging Spectrometer (MERIS), MODIS, Sea-viewing Wide Field-of-view Sensor (SeaWiFS), VIIRS and Ocean and Land Colour Instrument (OLCI) measurements [70].

Components of the sensor measurements
The total reflectance ρ sensor t (λ) measured by a sensor can be expressed as [71]; where t gas is the transmittance due to atmospheric gas absorption, ρ atm+sfc is the contribution from the atmosphere and surface; the subscripts, r, a, g, wc, and c denote reflectance ρ at sensor level contributed from Rayleigh scattering, aerosol scattering, direct sun glint, white caps, and coupling between molecules, aerosols, and sun glint. ρ + w is the water leaving signal just above the water, t u is the upward radiance transmittance from water surface to sensor, and t d is the total downward irradiance transmittance from TOA to water surface for atmospheric scattering. Atmospheric correction is a process to obtain ρ + w (λ) from ρ sensor t . The water leaving reflectance at the sensor level (ρ sensor w ) can be written as; where L + w is the water leaving radiance just above the water surface. Remote sensing reflectance defined as R rs = (L + w )/(F + d ) is commonly used in ocean color community to represent water leaving signal where F + d is the downwelling irradiance just above the water surface. Different instruments (RSP, SPEX airborne, AERONET, and VIIRS) have different viewing and solar geometries which hinders the ability to inter-compare the water leaving signals. In order to compare the results from different instruments, water leaving radiance is converted to exact normalized water leaving radiance [L w ] ex N using the BRDF correction. [L w ] ex N corresponds to a hypothetical radiance that would be measured at nadir if the Sun were at zenith and Earth at its mean distance from the Sun in the absence of atmosphere. In the following we present R rs = [L w ] ex N /F 0 . A detailed mathematical overview of our BRDF correction is given in the appendix.

POLYMER
The POLYMER algorithm initially corrects for ozone gas absorption, Rayleigh scattering, sun glint, and coupling between Rayleigh and sun glint using radiative transfer simulations with the isotropic Cox-Munk model for ocean surface. The accuracy of this correction depends on the wind speed adopted from European Centre of Medium-Range Weather Forecasts (ECMWF). When the wind speed estimation has a large error, the initial correction becomes ineffective to fully correct the sun glint. The remaining part of the total reflectance (ρ ′ ) after the initial correction becomes, where t oz is the transmittance of ozone and ∆g denotes the residue of the sun glint due to inaccurate estimation of wind speed and the fact that aerosol extinction was not considered in the initial estimation of sun glint using radiative transfer simulation. ρ ag is then modelled as; where c 0 , c 1 , and c 2 are the fitting coefficients to be estimated by least squares fitting of the spectral measurements and the transmittance T 0 (λ) can be expressed as [44], where Rayleigh optical thickness τ r ≃ 0.00877λ −4.05 with λ in nm, and µ s and µ v are the cosines of solar zenith and viewing zenith angles, respectively. For ρ g >0.02, T 0 becomes the direct transmission and for ρ g <0.02, T 0 becomes the diffuse transmittance. The aerosol transmittance is assumed to be negligible in Eq. (9). The first term c 0 in Eq. (8) represents the spectrally flat components of the reflectance (coarse mode aerosols, clouds, etc.); the second term models the aerosol contribution with the Angstrom exponent assumed as 1; and the last term stands for the coupling between spectrally flat components and Rayleigh scattering. For more information regarding the POLYMER algorithm including its water reflectance model and minimization procedures, readers are referred to [44].

POLYAC
Improved atmospheric correction was a primary motivation for including two polarimeters, HARP2 and SPEXone, in the PACE mission [35,36]. The POLYAC algorithm (Fig. 3) takes advantage of the rich information content of the MAP data to assist atmospheric correction of co-located hyperspectral radiometer measurements to estimate R rs . POLYAC adopts the POLYMER polynomial (Eq. (8)) to parametrize atmospheric path radiance ρ atm+sfc . It utilizes the retrieved aerosol and ocean color information from MAPs with the MAPOL algorithm [28,32,43], which provides wind speed, aerosol size distribution and volume concentration, and spectral complex refractive indices at discrete wavelengths of MAP. The retrieved wind speed and aerosol information is used to calculate ρ atm+sfc (λ) with our radiative transfer model [65] where the atmospheric gas absorption is considered. Due to this, ρ atm+sfc (λ) is later corrected for atmospheric gas absorption based on the Beer-Lambert law, which is then used to fit the gas absorption corrected hyperspectral atmospheric path reflectance ρ ′ atm+sfc (λ) with the following polynomial adapted from Eq. (8): where the aerosol Angstrom exponent is relaxed as a general value m here. The fitting parameters c ′ 0 , c ′ 1 , c ′ 2 , and m are determined through nonlinear least squares fitting of spectral ρ atm+sfc . c ′ 0 represents all the spectrally flat components, whereas the last polynomial term represents Rayleigh scattering and associated coupling terms. The PACE OCI is a tilting sensor designed to minimize sun glint [36] so that we have neglected sun glint contribution in Eq. (7). The atmospheric gas absorption is based on the vertical molecule density profiles specified by the 1976 US standard atmosphere [72] for ozone, oxygen, water vapor, nitrogen dioxide, methane, and carbon dioxide.The spectral ρ ′ atm+sfc (Eq. (10)) can be used to model the hyperspectral atmospheric path reflectance at any given wavelengths. In order to estimate the accuracy and feasibility of Eq. (10), we have applied it to the co-located RSP and SPEX airborne measurements (rescaled with respect to the error function) by considering SPEX airborne as a hyperspectral radiometer operating in the wavelength range of 410 -670 nm, along with RSP aerosol retrievals based on MAPOL. The comparison of ρ atm+sfc from RSP polarimetric retrieval with the parameter fitting of ρ ′ atm+sfc in Eq. (10) is shown in Fig. 4(a). The percentage error between RSP and the POLYAC scheme for each RSP wavelength is shown in Fig. 4(b) which indicates agreement less than measurement uncertainty at the four selected RSP wavelengths. The hyperspectral water leaving reflectance for SPEX airborne measurements is estimated from Eq. (4). MAPOL retrieval also provides the total upward radiance transmittance and total downward irradiance transmittance, from which the respective transmittances due to gas absorption are removed to obtain t d and t u . They are then used to fit the hyperspectral transmittance using the following equations: The fitting parameters c ′ 3 and c ′ 4 are obtained through nonlinear least square fitting of spectral t d and t u . F + d = F 0 µ 0 t d with F 0 provided by the hyperspectral irradiance measurements of solar irradiance. [73].

Results
Polarimetric retrievals were carried out with MAPOL over nine different sites for RSP and six for SPEX airborne, respectively. Among these, five co-located cases of RSP and SPEX airborne were found and used to evaluate the POLYAC algorithm, which is summarized in Table 1. For each case MAPOL was run 40 times with different initial retrieval parameter values to ensure the global minimum is achieved. The 40 retrieval results were sorted according to the χ 2 values and the cumulative probability (CP) is calculated for each χ 2 . The minimum and maximum values of χ 2 are denoted as χ 2 min and χ 2 max , respectively. Note that CP( χ 2 min )=0 and CP( χ 2 max )=1. Furthermore we introduce a χ 2 cut-off value χ 2 cut such that CP( χ 2 cut )=70%. The final retrieval products were then obtained by averaging all the retrieval outputs of χ 2 with χ 2 min < χ 2 < χ 2 cut . For the cases studied in this paper, χ 2 min is in the range of 0.563 -1.508, and χ 2 cut within 1.602 -2.06, respectively. We choose CP( χ 2 cut )=70% because CP=70% corresponds approximately to the 1σ uncertainty [28,32,43]. For the SPEX airborne retrievals, χ 2 min ranges within 0.91 -1.24 and χ 2 cut within 1.30 -1.34.

Polarimetric retrievals of AOD
There are no AERONET AOD measurements available for comparison with all the cases in this study, thus the AOD values obtained from the RSP and SPEX airborne polarimetric retrievals were compared with the HSRL-2 532 nm AOD. There are two types of plots added for comparison, scatter plot ( Fig. 5(a)) and Bland-Altman difference plot [74] (Fig. 5(b)). The overall comparison of 532 nm AOD retrieved from RSP/SPEX airborne is satisfactory with an RMSE of 0.0041/0.0067, Pearson correlation coefficient of 0.908/0.883 and Spearman's rank coefficient of 0.849/0.943. The differences between the HSRL-2 and RSP/SPEX airborne AOD can be explained by different sampling volumes and viewing geometries of the instruments, low aerosol loadings, information content associated with different viewing angles and wavelength range and contamination by nearby clouds.

Polarimetric retrievals of R rs
The spectral R rs obtained from the polarimetric retrievals were compared with the data product from the AERONET USC_SeaPRISM site, which was 1.5 km and 1 km away from the case 1 and case 3 locations, respectively. No co-located AERONET-OC measurements are available for the other three cases. Hence the co-located VIIRS products were used. The AERONET measurements averaged within one hour from the RSP scanning time were considered. R rs estimated from RSP and SPEX airborne show similar spectral shapes (Fig. 6) except for 410 nm, which can be explained by the retrieval uncertainty associated with larger reflectance values.The differences between R rs obtained for RSP and other data products can be attributed to the time and footprint differences, small magnitudes of water leaving signals, the differences in the VIIRS and MAPOL atmospheric correction algorithms and less information content of the VIIRS measurements. For SPEX airborne, partially corrected measurements might also have affected the retrieval accuracy as the water leaving radiance is highly sensitive to the measured total spectral radiance especially within shorter wavelengths where Rayleigh scattering is dominant. Fig. 6. The comparison of RSP and SPEX airborne estimated R rs (Sr −1 ). Vertical bars indicate the retrieval uncertainties of RSP and SPEX airborne retrievals and measurement uncertainties in AERONET measurements. SPEX corrected retrieval stands for R rs obtained from polarimetric retrieval with corrected measurements. Time is given in UTC.
Overall, R rs and AOD products obtained from RSP and SPEX airborne instruments are in good agreement with each other and with AERONET and VIIRS data products. This indicates the validity of MAPOL algorithm for further use in the POLYAC atmospheric correction algorithm.

POLYAC retrieval of hyperspectral R rs
The hyperspectral R rs for SPEX airborne within 470 -675 nm was estimated using the POLYAC scheme (Fig. 7). Extrapolation of atmospheric correction in the UV spectral range (380 nm -410 nm) was not considered as the MAPOL retrieval with the RSP wavelengths (>410 nm) cannot provide accurate aerosol properties in the UV spectral range [35]. Wavelengths less than 450 nm were also disregarded due to the larger discrepancy between the SPEX airborne and RSP measurements with respect to the hyperspectral error function.
As we have described earlier, rescaled SPEX airborne measurements based on the hyperspectral error function are used to estimate hyperspectral R rs . The R rs obtained through POLYAC for the corrected SPEX airborne measurements agrees well with that from the polarimetric RSP and SPEX airborne retrievals. Fig. 7. Comparison of R rs (Sr −1 ) obtained from the RSP and SPEX airborne retrievals with hyper spectral R rs obtained with POLYAC scheme. The dashed line represents R rs obtained from POLYAC scheme for rescaled SPEX airborne measurements (SPEX), while the solid lines are R rs obtained from the POLYAC scheme for rescaled SPEX measurements. Stars indicate R rs obtained from SPEX airborne retrievals with rescaled measurements. R rs obtained from RSP retrievals are represented in filled circles In the retrieved R rs spectrum, there are some gas absorption features which are not completely removed, particularly noticeable at the O 2 γ band around 629 nm. Accurate gas absorption correction is challenging due to the interaction between multiple scattering and gas absorption and lack of information on volume density profiles of water vapor and ozone, aerosol height and exact instrument line shape function of the instrument.
The comparison between R rs obtained for POLYAC scheme and RSP retrieval under different wavelengths is shown in Fig. 8. The RMSE is less than 0.00012 for all the three wavelengths considered in the comparison.
To understand the spatial variability of R rs , we applied POLYAC to the SPEX airborne measurements over a region surrounding the case 1 retrieval site in Fig. 7. Figure 9(a) shows the locations of the RSP and SPEX airborne measurements used in this work. Co-located VIIRS data are used for comparison. We assumed the aerosol properties are uniform over this small region [75] and applied the atmospheric path radiance and transmittance estimated from RSP MAPOL retrieval for case 1 to the retrieval of R rs . Due to the non-negligible differences between RSP and SPEX airborne radiance measurements, we rescaled all the SPEX airborne radiance measurements across the region using the same scaling factor we used for case 1 retrieval site.
The comparison of R rs obtained from POLYAC for SPEX airborne measurements with VIIRS across the region is given in Fig. 9(b). The R rs spectra from the AERONET-OC and MAPOL retrieval for the case 1 retrieval site are also shown. Within the selected region the mean value of R rs at 550 nm for SPEX airborne is 0.0015 Sr −1 , slightly higher than that of the VIIRS mean value 0.0013 Sr −1 , which is consistent with the observed difference between AERONET and VIIRS. At 670 nm the mean values are 8.5×10 −5 and 6.8×10 −5 , respectively, for SPEX airborne and VIIRS. The standard deviation of the R rs spectra for SPEX airborne and VIIRS at 550 nm are 1.2×10 −4 and 7.3×10 −5 , respectively. At 670 the standard deviations are 1.6×10 −4 and 3.9×10 −5 , respectively, for SPEX airborne and VIIRS. At both wavelengths the variabilities of VIIRS product is smaller than those of SPEX airborne data within the same region, mainly due to the smaller range of viewing angles among these pixels, as VIIRS is space-borne and far away from the surface.
The results indicate that POLYAC successfully applies the aerosol and surface properties obtained from the MAP retrievals to assist atmospheric correction of hyperspectral radiometers. Challenges remain in the consistency of the radiometric calibration among different instruments and gas absorption correction. Given the rigorous calibration requirements assigned to all instruments onboard the PACE observatory, the POLYAC scheme is expected to result in more accurate performance when considering OCI and the co-located HARP2 and SPEXone datasets.   9. (a) Location map: SPEX airborne measurements (filled cyan circles). Region selected to apply POLYAC (red crosses). RSP track (green dots). RSP measurements used in the MAPOL retrieval/ case 1 retrieval site (blue star). VIIRS data selected for comparison (filled gray circles). (b) The comparison of R rs obtained from POLYAC for SPEX airborne, and VIIRS across the selected region

Conclusions
We have retrieved aerosol and ocean properties using RSP and SPEX airborne for five cases. The MAPOL retrieval algorithm has been extended to be applied to the SPEX airborne measurements and its performance is validated. The retrieved AOD was compared with AOD from HSRL-2 instrument and R rs was compared with the AERONET and VIIRS products. The comparison of retrieved products shows a good agreement for most of the cases, which demonstrates satisfactory performance of the MAPOL algorithm.
We have also presented the POLYAC scheme to perform the hyperspectral atmospheric correction using the MAP retrievals. The performance and the accuracy of POLYAC scheme have also been validated using the RSP and co-located SPEX airborne measurements. The comparison between RSP and SPEX airborne retrieved R rs and POLYAC estimated R rs for corrected SPEX airborne measurements shows a good agreement. We have also applied POLYAC over a small region and the retrieved R rs was compared with that from VIIRS. We demonstrated that it is important to have highly accurate and consistent radiometric calibrations between MAP and spectroradiometers in order for POLYAC to retrieve R rs well. Future development of the POLYAC scheme requires evaluation of existing assumptions and accurate handling of hyperspectral gas absorption correction. The MAPOL algorithm together with POLYAC scheme can therefore be used as a computationally effective and robust atmospheric correction scheme for hyperspectral radiometers such as OCI, which will be in synergy with the two MAPs (SPEXone and HARP2) on board PACE mission.

Appendix: mathematical representation of R rs
The downward direct irradiance transmission is formulated by, where F + d is the downwelling irradiance just above the water surface calculated from the forward model. Transmittance of the upwelling water leaving radiance which depends on the viewing direction is given by, where the superscript sensor and + denote the sensor level and just above the water, respectively.
(atm + sfc only) represents cases where no water column contribution is considered except the water surface, i.e. corresponding to atmospheric correction. Normalized water leaving radiance can be written as [76], where r 0 refers to the mean distance between sun and earth. The BRDF correction is used to convert [L w ] N into exact normalized water leaving radiance [L w ] ex N [76].
where R is a dimensionless factor which depends on viewing direction and wind speed, and accounts for reflection and refraction effects when light propagates through wavy interfaces. For viewing angles <15 • , R=0.53. (In this study we used nadir as the viewing angle, hence R=R 0 =0.53). Q is the bidirectional function and f is a dimensionless coefficient that accounts for irradiance that varies with water optical properties and illumination conditions [76].
where F − u and F − d are the upward and downward irradiances just below the water surface, respectively; L − u is the upward radiance just below the water surface; a and b b are the total absorption and back scattering coefficients of the ocean body, respectively. f 0 is the factor f for sun at the zenith and Q 0 is Q for sun at zenith and viewed at nadir.
For nadir viewing, R rs becomes; For hyperspectral [L w ] ex N calculations, the hyperspectral BRDF correction was obtained using a third order polynomial fit with the known BRDF correction factors retrieved from the MAPOL retrieval.