Efficient multi-angle polarimetric inversion of aerosols and ocean color powered by a deep neural network forward model

NASA’s Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission, scheduled for launch in the timeframe of 2023, will carry a hyperspectral scanning radiometer named the Ocean Color Instrument (OCI) and two Multi-Angle Polarimeters (MAP): the UMBC Hyper-Angular Rainbow Polarimeter (HARP2) and the SRON Spectro-Polarimeter for Planetary EXploration one (SPEXone). The MAP measurements contain rich information on the microphysical properties of aerosols and hydrosols, and therefore can be used to retrieve accurate aerosol properties for complex atmosphere and ocean systems. Most 5 polarimetric aerosol retrieval algorithms utilize vector radiative transfer models iteratively in an optimization approach, which leads to high computational costs that limit their usage in the operational processing of large data volumes acquired by the MAP imagers. In this work, we propose a deep neural network (NN) forward model to represent the radiative transfer simulation of coupled atmosphere and ocean systems, for applications to the HARP2 instrument and its predecessors. Through the evaluation of synthetic datasets for AirHARP (airborne version of HARP2), the NN model achieves a numerical accuracy 10 smaller than the instrument uncertainties, with a running time of 0.01s in a single CPU core or 1 ms in GPU. Using the NN as a forward model, we built an efficient joint aerosol and ocean color retrieval algorithm called FastMAPOL, evolved from the well-validated Multi-Angular Polarimetric Ocean coLor (MAPOL) algorithm. Retrievals of aerosol properties and water leaving signals were conducted on both the synthetic data and the AirHARP field measurements from the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign in 2017. From the validation with the synthetic data and the collocated 15 High Spectral Resolution Lidar (HSRL) aerosol products, we demonstrated that the aerosol microphysical properties and water leaving signals can be retrieved efficiently and within acceptable error. Comparing to the retrieval speed using conventional radiative transfer forward model, the computational acceleration is 10 times faster with CPU or 10 times with GPU processors. The FastMAPOL algorithm can be used to operationally process the large volume of polarimetric data acquired by PACE and other future Earth observing satellite missions with similar capabilities. 20

instrument uncertainties, with a running time of 0.01s in a single CPU core or 1 ms in GPU. Using the NN as a forward model, we built an efficient joint aerosol and ocean color retrieval algorithm called FastMAPOL, evolved from the well-validated Multi-Angular Polarimetric Ocean coLor (MAPOL) algorithm. Retrievals of aerosol properties and water leaving signals were conducted on both the synthetic data and the AirHARP field measurements from the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign in 2017. From the validation with the synthetic data and the collocated High Spectral Resolu- 15 tion Lidar (HSRL) aerosol products, we demonstrated that the aerosol microphysical properties and water leaving signals can be retrieved efficiently and within acceptable error. The FastMAPOL algorithm can be used to operationally process the large volume of polarimetric data acquired by PACE and other future Earth observing satellite missions with similar capabilities.
where the 670 nm band will measure 60 viewing angles and the other bands 10 viewing angles, with a swath of 1,556 km at nadir on the Earth surface. To facilitate cross calibrations and validations, a PACE Level-1C common data format has been developed, with the purpose of projecting all three PACE instruments onto an uniform spatial grid (Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission, 2020). The PACE instruments will provide an unprecedented opportunity to improve the characterization of the atmosphere and ocean states (Remer et al., 2019a, b;Frouin et al., 2019). 5 To retrieve the aerosol information from polarimetric measurements over oceans, several advanced aerosol retrieval algorithms have been developed for both airborne and spaceborne MAPs, such as POLDER/PARASOL (Hasekamp et al., 2011;Dubovik et al., 2011Dubovik et al., , 2014Li et al., 2019;Hasekamp et al., 2019b;Chen et al., 2020), the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) (Xu et al., 2016, SPEX Airborne (the airborne version of SPEXone ) (Fu and Hasekamp, 2018;Fu et al., 2020;Fan et al., 2019), the Research Scanning Polarimeter (RSP) (Chowdhary et al., 2005;Wu et al., 2015;10 Stamnes et al., 2018;Gao et al., 2018Gao et al., , 2019Gao et al., , 2020, the Directional Polarimetric Camera (DPC)/GaoFen-5 (Wang et al., 2014;Li et al., 2018). The retrieval algorithms are mostly based on iterative optimization approaches that utilize vector radiative transfer (RT) models as the forward model. The high computational cost of the RT simulations pose great challenges in the operational processing of the large data volumes acquired by the MAP imagers. To alleviate this issue, the SPEX team represented the polarimetric reflectance for an open ocean system using a deep neural network (NN) and coupled it with a radiative 15 transfer model for the atmosphere (Fan et al., 2019). This hybrid forward model avoids the direct calculation of the scattering and absorption properties inside the ocean, and still maintains high accuracy, therefore enabling sufficient efficiency for SPEXone data retrieval. For coastal waters, Mukherjee et al. (2020) developed a NN model to predict the polarimetric reflectance associated with complex water optical properties. This NN model can be combined with a flexible atmosphere model for MAP aerosol retrievals over complex waters. 20 For non-polarimetric remote sensing studies, several NN approaches have been developed to derive aerosol and ocean properties simultaneously (Fan et al. (2017); Shi et al. (2020) and reference within). Fan et al. (2017) developed NN models to directly invert the aerosol optical depth (AOD) and remote sensing reflectance R rs (λ) (sr −1 ) from the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) measurements. Shi et al. (2020) developed a NN radiative transfer scheme for coupled atmosphere and ocean systems including both open and coastal waters, which is then applied in an optimal estimation 25 algorithm for the Cloud and Aerosol Imager-2 (CAI-2) hosted on the Greenhouse gases Observing Satellite-2 (GOSAT-2).
A number of NN models have been developed to directly invert the aerosol microphysical properties from MAP measurements. Di Noia et al. (2015) discusses the NN employed to retrieve aerosol refractive index, size, and optical depth (AOD) from GroundSPEX (a ground version of SPEX instrument) measurements. Di Noia et al. (2017) developed a NN inversion method for airborne MAP measurement over land from RSP. In both works, the results from the NN inversion are further used as initial 30 values for iterative optimization, and both efficiency and the retrieval accuracy are shown to be improved. Using NN to conduct direct inversion is efficient, but it is often viewed as a black box and it is difficult to account for measurement uncertainties.
The combination of a NN inversion with an iterative optimization method shows promise for MAP retrievals.
Even with such ample progress, it is still challenging for current state-of-the-art algorithms to process MAP data operationally through iterative optimization. In this work, we present a joint retrieval algorithm for aerosol properties and water 35 leaving signals that uses a deep NN model to replace the radiative transfer forward model for simulation of the polarimetric reflectances. This approach is one step further than (Fan et al., 2019), as both the atmospheric and oceanic radiative transfer processes are represented by the NN. The NN forward model is then used in an iterative retrieval algorithm that is significantly more computationally efficient than approaches that use traditional radiative transfer. The benefits using a NN model as the forward model in retrieval algorithms can be summarized as follows with details provided in later sections: 5 -Fast: NN models mostly involve matrix operations that can be evaluated efficiently.
-Accurate: Given sufficient training data volumes and accuracies, NN models can be trained with high precision.
-Differentiable: The Jacobian matrix of NN models can be represented analytically and therefore further improves efficiency and accuracy in retrievals.
-Transferable: The parameters of a NN can be exported and implemented into existing retrieval algorithms. 10 The retrieval algorithm we developed is called FastMAPOL, which is evolved from the well-validated Multi-Angular Polarimetric Ocean coLor (MAPOL) algorithm (Gao et al., 2018(Gao et al., , 2019 by replacing its forward model with NN models. To validate the retrieval algorithm, we applied FastMAPOL to both synthetic and field measurements from AirHARP (the airborne version of HARP2 and HARP CubeSat) for the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign in 2017 . The synthetic AirHARP data is a supplement of the field measurements with a wider range 15 of aerosol and ocean optical properties, and solar and viewing geometries. The AOD derived from coincident High Spectral Resolution Lidar (HSRL, Hair et al. (2008)) and Aerosol Robotic Network (AERONET, Holben et al. (1998)) measurements are used to evaluate the performance of the AOD retrieval from the AirHARP field measurements. Using the retrieved aerosol properties, atmospheric correction is applied to the AirHARP measurements to derive the water leaving signal at four AirHARP bands. The retrieved aerosol products from MAP can also assist hyperspectral atmospheric correction on instruments such as 20 PACE OCI as previously demonstrated using the aerosol properties retrieved from RSP and hyperspectral measurements from SPEX Airborne Hannadige et al., 2020). Retrieval uncertainties of both aerosol and water leaving signals under various aerosol loadings are also discussed in this study. The retrieval algorithm powered by the NN forward model provides a practical approach for operational applications of polarimetric aerosol and ocean color retrieval for PACE, and other satellite missions that utilize polarimeters in the retrieval of geophysical properties from Earth observations. 25 The paper is organized into seven sections:, Sect. 2 reviews the retrieval algorithm and its radiative transfer forward model, Sect. 3 discusses the training and accuracy of the NN forward model, Sec 4. applies the NN forward model to aerosol and water leaving signal retrievals from the synthetic AirHARP data, Sect. 5. discusses the retrievals on AirHARP field measurements from ACEPOL campaign, Sect. 6 and 7 provide discussions and conclusions.
2 Joint aerosol and ocean color retrieval algorithm 30 In this section, we will discuss the MAPOL retrieval algorithm based on multi-angle polarimetric measurements and the associated radiative transfer forward model. The retrieval algorithm have been validated using both synthetic data (Gao et al., will first discuss the AirHARP instrument characteristics. AirHARP measures the total and linearly polarized radiance at 60 viewing angles at the 660 nm band, and at 20 viewing angles at the 440, 550, and 870 nm bands. Different from AirHARP, HARP2 reduces the number of viewing angles to 10 at 440, 550, and 870 nm, and maintains 20 viewing angles at 660 nm, in order to fulfill the bandwidth requirement and preserve 5 information content as much as possible. HARP instruments (AirHARP, HARP CubeSat, and HARP2) use a modified threeway Phillips prism located after the front lens to split the incident light into the three orthogonal linear polarization states (0 • , 45 • , and 90 • ), which can be recombined to obtain the Stokes parameters L t , Q t , and U t at the aircraft altitude (Puthukkudy et al., 2020). Circular polarization (Stokes parameter V) is not measured by any of the polarimeters in ACEPOL as it is negligible for atmospheric studies (Kawata, 1978). We use the total measured reflectance (ρ t (λ)) and DoLP (P t (λ)) at the 10 height of the aircraft with spectral dependencies hereafter implied, which are defined as where F 0 is the extraterrestrial solar irradiance, µ 0 is the cosine of the solar zenith angle, r is the Sun-Earth distance correction factor in astronomical units. 15 Based on the MAP measurements, the MAPOL retrieval algorithm is developed to derive both the aerosol properties and the water leaving signal simultaneously. The retrieval algorithm minimizes the difference between the MAP measurements and the forward model simulations computed from vector radiative transfer simulations (Zhai et al., 2009(Zhai et al., , 2010. By assuming the measurement and modeling uncertainties follow Gaussian statistical distributions, the retrieval parameters can be estimated through Bayesian theory using the cost function χ 2 to quantify the difference between the measurement and the forward model 20 simulation (Rogers, 2000): where ρ t and P t are the measured reflectance and DoLP as defined in Eqs. (1) and (2), and ρ f t and P f t are the corresponding quantities computed from the forward model. The state vector x contains all retrieval parameters, such as the aerosol size, refractive indices; the subscript i stands for the index of the measurements at different viewing angles and wavelengths; and 25 N is the total number of the measurements used in the retrieval. For AirHARP measurements, the maximum value of N is 120 considering all the viewing angles from the four bands. The total uncertainties of the reflectance and DoLP used in the algorithm are denoted as σ ρ and σ P , which are contributed by both the measurement uncertainties σ m and the forward model uncertainties σ f (more details in Sect). 3.3: 5 https://doi.org /10.5194/amt-2020-507 Preprint. Discussion started: 9 February 2021 c Author(s) 2021. CC BY 4.0 License.
One important component of σ m is the calibration uncertainty. AirHARP was calibrated in the lab with an accuracy of 3-5 % for reflectance, and 0.005 for DoLP (McBride et al., 2019). In-flight uncertainty for the AirHARP DoLP is conservatively estimated to be at most 0.01 without an onboard calibrator. In this study, we adopted the calibration uncertainty for reflectance as σ ρ,cal = 3%ρ t and for DoLP as σ P,cal = 0.01 for all four bands. The accuracy of the HARP CubeSat and HARP2 measurements can be further improved through onboard calibration (McBride et al., 2020;Puthukkudy et al., 2020). AirHARP 5 conducted high spatial resolution measurements with a grid size of 55 m. We averaged every 10x10 pixels together (a box of 550m × 550 m). The standard deviations of the pixels within the box are used to estimate the random noise and the spatial variability of the geophysical properties (σ avg ). To account for the total measurement uncertainties, we considered the contributions from both the calibration (σ cal ) and pixel averaging (σ avg ) as for both reflectance and DoLP. As observed by AirHARP (Puthukkudy et al., 2020) and RSP measurements , the sunglint angular pattern cannot be well modeled by an isotropic Cox-Munk model. To minimize the impact of sunglint in our discussions, we removed the signals within an angle range of 0 • to 40 • relative to the solar specular reflection direction.
Furthermore, noise correlation is an import influence on the retrieval accuracy (Knobelspiesse et al., 2012) that is ignored in this study due to the lack of knowledge on this characteristic for AirHARP. 15 The forward model uncertainties σ f include the numerical accuracy of the radiative transfer calculation, and can include any estimation of uncertainties due to the incompleteness of the model to describe the system. For convenience, as discussed in the next section, we will only consider the uncertainty of the NN forward model (σ N N ) and the numerical accuracy of the radiative transfer simulation used for generating the NN training data (σ RT ): 20 To fully utilize the information contained in the AirHARP measurements, the forward model needs to achieve an accuracy level much better than the measurement uncertainties. This becomes the goal of the NN training in the next section, to reproduced the forward model within an error that is much less than the measurement uncertainty. Detailed comparisons of the forward model uncertainties and the measurement uncertainties will be provided in the next section. To minimize the cost function defined in Eq.
(3), we use an optimization method, called the Subspace Trust-region Interior Reflective (STIR) ap- 25 proach (Branch et al., 1999) as implemented in the Python SciPy package (Virtanen et al., 2020), to solve the state parameters x iteratively. The method is based upon the Levenberg-Marquart method (Moré, 1978) and shows good stability for the boundary constraints.

Forward model
We used a vector radiative transfer model based on the successive order of scattering method for coupled atmosphere and ocean 30 systems (Zhai et al., 2009(Zhai et al., , 2010 to model the measured reflectance and DoLP. The atmosphere is configured as three layers: a top molecular layer above the aircraft with only trace gas presented, a molecular layer below the aircraft in the middle, and an 6 https://doi.org/10.5194/amt-2020-507 Preprint. Discussion started: 9 February 2021 c Author(s) 2021. CC BY 4.0 License.
aerosol and molecular mixing layer on the bottom (assumed uniformly distributed within 2 km from the ground) (Gao et al., 2019) as shown in the left panel of Fig. 1. The same vertical structure of the atmosphere was successfully used in the inversion of RSP data (Gao et al., 2019.
Aerosols are diverse in size, composition, and morphology. To capture their variation in the atmosphere, we modeled the size and refractive index for both fine and coarse modes. The aerosol size is represented by the volume density distribution as 5 a combination of five lognormal distributions: where V i is the column volume density for each submode, the mean radius r i and standard deviation σ i are fixed with values of 0.1, 0.1732, 0.3, 1.0, 2.9 µm, and 0.35, 0.35, 0.35, 0.5, 0.5 respectively (Dubovik et al., 2006;Xu et al., 2016). The first three submodes are categorized as the fine mode aerosol, while the last two submodes are the coarse mode. All aerosols are assumed 10 to be spherical in the current forward model. The nonspherical particle shape is important in the aerosol model (Dubovik et al., 2006), and will be considered in future studies. The aerosol refractive index spectra for the fine and coarse modes are represented by the principal component analysis in MAPOL (Wu et al., 2015;Gao et al., 2018) as where p 1 (λ) is the first-order principal component computed from the aerosol refractive index dataset including water, sea salt, 15 dust-like particles, biomass burning, soot, sulfate, water-soluble, and industrial aerosols (Shettle and Fenn, 1979;d'Almeida et al., 1991). m 0 and α 1 are two coefficients to determine the spectrum. For the application to AirHARP bands, p 1 (λ) for the real part of the refractive index is approximately spectral flat for both the fine and coarse mode aerosols within the visible spectrum. We further assume the spectral shape for the imaginary refractive spectra is also flat. Therefore, the two parameters can be combined into one to represent the refractive index. In this study, only four independent parameters are used to determine 20 the real and imaginary refractive index spectra for the fine and coarse modes. With the aerosol size and refractive index, the polarimetric single scattering properties are modeled by the Lorenz-Mie theory and computed by the code developed by Mishchenko et al. (2002).  Table 1. The molecular absorption properties are computed by a hyperspectral radiative transfer simulation (Zhai et al., 2009(Zhai et al., , 2010, including contributions from ozone, oxygen, water vapor, nitrogen dioxide, methane, and carbon dioxide under the US standard atmospheric constituent profiles (Anderson et al., 1986). Ozone is the most important gas that influences the absorption transmittance at the AirHARP bands of 550nm and 660nm. For the application to AirHARP measurements in ACEPOL, we use the ozone column density as a free parameter with values from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) developed by NASA's Global Modeling and Assimilation Office (Gelaro et al., 2017) to rescale the molecular absorption optical depth calculated under the above-mentioned standard atmospheric profile.
For the ocean layer shown in Fig. 1, two ocean bio-optical models are implemented in the forward model of MAPOL: one with chlorophyll-a concentration (Chla; mgm −3 ) as the single parameter applicable to open ocean optical properties, and 5 the other with seven parameters more suitable to fully describe complex coastal waters (Gao et al., 2019). Since the waters are mostly clear within the ocean scenes in this study , open ocean model is used for both NN training and retrievals. The optical properties of open ocean waters include contributions from pure seawater, colored dissolved organic matter (CDOM), and phytoplankton, where the CDOM and phytoplankton absorption coefficients, and phytoplankton scattering coefficient and phase function are parameterized by Chla (Gao et al., 2019). Complex costal water model for NN trainings 10 will be investigated in future studies. The ocean surface roughness is modeled by the isotropic Cox-Munk model with a scalar wind speed.
In summary, the parameters used to represent the forward model include five volume densities (one for each submode), four independent parameters for the refractive indices of fine and coarse modes, one parameter for wind speed, and Chla. Three additional geometric parameters are used to set up the system, including the solar zenith angle, viewing zenith angle, and 15 relative viewing azimuth angle. Therefore, it requires a total of 15 parameters to conduct the radiative transfer calculation, with a total of 11 independent state parameters that can be retrieved from optimizing the cost function as defined in Eq. (3).

Remote sensing reflectance
An important task for the joint retrievals is to obtain the water leaving signal, which is often represented in ocean color studies by the spectral remote sensing reflectance defined as R rs = L + w /E + d where E + d is the downwelling irradiance and L + w is the 20 water leaving radiance just above the ocean surface (Mobley et al., 2016). The remote sensing reflectance can be derived from the water leaving reflectance reaching the sensor (ρ w ) via: where θ 0 and θ v are the solar and viewing zenith angles. ρ w represents the signals originating from scattering in the ocean that reached the sensor, which can be derived from the atmospheric correction process as where ρ t is the measured total reflectance as defined in Eq.
(1) and ρ f t,atms+sf c is the reflectance from a system with only atmosphere and ocean surface (Mobley et al., 2016) as represented in the right panel of Fig. 1. The same formalism has been used to derive R rs from RSP measurements Gao et al. (2019Gao et al. ( , 2020. The downwelling irradiance transmittance T f d is for the solar irradiance from TOA to the surface, and the upwelling radiance 30 transmittance t f,+ u is for the water leaving radiance from BOA to the sensor (Gao et al., 2019). Both T f d and t f,+ u are denoted in Fig. 1 and represented as follows: where ρ f,+ t , ρ f,+ t,atm+sf c are reflectance just above ocean surface also denoted in Figure 1 and Table 1. To remove the dependency of R rs on the solar and viewing geometries, a BRDF correction C BRDF is applied to adjust R rs 5 to the observation with a zenith sun and a nadir viewing direction as defined by Morel et al. (2002): where R o /R accounts for reflection and refraction effects when light propagates through the ocean interface. C BRDF in its original form is defined using the radiance and irradiance just below the ocean surface (Morel et al., 2002), here we have converted all quantities into the radiance reflectance ρ f,− t and the irradiance transmittance (1) and 10 (12). In our study, of AirHARP retrievals of remote sensing reflectance, we only consider the minimum viewing angle at each wavelength which is less than 1 • , the contribution of the R factor is ignored due to its small angular variations (Morel and Berthon, 1989). All quantities denoted in Fig 1 need to be determined for the forward model and the calculation of remote sensing reflectance, and will be represented by NN models.
3 Neural network for forward model 15 Deep NN models are developing rapidly due to the advancement in machine learning infrastructure and demands in broad applications (Goodfellow et al., 2016), and are demonstrated to be efficient in approximating physical functions (Lin et al., 2017). In this study, we employed the deep feed-forward NN (Goodfellow et al., 2016) to represent the MAP measurements. In this section, we will discuss the procedures to train the NN forward model for AirHARP measurements, with its performance evaluated.

Training data
To train a NN that can represent the forward model accurately for the AirHARP measurements from the ACEPOL field campaign, we generated the training data according to the average aircraft height of 20.1km on the day of 10/23/2017 from ACEPOL. We simulated 21,000 cases according to the forward model as discussed in the previous section by considering general aerosol and ocean properties, as well as a large range of solar and viewing geometries with the minimum and maximum 25 values of all parameters summarized in Table 2. The range of solar zenith angle θ 0 , viewing zenith angle θ v , and relative viewing azimuth angle φ v are from 0 • to 70 • , 60 • and 180 • , respectively. The reflectance and DoLP with a viewing azimuth angle larger than 180 • can be evaluated by the corresponding value less than 180 • due to symmetry with respect to the principal plane (defined by φ v = 0 • and φ v = 180 • ). For each solar zenith angle, the polarized reflectance is calculated for all viewing angles within the aforementioned ranges with an angular resolution of 1 • . The solar zenith angle, ozone column density, refractive index, and wind speed are randomly sampled in a linear scale. Chla is randomly sampled in a log scale. The fine mode volume fraction is sampled uniformly within [0, 1], which is then randomly partitioned to each submode. To maintain a uniform distribution of the total AOD, we sampled the AOD at 550nm within [0,0.5] in a linear scale. The volume density V i of each submode is determined by the total aerosol optical depth and volume fraction for each mode. Fig. 2 shows one example 5 simulation dataset for the angular distribution of reflectance and DoLP.  We randomly selected 20,000 cases out of the total 21,000 simulated cases for the training and validation processes, and the remaining 1000 random cases will be used as test cases to evaluate the NN accuracy, which will be discussed in the next section. To enable the NN to predict reflectance and DoLP at any given viewing geometry, for each case, we sampled 100 random pairs of viewing zenith and azimuth angles. If the sampled angles fall outside of the pre-defined angular grids, values from spline interpolation are used. The sunglint angles within an angle of 40 • to the solar specular reflection direction are 5 removed. Approximately 1 million data points are obtained for each wavelength for training.

Neural network training
A feed-forward NN can be defined recursively with one input layer, one output layer, and k hidden layers (Aggarwal, 2018): where x is the input parameter vector including all 15 parameters needed to define the forward model as listed in Table 2.
Here x not only contains the retrieval parameters in the state vector defined in Eq.
(3) but also include additional non-retrieval parameters such as the solar zenith angle, viewing zenith and azimuth angles and the ozone column density. y is a fourdimensional output vector for reflectance or DoLP at the four AirHARP bands. The weight matrix W p+1 connects the p-th 15 and (p+1)-th NN layers. The bias vector for the (p+1) layer is defined as b p+1 . The output of each layer h p+1 becomes the input of the next layer as shown in Eq. 16. k is the number of hidden layers, k+1 refers to the output layer. In this study, we tested several NN architectures and eventually chose three hidden layers with the number of nodes of 1024, 256 and 128 as shown in Table 3. The nonlinear activation function Φ used in this model is the LeakyReLU function. For each vector element, it is defined as The training process is to minimize the cost function defined as the mean square error between the training data generated from radiative transfer simulations and the NN predicted values (Aggarwal, 2018). All parameters in the neural network 5 weight matrices and bias vectors, over 670,000 numbers, need to be trained. With this large number of parameters, it is a challenging task to avoid overfitting where the model works well for the training dataset but poorly for the dataset not used in the training process. Several training procedures are performed for reflectance and DoLP data to avoid overfitting and improve NN performance: 1. Both input and output data are normalized before training. We normalize the input data into the range of [0,1] using the 10 minimum and maximum values from the datasets as listed in Table 2. The reflectance and DoLP in the output layers are normalized by dividing their standard deviation of the training data at each wavelength. 3. The learning rate determines the step size in the parameter update. We use an exponential decay schedule to reduce the learning rate: we start with a learning rate of 0.005 and reduce the learning rate by a factor of 10 every 200 epoches.
4. To monitor over-fitting in the training process, we split the data into 70% for training and 30% for validation. We conduct 20 the optimization based on the training dataset, and meanwhile, we monitor the performance of training by applying the NN model on the validation dataset. To avoid overfitting, the early-stopping approach is employed where the training is stopped when the cost function on the validation dataset stops to reduce for a threshold of 50 epochs.
The machine learning Python library Pytorch is used for the training (Paszke et al., 2019). The trained NN model is used to replace the radiative transfer model to compute the reflectance and DoLP in the retrieval algorithm. The Jacobian matrix used 25 in the optimization is computed by the finite difference approximation of the partial derivatives of reflectance and DoLP with respect to the retrieval parameters. Here central difference method is used. Note that the Jacobian matrix can also be computed analytically from the NN model using the automatic differentiation techniques based on the chain rule of differentiation (Baydin et al., 2018). This will be a topic in our future studies.

Neural network accuracy 30
After training the NN model, we evaluated its accuracy using synthetic AirHARP measurements generated from the 1000 simulation cases which have not been used in the training and validation process. Each simulation dataset includes polarized reflectance on regular viewing angle grids, which are interpolated to the viewing geometry of AirHARP to create synthetic measurement data and compare with the NN predictions. Glint angles are excluded from the comparison because the NNs are not trained over these angles. As one example shown in Fig. 4, both the reflectance and DoLP are in good agreement between the synthetic data and the NN results, where the maximum absolute differences for reflectance and DoLP are within 0.001 and 0.0025. This translates to a difference for both reflectance and DoLP mostly less than 1%. The maximum percentage difference 5 can be as large as 3% for 870nm bands due to the small reflectance magnitude.  The comparison with all 1000 synthetic datasets and their NN predictions are shown in Fig. 4. The mean absolute error (MAE), and the root mean square error (RMSE) between the simulation data (T i ) and the NN predicted data (R i ) shown in Fig. 4 are defined as Both MAE and RMSE are useful metrics, where MAE has less dependency on outliers comparing with RMSE. Analysis shows that the statistics of the differences between the NN prediction and the RT simulations as shown in Fig.   4 can be well modeled by Gaussian distributions and characterized by RMSE. Therefore the RMSE is used to represent the NN uncertainties for both reflectance (σ ρ,N N ) and DoLP (σ ρ,N N ) and will be incorporated into the total uncertainties in the cost function. Table. 3 summarizes the uncertainties of the NN models. The σ ρ,N N at 440 nm is 0.0006, which decreases to 0.0004 at 870 nm. However, due to the smaller reflectance magnitude at 870nm, the corresponding RMSE for the percentage 5 reflectance difference as shown in Fig. 4 is increased from 0.4% at 440nm to 1.0% at 870nm. For DoLP, the maximum σ P,N N is 0.003 at 870 nm which decreases to 0.0016 at 440 nm. The uncertainties can be further improved with more training data points. For the readers' information, RMSE of the NN model trained with 20,000 cases (1 million data points) decreases by a factor of √ 2 in comparison with the one using 10,000 cases (0.5 million data points). It takes 0.01s in a single-core CPU (AMD EPYC Processor) or 1 ms in a GPU (GeForce GTX 1060) to predict all 120 angles for both reflectance and DoLP in the 10 NN forward model. The assessment of the NN accuracy is relative to the synthetic measurements simulated by the vector radiative transfer simulations. To account for the modeling uncertainties of the forward model σ f , we consider both the NN accuracy σ N N and the numerical accuracy of the radiative transfer simulations σ RT for reflectance and DoLP, respectively. Uncertainties due to incomplete assumptions in the forward model are not considered. Several internal parameters determine the numerical 15 accuracy of the radiative transfer simulations. In the framework of the successive order of scattering (Zhai et al., 2008(Zhai et al., , 2009, these parameters include the number of scattering orders (Ns), the number of Gaussian quadratures for discretizing the viewing zenith angle in the atmosphere (P a ) and ocean (P o ), and the order of Fourier decomposition (M ) for the viewing azimuth angle, and the order of Legendre expansion (L) of the single scattering phase function. In this study, we chose N = 20, P a = 32, P o = 64, M = 32, and L = 32, which has a higher accuracy than the radiative transfer forward model directly used 20 in our previous retrieval studies Gao et al. (2020).
To quantify the accuracy of the radiative transfer calculation used for generating training data (σ RT ), we simulated an additional 1000 synthetic AirHARP datasets with all internal parameters doubled as the most rigorous calculations, and the viewing angular resolution was reduced from 1 • to 0.5 • in order to reduce interpolation errors. The resultant reflectance and DoLP values are compared between these two sets of radiative transfer calculations. The RMSE for each band can be used as a measure of the accuracy for the radiative transfer calculation used to generate the training data (σ RT ). The uncertainties σ RT 5 for reflectance and DoLP are summarized in Table 4, with reflectance uncertainties less than 0.00015 and DoLP uncertainties less than 0.0007 for all AirHARP bands. σ ρ,RT is about four times smaller than the NN uncertainties, and σ P,RT is about 4 to 10 times smaller. Therefore, NN uncertainties are not dominated by the uncertainties of the RT simulation. The measurement uncertainties from calibration (σ cal ) and pixel averaging (σ avg ) as discussed in Section 2 are also summarized in Table 4, which shows the total forward model uncertainties σ f = σ 2 RT + σ 2 N N as approximated in Eq. (7) are much smaller than the 10 total measurement uncertainties σ m = σ 2 cal + σ 2 avg as defined in Eq. (6). The overall uncertainties used in the retrieval cost function in Eq (3) are dominated by the measurement contributions as we expected. Table 4. Comparisons of the uncertainties for reflectance (ρ) and DoLP (P) for both measurement and forward model including calibration uncertainty (σ cal ), the uncertainty from averaging AirHARP pixels into a 550 m × 550 m box (σavg), the radiative transfer simulation uncertainty (σRT ), and the NN uncertainty (σNN ). Same σρ,NN and σP,NN have been shown in Table 3, and are repeated here for comparisons.
(Same as Furthermore, in this study higher accuracies from the radiative transfer simulations are used for the NN training for comparison with the accuracies from the radiative transfer model directly used in our previous retrieval algorithm. Since the simulations of the training data can be conducted independent to the retrieval algorithm, higher computational costs can be accommodated 15 to improve NN forward model accuracy. After the NN model is trained, the model can be applied to the retrieval algorithm through efficient matrix operations.

Neural network model for remote sensing reflectance
As discussed in Sect. 2.2, the water leaving signals are represented by the remote sensing reflectance as defined in Eq. (10) ( Mobley et al., 2016). To conduct the atmospheric correction in Eq. (11), we need to determine the reflectance ρ f t,atmos+sf c at the aircraft level, transmittance t f u and T f d , and the BRDF correction coefficient C BRDF . Based on Eq. 10, we combined T f,+ d , t f,+ u , and C BRDF into a single term denoted as [C BRDF /T d t u ]. To efficiently determine R rs , two NNs need to be trained to 5 represent ρ f t,atmos+sf c and [C BRDF /T d t u ], respectively, Following similar NN training schemes as discussed previously, we conducted 10,000 simulations to determine the reflectance at aircraft altitude ρ t,atmos+sf c from a system with only atmosphere and ocean surface (right panel of Fig. 1), and trained the NN for ρ t,atmos+sf c in the same way as ρ f t . Since this system only includes atmosphere and ocean surface but without ocean body, there are total 14 input parameters (without Chla). To train a NN for [C BRDF /T d t u ] with T f,+ d , t f,+ u , and 10 C BRDF defined in Eqs. (12), (13) and (14), we obtained five additional quantities corresponding to the above-mentioned 10,000 cases with and without ocean body: for the fully coupled system with atmosphere, ocean surface and ocean body (left panel of Fig. 1), we computed the reflectance just above and below the ocean surface (ρ f,+ t and ρ f,− t ), and irradiance transmittance just above and below the ocean surface (T f,+ d and T f,0 d ); for the system without ocean body but with ocean surface (right panel of Fig. 1), we computed the reflectance just above the ocean surface (ρ f,+ t,atms+sf c ). The accuracy of the NNs for ρ f t,atmos+sf c and

15
[C BRDF /t u T d ] are evaluated and shown in Table 3, which are in the same order of uncertainty magnitudes in percentage as other quantities.
To evaluate the overall accuracy for the R rs after the BRDF correction, we conducted radiative transfer simulations with a zenith sun and a nadir viewing direction, and obtained the truth remote sensing reflectance using the upwelling radiance and downwelling irradiance just above the ocean surface as examples shown in Fig. 5. The predicted R rs were computed from Eq. 20 (10) after the application of two NNs. The RMSE of the difference between the simulated and NN predicted R rs are shown in Table 3

Joint retrieval results on synthetic AirHARP measurements
The NN forward models for reflectance (ρ f t ) and DoLP (P f t ) are used in the FastMAPOL retrieval algorithm as discussed in Sect. 2. To evaluate the performance of the retrieval algorithm, we conducted retrievals on the synthetic AirHARP data. The creation of the synthetic data is discussed in Sect. 3.3. To account for the measurement uncertainties, random noise is added to the simulated data according to the calibration uncertainties as listed in Table 4. The total uncertainties in the cost function 5 include contributions from calibration (σ cal ), radiative transfer simulation (σ RT ), and NN model (σ N N ). Uncertainties from pixel averaging (σ avg ) for the AirHARP field measurements are not considered in the synthetic dataset.
Using the initial values as listed in Table 2, a total of 1000 synthetic AirHARP cases are retrieved with the cost function values(χ 2 ) summarized in Fig 6. Retrievals with χ 2 < 1.5 are chosen in our following discussion, which includes 96% of all retrieval cases. Gao et al. (2020) showed that the retrieval results depend on the initial values. Testing with several random sets  Table 3 with a total of 1000 cases. The most probable χ 2 is 0.82. A threshold of χ 2 < 1.5 is used in the discussion.
With the directly retrieved aerosol refractive index and volume densities (see Table 2) as inputs, the aerosol optical depth (AOD) and single scattering albedo (SSA) for both the fine and coarse modes, were computed using additional NNs to represent the Lorenz-Mie calculations in Appendix A. The retrieved total AOD, SSA, wind speed, and Chla are compared with the truth values as shown in Fig. 7. Total AOD indicates the summation of the fine and coarse mode AODs and total SSA is the ratio of the total scattering and extinction cross-sections, both are specified in in Appendix A. For fine aerosol, the AOD, SSA, 5 refractive index (m r ), effective radius (r ef f ) and variance (v ef f ) are shown in Fig. 8. The color plots indicate the data point density (normalized by its maximum value) approximated by a kernel density estimation method (Silverman, 1986).
In order to quantify the variation of the retrieval uncertainties with respect to different aerosol loadings, we computed the fine mode volume fraction are uniformly sampled for the simulated data, therefore, there is an equal mixing fraction of fine and coarse mode aerosol for each AOD range. The retrieval uncertainties for aerosols are shown in Figure 9 with the corresponding ranges indicated by AOD values from 0.1 to 0.5. All discussions regarding the AOD and SSA are for wavelength of 550nm in this section.
As shown in Fig. 7 and Fig. 9, the errors of the retrieved total AOD increase with aerosol loadings: the uncertainty (  For wind speed retrievals as shown in Fig. 7 the agreements between the truth and retrievals depend strongly on the wind speed value itself: when the wind speed is small, there is less retrieval sensitivity due to the removal of glint; for larger wind speed, the agreements are improved, likely due to the larger range of angles influenced by wind speed. The retrieval uncertainties are shown in Fig. 9, for wind speed (WS) smaller than 3 m/s, the uncertainty increases from 1.5 to 2.1 m/s for AOD ranges from [0.1,0.2] to [0.4, 0.5]. While for wind speed larger than 3 m/s, the retrieval uncertainty is almost an constant 5 of 1.2 m/s with a small increase less than 0.1 m/s. The retrieved and truth Chla is compared in Fig. 7 where the MAE in log scale, or MAE(log), is used with definition as: 22 https://doi.org/10.5194/amt-2020-507 Preprint. Discussion started: 9 February 2021 c Author(s) 2021. CC BY 4.0 License.
where R i and T i denote the retrieval and truth values. MAE(log) is recommended by Seegers et al. (2018) as a better metrics for Chla, which indicates the averaged ratio between the retrieval and truth values. The dependency of the MAE(log) for Chla with the aerosol loadings is shown in Fig. 9. The Chla retrieval performance depends on the magnitude of the Chla.
In this work, we chose four ranges of Chla according to the trophic regions discussed in (Seegers et al., 2018). Note that Chla from in-situ measurements is typically lager than 0.01mg/m 3 and we chose a broader range of Chla with its minimum 5 value of 0.001mg/m 3 as listed in Table 2 for sensitivity studies. For 0.01mg/m 3 < Chla < 0.1mg/m 3 and 0.1mg/m 3 < Chla < 1mg/m 3 , Chla retrieval uncertainties vary within 1.3 to 1.6 when AOD<0.3, and then increase to 2.3 at AOD range [0.4,0.5]. For Chla > 1mg/m 3 and Chla < 0.01mg/m 3 , the uncertainties are generally larger with a value around 2 to 3.
With the retrieved aerosol and ocean properties, the atmospheric correction procedures can be applied to compute the remote sensing reflectance as discussed in Section 3.4. The comparison of the retrieved R rs with the truth data is shown in Fig. 5. To 10 account for the various solar geometries, the BRDF correction has been applied on the retrieved R rs as discussed in Section 3.4. The truth R rs was computed with a zenith sun and a nadir viewing direction, emphasizing the need for the latter correction to the MAP observations. Overall R rs uncertainties for the four bands are 0.007, 0.0004, 0.0002 and 0.0002 as shown by the RMSE in Fig. 10. MAE showed values of 0.0006, 0.0003, 0.0002 and 0.0001, which are less sensitive to the outliers. Note that the atmospheric correction is applied on the synthetic measurements without adding additional random noise in order to 15 evaluate the impacts on R rs uncertainties from only aerosol and ocean surface properties retrievals. The retrieval uncertainties for R rs for each AirHARP bands are shown in Fig. 10 depending on the aerosol loadings: larger uncertainties are found with larger aerosol optical depth.
The PACE accuracy requirements on ocean color are specified in terms of the water-leaving reflectance, which can be converted to those of R rs by dividing them by a factor of π. The resultant requirements in terms of R rs are 0.0006 sr −1 or 5% 20 from 400 to 600 nm, and 0.0002 sr −1 or 10% from 600 to 710 nm (Werdell et al., 2019). As shown in Fig. 11, R rs at 550nm are within the requirement of 0.0006 sr −1 for all AOD ranges. For 440 nm band, the R rs retrieval uncertainties are larger than the requirement when AOD is larger than 0.3. R rs at 670 and 870 nm varies in a very small dynamical range and has less impacts by the aerosol retrievals. R rs uncertainties at 670 and 870 nm are slightly larger than the requirement of 0.0002 sr −1 when AOD(550nm) is larger than 0.4 and 0.3 respectively. Further work is needed to understand how the uncertainties of the 25 retrieved aerosol properties influence the retrievals. Note that from Table 3, the uncertainties of the R rs computed using NNs have an uncertainty of 0.0004 to 0.0001 from 440nm to 870nm, which may be further reduced with better training and help the reduction of the R rs retrieval uncertainties.

Joint retrieval results on AirHARP measurements from ACEPOL
The ACEPOL field campaign, conducted from October to November of 2017, included a total of six passive and active instruments on the NASA ER2 high-altitude aircraft  with four MAPs: AirHARP (McBride et al., 2020)), AirMSPI (Diner et al., 2013), SPEX Airborne  and the RSP (Cairns et al., 1999), and two lidars: HSRL-2 (Burton et al., 2015) and CPL (the Cloud Physics Lidar) (McGill et al., 2002). Aerosol retrieval algorithms have been applied 5 for all four MAPs Puthukkudy et al., 2020;Gao et al., 2020). The measurement datasets are available from the ACEPOL data portal . In this work, we focus on the study of the AirHARP measurements over ocean scenes as shown in Fig. 12 on Oct 23, 2017. The viewing angles are within ±57 • along-track, and ±47 • cross-track as shown in the polar plots in Fig. 12. Fig. 13 shows the RGB images (670, 550 and 440nm) for the three scenes at near nadir viewing direction.  The HSRL-2 instrument from ACEPOL provided useful aerosol optical depth ground truth at 355 and 532 nm (Hair et al., 2008;Burton et al., 2016), which was used for the validation of the aerosol retrieval algorithm using the AirHARP data.
The HSRL-2 measures the pixels along the track as shown in Fig. 12, where an assumed lidar ratio of 40 sr is multiplied by the aerosol backscatter coefficient derived from the HSRL technique to produce aerosol extinction and AOD at 532 nm. For the low aerosol loading cases considered in this study, the assumed lidar ratio approach produces a systematic uncertainty 5 of ±50% . In scene 3, the aircraft also flew over an AERONET USC_SEAPRISM site, which is equipped with a CIMEL-based system called the Sea-Viewing Wide Field-of-View Sensor (Sea-WiFS) Photometer Revision for Incident Surface Measurements (SeaPRISM) that collects radiances at eight wavelengths of 412, 443, 490, 532, 550, 667, 870, and1020 nm (Zibordi et al., 2009). AOD from the AERONET data product version 3 level 2.0 data was used in this study, which is also consistent with the HSRL-2 AOD at 532nm as shown latter in Fig. 20. The estimated AERONET AOD uncertainty is 10 from 0.01 to 0.02 with the maximum uncertainty in the UV channels (Giles et al., 2019). We compared AOD from AirHARP retrievals with those from both HSRL and AERONET. Furthermore, to validate the atmospheric correction procedure in the retrieval algorithm, we compared the retrieved remote sensing reflectance with the AERONET ocean color (OC) products as reported by the SeaPRISM measurements at the USC_SEAPRISM site.
Here we applied the FastMAPOL algorithm to the AirHARP field measurements from ACEPOL. The solar and viewing 15 geometries as shown in Fig. 12, and the ozone column density from MERRA2 are the inputs to the forward model. As discussed in Section 2, the total uncertainties are modeled as σ 2 = σ 2 cal + σ 2 avg + σ 2 RT + σ 2 N N for reflectance and DoLP respectively, with all values listed in Table 4.
The histograms of χ 2 for all pixels retrieved in each scene are shown in Fig. 14. Comparing with the histogram in synthetic data retrievals in Fig. 6, the histograms in Fig. 14 show a longer tail with larger χ 2 value. This may due to that some pixels are can be seen as a stripe, due to the removal of angles influenced by water condensation in the lens, which can also be observed 5 in the polar plots in Fig. 12 with the nadir region removed. All three figures show that the number of viewing angles are smaller at the edges parallel to the flight track, where small χ 2 can be achieved but may be less reliable.
For pixels with large χ 2 , the forward model cannot fit reflectance or DoLP to the measurements, which may be due to the contamination by cloud (Stap et al., 2015), land, or residuals of glint. The large χ 2 values are also correlated with the large retrieved AOD(550nm) values, for instance, the central region in scene 1. We have excluded retrievals with less confidence and 10 only discussed the retrievals simultaneously satisfying the two criteria with N v > 10 and χ 2 < 1.5. Overall about 40% to 60% pixels are available after applying these two criteria.   Figure 17. Same as 15 but for scene 3, 94% pixels with Nv > 10, 47% pixels with χ 2 > 1.5, and 43% pixels under both conditions. The pixels with large χ 2 are mostly influenced by the land (upper region) and island (lower left).
The AOD at 532nm from HSRL data product in all three scenes is in the range from 0.03 to 0.1 as shown in Fig. 18.
After applying the two criteria in N v and χ 2 , many pixels from AirHARP retrievals along the HSRL track are not available.
To compare with the HSRL AOD along the track direction, the retrieved AOD is averaged in the cross-track direction, the averaged values and its standard deviations are shown in Fig. 18. The averaged and standard deviation of all pixels satisfying the criteria are also shown in the plots. For scenes 1-3, the averaged AOD(550nm) are 0.076, 0.066 and 0.052. The averaged HSRL AOD is 0.079, 0.072 and 0.038. The average AOD values retrieved from FastMAPOL and those from HSRL agree well for both scene 1 and scene 2. For scene 3, the average AOD value of the FastMAPOL algorithm is larger than the averaged HSRL AOD by 0.014. This may be due to complex water properties not well represented by the open water bio-optical model used in the simulation (Gao et al., 2019). Furthermore, the pixels near the coast are potentially impacted by the adjacency effect of land pixels. However it is challenge to investigate due to the small aerosol loadings. Furthermore, thin cirrus clouds were  for scene 1 and 2 show a larger value of 0.005 to 0.007 at 440nm but reduced to 0.0025 for scene 3, which is closer to the coast.
The decrease of R rs may be due to the increase of CDOM as its spectra demonstrated in Fig. 5. The AERONET AOD is compared with the FastMAPOL retrieved AOD at 550nm in Scene 3 of Fig. 18, and the whole AOD spectra are compared in Fig. 20. The uncertainties of the AOD and R rs from retrieval are estimated by the average of 2×2 pixels (in 1.1km×1.1 km box) around the AERONET site. The retrieved AODs at 440nm and 870nm agree well with the 5 AERONET results. The difference between the HSRL AOD at 532nm and the AirHARP retrieved AOD near the AERONET site and is smaller than the one shown in the Scene 3 of Fig. 18, which suggests larger differences are contributed by the pixels further away. The retrieved R rs agrees well with the AERONET R rs with a value slightly larger than the AERONET results at 410nm and 550nm by 0.0005 and 0.0006 sr 1 . The retrieved R rs has a small negative value at 670nm, lower than AERONET by 0.0007 sr 1 . Their difference at 868 is 0.0002 sr 1 . There are larger uncertainties when larger spatial box is considered. Note 10 that this study is done with the AirHARP measurement with 3% uncertainties in the reflectance measurement, which may impact the accuracy in the comparison between the retrieval and in-situ measurements.

Discussion
The NN model greatly improved the speed of the forward model used in the iterative optimization approach and boosted the efficiency of the FastMAPOL retrievals. The average retrieval speed for one pixel with FastMAPOL is approximately 2.6 15 seconds with a single CPU core, or approximately 0.3 seconds with a GPU using the same hardware as mentioned in Sect 3.3. Meanwhile, the NN model maintains a high accuracy as shown in Table 3 and Table 4. For retrieval algorithms running radiative transfer simulation, the accuracy is often tuned down to reduce simulation time. By training a NN, however, the high accuracy of the radiative transfer model simulation can be achieved, as has been demonstrated in this work. Thus, despite the one time high computational costs in generating the training datasets and conducting the training, the trained NN can be applied efficiently to large observational datasets in the retrieval algorithm.
In the retrieval of the AirHARP field measurement, each pixel has multiple viewing angles that are aggregated from measurements at different times with slightly different solar zenith angles. The NN used in FastMAPOL computes the polarimetric measurement for specific solar zenith angles for each viewing direction, and therefore, the variation of the solar zenith angle 5 can be captured. These effects may be small for AirHARP measurement in ACEPOL, with a maximum solar zenith angle difference within 0.6 degrees. However, this capability can help to minimize the impacts of the solar angle for HARP2 in space-borne measurements, which can reach to a maximum difference of 1.5 • for HARP2 observations.
With the efficient retrievals from FastMAPOL, we have discussed the retrieval performance and uncertainties for the aerosol properties, including AOD, SSA, refractive index, and particle sizes. Since the AirHARP measurements share many similar 10 characteristics with HARP2 as planned for the PACE mission, the knowledge from the retrieval analysis can help to understand the retrieval performance for the HARP2 instrument in space-borne measurements. Note that HARP2 is likely to have high accuracy due to the onboard calibration capability and the potential to conduct cross-calibration with the OCI instrument.
For the development of the NN forward model for space-borne measurements, similar training procedures can be applied with the sensor altitude at the top of the atmosphere instead of the aircraft altitude used in this study. Due to the impact of 15 retrieval capability by geometry (Fougnie et al., 2020), solar and viewing geometries according to the PACE orbits need to be considered.
The water leaving reflectance is obtained from the atmospheric correction process using the aerosol and ocean properties retrieved from the AirHARP measurements, and a similar procedure can be applied to HARP2 retrievals. Since the hyperspectral OCI in PACE will provide high accuracy measurements, the retrieved information can be applied to OCI and therefore 20 assist hyperspectral atmosphere corrections as demonstrated by Gao et al. (2020); Hannadige et al. (2020). However, aerosol retrieval and atmospheric correction are challenging in the UV spectral range (Remer et al., 2019a). For the ocean bio-optical model in this study, the water properties are modeled as open ocean waters parameterized by a single Chla value. For complex coastal water, complex bio-optical models are preferred in the retrieval of both accurate aerosol properties and water leaving signals as demonstrated by Gao et al. (2019). 25

Conclusions
We have demonstrated the application of a NN for highly accurate forward model calculations of polarimetric measurements for AirHARP. Additional NN models were used to conduct atmospheric correction. These models are used in the FastMAPOL joint retrieval algorithm to conduct simultaneous aerosol property and water leaving signal retrieval. Applications to both the synthetic AirHARP data and field measurements from ACEPOL are discussed. The uncertainties of the retrieved aerosol prop-30 erties and remote sensing reflectance are discussed for different aerosol loadings. These results from AirHARP retrievals can help to evaluate the retrieval capabilities for the HARP2 instrument on PACE. In application to field measurements from ACE-POL, the impacts of the number of viewing angles and the value of cost function to the retrieval quality are discussed. Further properties such as AOD and SSA as discussed in Section 4, we developed additional four NNs to represent the AOD and SSA for both fine and coarse modes, respectively. These NNs are only used to analyze the retrieved aerosol properties and are not used in the retrieval process. The NN architectures and accuracy are shown in Table A1. The input parameters for the fine mode SSA and AOD are the three submode volume densities, and the real and imaginary parts of refractive index, with a total of 5 parameter. For coarse mode aerosols, there are a total of 4 parameters with only two submodes used. The outputs are the AOD 5 and SSA at the four AirHARP bands.
A total of 10,000 training data points are generated in the same way as in Section 3.1 using the Lorenz-Mie code discussed in Section 2. The NN model accuracy is evaluated with additional 1000 data points not used in the training. As shown in Table   A1, the accuracy is much smaller than the retrieval uncertainties shown in Figure 9, therefore the NNs for AOD and SSA provide sufficient accuracy to evaluate the aerosol single scattering properties. 10 With the fine and coarse mode AOD and SSA evaluated, the total AOD and SSA can be derived. The total AOD (τ t ) is the summation of the fine and coarse mode AODs as where τ f and τ c are the fine and coarse mode AODs. The total (or averaged) SSA (ω t ) is defined as the ratio of the total scattering cross-section and the total extinction cross-section for both fine and coarse modes, which can be computed as where ω f and ω c are the fine and coarse mode SSA.