Sensitivity of Ocean Color Atmospheric Correction to Uncertainties in Ancillary Data: A Global Analysis With SeaWiFS Data

Atmospheric correction (AC) algorithms for ocean color (OC) data processing usually rely on ancillary data documenting the atmosphere and the sea state to help the calculation of the remote sensing reflectance <inline-formula> <tex-math notation="LaTeX">$R_{\text {RS}}$ </tex-math></inline-formula> from the radiance measured by a space sensor. This study aims at assessing the impact that the uncertainties associated with these ancillary data have on the AC outputs. For this objective, a full year of global Sea-viewing Wide Field-of-view Sensor (SeaWiFS) imagery is processed with the standard AC algorithm <italic>l2gen</italic> of the National Aeronautics and Space Administration with different sets of ancillary data, the reference case with National Centers for Environmental Prediction (NCEP) Reanalysis-2 meteorological data and satellite ozone products, as well as with ten ensemble members from the European Centre for Medium-Range Weather Forecast (ECMWF) CERA-20C data. The spread within the ensemble data and the differences with respect to the reference case are taken as a measure of the uncertainties associated with ancillary data. The impact on <inline-formula> <tex-math notation="LaTeX">$R_{\text {RS}}$ </tex-math></inline-formula> of perturbations in ancillary variables vary in space, the variables having the largest effects being wind speed and relative humidity, and ozone at bands where ozone absorption is largest, while sea-level pressure and precipitable water have the smallest effect. Sensitivity coefficients quantifying the relationship between perturbations in ancillary variables and effects on <inline-formula> <tex-math notation="LaTeX">$R_{\text {RS}}$ </tex-math></inline-formula> change with variable and wavelength. At the global scale, the variations found on <inline-formula> <tex-math notation="LaTeX">$R_{\text {RS}}$ </tex-math></inline-formula> when ancillary data are perturbed are usually small but not negligible and should be considered in the ocean color (OC) data uncertainty budget.


I. INTRODUCTION
O CEAN color (OC) provides access to marine biological quantities, foremost the concentration of chlorophyll-a [1], a major and universal phytoplankton pigment. Other in-water quantities can be determined by Manuscript  OC remote sensing, e.g., in relation to sediment, dissolved organic matter or indicators of water quality [2]. Ultimately all OC in-water products are derived from the spectrum of reflectance characterizing the water body, expressed as remote sensing reflectance R RS or water-leaving radiance L w . This spectral quantity is recognized as an essential climate variable (ECV) by the Global Climate Observing System [3], as is the concentration of chlorophyll-a (Chl-a).
As with any geophysical product, R RS data need to be accompanied by uncertainty estimates to allow an informed use of these data by the user community. Even though there is growing emphasis on the provision of uncertainties for satellite data [4], [5], this practice is in its infancy in the field of OC, which can be partly explained by the fact that this is particularly challenging: OC remote sensing is affected by a large and complex ensemble of error sources [6], including errors associated with the top-of-atmosphere (TOA) signal and the numerous assumptions and approximations needed to solve the remote sensing problem and find a solution for R RS . This requires a proper metrological treatment of the OC data processing to model sources of uncertainties and their propagation in the various processing steps [7].
The remote sensing reflectance R RS (or equivalently L w ) is obtained from the TOA signal L t collected by the sensor in space through a process termed atmospheric correction (AC). The AC process is usually supported by ancillary (or auxiliary) information that documents the state of the atmosphere and sea surface as seen by the sensor. Appropriate quantities are typically provided by satellite products or reanalysis meteorological data given multiple times per day. These data have their own uncertainties that propagate through the AC algorithm and affect R RS . This impact has been largely overlooked by the OC community, assuming that it was small. However, more solid knowledge is required to construct full uncertainty budgets for OC products. This issue is addressed here by exploiting the existence of ancillary data distributed by independent organizations, including ensemble datasets. In the absence of uncertainty estimates associated with each ancillary datum, it is assumed that variations within ensemble simulations and differences between products from independent organizations can be used as representations of uncertainties for testing purposes (as also discussed in [8]). In practice, the present study relies on the processing by a standard AC of Sea-viewing Wide Field-of-view Sensor (SeaWiFS) imagery at a global scale repeated for various sets of ancillary data to evaluate This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the resulting variations in terms of OC R RS products. After having introduced the datasets and approach adopted for the analysis, global results are presented for various test cases. Then the sensitivity of the AC process to individual ancillary variables is discussed before conclusions are given.

II. DATA AND METHODS
Besides being introduced in the text, the main terms used in this article are listed in Table I.

A. Satellite Data
SeaWiFS level-1 global area coverage (GAC) data were obtained from the Ocean Biology Distributed Active Archive Center (OB.DAAC) of the National Aeronautics and Space Administration (NASA), Greenbelt, MD, USA, for the year 2003. They were processed to level-2 with the standard configuration of the l2gen AC of the SeaWiFS data analysis system (SeaDAS) [9] and remapped on a 12thdegree sinusoidal global grid (∼9 km in grid cell size) to get daily global maps. A specific dataset (subsequently referred to as the pixel dataset) is also assembled by accumulating all valid pixels (i.e., with valid AC outputs) of all level-2 images for four representative days in the year (15th of March, June, September, and December). Even though smaller than the 2003 daily level-three data (but with 3.7 million pixels distributed globally and across four seasons), it serves for pixel-based analyses that need to compare data before any averaging takes place (as happens with level-3 data).

B. Atmospheric Correction Algorithm
The l2gen standard AC algorithm [10], [11] requires a set of ancillary variables, including surface wind speed (WS), mean sea-level pressure (SLP), precipitable water vapor (PW), relative humidity (RH), total ozone concentration ([O 3 ]), and nitrogen dioxide concentration ([NO 2 ]). Both zonal and meridional components of WS are actually distinguished but only the modulus of the wind vector is used in the AC. Data of sea surface temperature and sea ice concentration are also requested but were not considered here as they do not directly impact the AC process. Local variations observed in [NO 2 ] can have a significant impact on the AC [12]- [14] but they are found mostly in coastal regions [15]. In addition, l2gen relies on satellite products for [NO 2 ], which are not available for the entire satellite OC record (and [NO 2 ] is not distributed by CERA-20C). For these reasons, the processing carried out in this study used climatology values for [NO 2 ] as the standard NASA processing. It is, however, acknowledged that a dedicated analysis of the impact of [NO 2 ] focusing on coastal regions should be encouraged.
In this section, the main characteristics of the l2gen AC are recalled to ease the interpretation of the results. First, the radiative transfer problem is expressed by the following where L t is a sum of contributions from L w and the radiance associated with air molecules (Rayleigh radiance, L r ), aerosols (L a ), white caps or sea foam (L f ), and glint (L g ) [10], [11]: where t d,v and T v are the diffuse and direct transmittance in the sensor-viewing direction [16], while absorption by gases (e.g., O 3 and NO 2 ) is written by the transmittance terms t g,v and t g,s for the paths from Sun-to-surface and surfaceto-sensor, respectively. Note that in this equation, L w is not the signal that would be actually measured in the field (as it does not include gaseous absorption occurring along the path Sun-to-surface). In the l2gen code implementation, this is corrected by having t g,s as denominator of L w in (1), so that the l2gen L w output (that can be saved in level-2 files) is comparable to that observed in situ. In (1), there is no term for correction of polarization effects as these are not considered for SeaWiFS, and corrections for out-of-band effects are not made explicit. For further use in bio-optical algorithms, L w goes through a process of normalization (sensu Gordon and Clark [17] and Morel and Gentili [18]) correcting for the path transmittance and bidirectional effects, and providing the normalized water leaving radiance L WN or remote sensing reflectance R RS where θ s and θ v are the solar and satellite zenith angles, respectively, and φ the relative azimuth between viewing and illumination direction; t d,s (λ) is the solar-path diffuse transmittance; d andd are the actual and average Earth-Sun distances, respectively. The correction for the bidirectional reflectance distribution function (BRDF) is written as following [18], Q being the ratio between irradiance and radiance just below the water surface, f the ratio between irradiance reflectance and b b /a (ratio of water backscattering and absorption coefficients), and being a function combining all reflection and refraction effects at the air-sea interface [18], [19]. In the l2gen code, the ratio f /Q is obtained from look-up tables (LUTs) [20] with a dependence on geometry and Chl-a, which means that it ultimately depends on all inputs to the code. is also dependent on geometry, and on WS for its above-to-below irradiance transmittance term, but this latter dependence is only seen for high solar zenith angles, above 60 • [19]. For open ocean waters, where L w in the near-infrared (NIR) can be assumed negligible following a black-pixel assumption [21], L a can be computed knowing the other terms of (1). The ratio of the aerosol reflectance (converted into its single-scattering equivalent) at 765 and 865 nm (noted ) can then be compared with tabulated values [10], [22]. The algorithm relies on 80 aerosol models [23] associated with eight values of RH (RH of 30%, 50%, 70%, 75%, 80%, 85%, 90%, and 95%) and ten values of the fraction (in volume) of the fine aerosol mode (from 0% to 95%). First, the two groups of aerosol models with RH bracketing the value obtained from the ancillary data are identified, and for each group, the two most appropriate models are selected on the basis of [10]. Ultimately, L a can be computed in all bands by weighted averages from the four selected models.
This serves as the first step of an iterative process by which a possible contribution of L w in the NIR is computed through a bio-optical model [24]. The contribution from L w is then subtracted from L t in the NIR bands leading to a new estimate of L a , and the process starts again until convergence (or a maximum number of iterations) is reached. Therefore, for any perturbation affecting (1) at one band, L w will not only be impacted at that band but also potentially at all other bands by this loop between visible and NIR wavelengths.
For the clarity of the discussion, it is worth introducing the various algorithmic modules where the ancillary data intervene.
1) WS: WS intervenes in various components of the algorithm. The radiance L f due to white caps is a cubic function of WS for values between 6.33 and 12 m s −1 [25]: with WS below 6.33 m s −1 , L f is considered negligible, while L f for WS above 12 m s −1 reaches saturation.
Since wind modifies the air-water interface geometry, it also has an effect on the glint radiance L g [26], and on the computation of the Rayleigh radiance L r [27]. As introduced earlier, WS is also found in the calculation of R RS from L w . 2) SLP: By modifying the Rayleigh optical thickness τ r (proportional to SLP), SLP has an impact on the transmittance terms t d,s , t d,v , and T v . The Rayleigh radiance L r also depends on SLP [28]. 3) PW: In the case of SeaWiFS, PW does not directly impact the gaseous transmittance terms. It has however an impact on the out-of-band corrections operating in the conversion between single and multiple scattering applied to aerosol reflectance. 4) RH: As explained earlier, RH is a major input to the algorithm as it modulates the selection of the aerosol models.

5) [O 3 ]:
The concentration of stratospheric ozone has an exponential effect on the gaseous transmittance t g .

C. Ancillary Data
Two independent sources of ancillary data were considered for the analysis. The first set relies on meteorological products (WS, SLP, PW, and RH) from National Centers for Environmental Prediction (NCEP), College Park, MD, USA, Reanalysis-2 data [29], [30] available every six hours on a 1 • rectangular grid, and on satellite-derived [O 3 ] given on a daily basis on a 1.25 • × 1 • grid. For 2003, the latter is provided mainly from the Total Ozone Mapping Spectrometer (TOMS) and secondarily from the TIROS Operational Vertical Sounder (TOVS) for a few days (11) when TOMS data are not available [31]. In the rest of this article, this source of [O 3 ] data will be referred to simply as TOMS.
An alternative source of ancillary data was provided by the CERA-20C dataset from a coupled atmosphere-ocean assimilation system of the European Centre for Medium-Range Weather Forecast (ECMWF), Reading, U.K. To favor the coupling between ocean and atmosphere, surface field observations of pressure and marine wind were assimilated in the system [32]. Interestingly, CERA-20C (simply called CERA thereafter) comes in the form of a ten-member ensemble given every 3 h (i.e., ten realizations are provided for each variable and time step). The CERA-20C ensemble generation accounts for errors in observations assimilated in the reanalysis as well as model errors, and its spread can be seen as a measure of uncertainty even though it is likely an underestimate, possibly by a factor of 2 for some variables [32]. CERA-20C directly provides all ancillary data required for l2gen, except RH that was computed from 2-m surface temperature and 2-m dewpoint temperature [33] (the result can, therefore, be considered equivalent to NCEP RH given at 1000 hPa).
The spread of the CERA ensemble and the differences existing between the CERA data and the NCEP or TOMS data are taken as representative estimates of the uncertainties characterizing such ancillary data (see typical values of uncertainties in ancillary data in [6] and discussion in [8]). To drive l2gen with the CERA data (obtained on a 1 • grid), the latter is selected at the NCEP 6-h time steps for the meteorological variables, and bilinearly interpolated on the TOMS grid and averaged for the day for [O 3 ] (the TOMS data being based on polar-orbiting satellite retrievals, a daily global map is not associated with a fixed time). Practically, NCEP or TOMS data were substituted in the files by the CERA equivalent data for easy ingestion within l2gen, where the values applied for each pixel are obtained by linear interpolation, in space (at the pixel location from the gridded ancillary data) and in time (from the ancillary data just before and after the satellite pass).
Besides the reference case based on the original NCEP and TOMS data, it was then possible to process the same SeaWiFS imagery using each member of an ensemble ancillary dataset, thus creating ensemble level-3 SeaWiFS data. In practice, global level-3 SeaWiFS data for 2003 were created 10 times in association with each ensemble member for the following cases. ] at its TOMS value. 7) ALL case: All five variables are substituted by the ten members of the CERA data. Adopting vocabulary typical of ensemble simulations, the ensemble variations of ancillary data that are input to l2gen will be termed perturbations.
For the cases where not all variables are perturbed, it is admitted that ancillary data are no longer strictly consistent (one variable coming from CERA and the others from NCEP/TOMS). However, the distributions provided by NCEP and CERA data are very similar (not shown) and the scope of the analysis is just to quantify the impact of small variations of each ancillary variable in terms of R RS output.

D. Statistical Analysis
For a given variable x j , average and standard deviation were computed over a year of data with the following expressions: x j (4) where N t is the number of time steps: when applied to ancillary data, N t was the number of 6-h intervals for the meteorological data (4×365) and the number of days for [O 3 ], whereas it was the number of valid daily level-3 values for OC products. While μ yr provides the average level of a quantity, σ yr is a measure of its natural variability. For a given time step j and grid cell, average and standard deviation among the N ens = 10 ensemble members (x i, j ) i=1,N ens (whether input ancillary data or output OC products) were computed as σ j is the ensemble spread and represents a measure of the difference between the ensemble members. In order to obtain average estimates of the ensemble spread σ j over the available time steps j of the year, a quadratic average σ ens was computed as Therefore, σ ens is an estimate of the average ensemble spread for the ancillary data or for the output level-3 SeaWiFS data.
Considering the NCEP/TOMS datasets as a reference x ref (as reflecting the standard processing but without implying a higher quality), the differences between the average of the ensemble cases and the reference case were documented with the root-mean-square (RMS) difference σ dif and the average difference δ dif computed over the N t available time steps Again, these statistics were computed for the input ancillary data or the output OC products. While σ ens represents the average spread associated with the ensemble, σ dif and δ dif document the distance between CERA and NCEP/TOMS data, or between their associated OC outputs.

A. Variations Within the CERA Ensemble
The impact of ensemble perturbations in ancillary data on the OC products is shown by Fig. 1 for the various cases where one variable is varied. These results pertain to the year 2003 but are very representative (similar results are obtained for other years as illustrated by the frequency distributions of the statistical quantities σ ens , σ dif and δ dif shown for the ancillary data in Fig. 1 in the Supplementary Material). Fig. 1 [34], [35] that provide field data for assimilation in the CERA system (that to the contrary does not assimilate satellite data); related faint features can be noticed for R RS (443). In general, σ ens for R RS is small when only PW is perturbed [ Fig. 1(e) and (f)]. There is a general similarity between σ ens of PW and R RS in the tropical region where σ ens is highest for PW; a local σ ens minimum is observed for R RS (443) around the globe slightly north of the Equator, associated with relatively low σ ens values for PW. The highest variations for R RS are observed in the southeast Atlantic, and along the coasts of India and northwest Africa, without a corresponding signature for PW in the latter two cases. These regions are very challenging for the AC process because of clouds and aerosols and are associated with scarce OC data coverage so that a few cases might have a disproportionate impact on the annual average.
When only RH is perturbed, σ ens for R RS (443) shows some similarity with RH's σ ens [ Fig. 1 The average ensemble spread σ ens in the ALL case (when all five ancillary data are perturbed) is given by Fig. 2 for the aerosol optical thickness τ a at 443 nm, Ångström exponent α (characterizing the slope of log(τ a ) between 510 and 865 nm), R RS at 443 and 555 nm. The observed patterns are similar for these four products (as well as for other bands) even though with different relative amplitudes. For instance, regions of highest variations are seen in the Southern Ocean for R RS with usually smaller values in the tropical oceans; to the contrary, high values of σ ens are comparable for τ a at low latitudes and in the Southern Ocean. Relatively low values are found for α in the high northern latitudes. Comparing Fig. 2 with Fig. 1, it can be seen that variations of WS contribute most in the northern Pacific and Atlantic Oceans and in the Southern Ocean, reinforced in that case by variations in SLP, whereas

B. Variations Between CERA and Reference Case
A similar analysis can be conducted about the impact on the OC products resulting from the differences between the CERA data and the reference dataset (NCEP and TOMS). As introduced in Section II-D, two metrics were used to document differences between datasets, the RMS difference σ dif and the average difference δ dif . As documented by Fig. 1 in the Supplementary Material, the differences between ancillary datasets, as quantified by σ dif , are usually higher than the CERA ensemble spread (σ ens ) (the σ dif global median is 1.64 m s −1 , 1.05 hPa, 5.1 kg m −2 , 8% and 20 dobs for WS, SLP, PW, RH, and [O 3 ], respectively). Consequently, the differences observed between OC products when processed with CERA and NCEP/TOMS data are usually higher than differences observed for the ensemble spread, as can be seen by comparing Figs. 2 and 4 (note the change in color scales). For the sake of brevity, this section does not repeat for σ dif the full analysis done on σ ens ; it merely shows σ dif obtained when changing all five ancillary variables (Fig. 4), whereas it focuses on the impact of systematic differences between ancillary datasets on the OC products.
The distribution of σ dif in general bears some resemblance with σ ens but significant differences can be noted beside the generally higher values. For τ a (443), patterns of high σ dif are seen in various northern seas (Okhotsk, Bering, Labrador Seas, and Hudson Bay), in the Mediterranean Sea, and along the west African coasts, while the Southern Ocean does not stand out as for σ ens . For the Ångström exponent α, σ dif is rather more homogeneous than σ ens [ Fig. 4(b)]. The largest σ dif for R RS (443) is seen at high latitudes and in some specific regions like the west African coasts [ Fig. 4(c)]. The impact of WS is predominant in the northern Atlantic and Pacific Oceans (but not close to the coasts) and in the Southern Ocean, while variations in RH have the largest effect elsewhere (not shown). While the median σ dif is 1.93 10 −4 sr −1 when all five variables are changed, it is 1.02 and 1.27 10 −4 sr −1 for the WS and RH cases, respectively. Despite the low resolution of the ancillary data, some local features can be noticed, like a high σ dif associated with the wind regime in the Gulf of Tehuantepec (Pacific coast of Central America) [36]. At 555 nm, variations in WS and RH contribute to σ dif in the same regions (see Differences between CERA and reference datasets have a systematic component at the scale of a year, the effect of which is seen as expressed by δ dif in Fig. 5 for the five cases of individual variables being changed. It is recalled that δ dif is positive (negative) when values associated with a CERA case are higher (lower) than the reference case. When only WS is substituted by CERA data, the resulting differences in R RS (443) follow fairly closely those seen for WS [ Fig. 5(a) and (b)] but there are noticeable exceptions. For instance, areas of positive δ dif for R RS (443) in the North Pacific are not obviously linked to patterns in WS. In the Indian sector of the Southern Ocean, δ dif for WS is positive only far from the Antarctic continent which is not so for R RS (443). The feature observed in σ dif in the Gulf of Tehuantepec is associated with a negative δ dif for both WS and R RS (443).
In the SLP case [ Fig. 5(c) and (d)], δ dif is usually small for R RS (443) with higher values in the Southern Ocean that appear Results in the RH case [ Fig. 5(g) and (h)] show higher values of δ dif . At low and mid-latitudes, δ dif for RH and R RS (443) are mostly negatively related. Positive δ dif for RH in the Indian Ocean, the tropical Pacific (except in its eastern part), or the Mediterranean Sea are associated with negative δ dif for R RS (443); along the Atlantic African shores, δ dif for RH is seen positive in the north and south and negative in the equatorial region, contrary to δ dif for R RS (443). Negative δ dif for RH in the mid-latitude ocean is also usually reversed for R RS (443). This inverse correlation is, however, not respected everywhere as can be seen in the Labrador Sea or the Hudson Bay. As for σ ens , the [O 3 ] case is illustrated at 555 nm [ Fig. 5(j)]: a very clear correspondence is seen in δ dif with a poleward gradient from negative to positive δ dif , reflecting a strong impact of the differences in [O 3 ] that exceed 20 dobs, particularly in the northern basins. Since these differences in [O 3 ] appear high, it is mentioned here that the CERA-20C simulations do not include [O 3 ] data in its assimilation scheme. Varying all meteorological data but not [O 3 ] (MET case) leads to much lower δ dif and a loss of the latitudinal gradient (see Fig. 2  ]. This is expected: for a given τ a at 865 nm, τ a at 443 nm will vary with α. For both aerosol products, the maps obtained for δ dif are very close to those obtained when only RH is changed to CERA values (not shown). Also for both quantities, δ dif at low latitudes is mostly positive (i.e., CERA-based data higher than the reference case) with some exceptions like the eastern Equatorial Atlantic and Pacific regions [ Fig. 6(a) and (b)], following the difference in RH [ Fig. 5(g)]. This relation is no longer obvious for high latitudes. In the open subtropical ocean, α is usually lower than 1 (e.g., [37]), corresponding to the candidate aerosol models of the AC having low fractions of the fine mode [23], and for these models, α tends to increase with RH, which is no longer true for higher fractions of the fine mode.
The signs of δ dif tend to be opposite for R RS (443) and τ a (443) (as variations in the aerosol radiance tends to compensate variations in R RS ) but this is not systematic [ Fig. 6(a)

C. Synthesis of Results
The preceding sections documenting the impact of variations in ancillary data on OC products show distinct patterns but their amplitude appears fairly small. In order to put these results in perspective, frequency distributions of various statistical indicators applied to R RS are given in Fig. 7: annual average μ yr , annual standard deviation σ yr , σ ens and σ dif , allowing the comparison of the latter two indicators with the average level of R RS and its average natural variability (see Fig. 1 in the Supplementary Material for statistics associated with ancillary data). First, the order of magnitude of the standard deviation of R RS amounts to approximately a fifth of its annual average except at 670 nm where σ yr is closer to μ yr [in the global ocean, R RS (670) is usually small and its variations consequently tend to be high in relative terms]. As anticipated earlier, σ dif is generally higher (by a factor ∼ 2) than σ ens . For completeness, it is mentioned that the OC products are stored as integers in level-2 files with a slope of 2 10 −6 sr −1 for R RS and 0.0001 for the aerosol products (τ a or α), which quantifies the noise due to digitization of the data.
Reasoning in orders of magnitude with the global median of the ratios between the various statistical indicators, σ ens is only ∼1.3% of μ yr between 412 and 490 nm, but increases to 2.8% and 7% at 555 and 670 nm, respectively (Table II). Compared with the natural variability σ yr , σ ens amounts to 6%-12%. Looking at σ dif instead, its median ratio with respect to μ yr is 2.5%-4% between 412 and 510 nm, increasing to 7.6% and 16.1% at 555 and 670 nm, respectively. These median ratios are five times larger with respect to σ yr , reaching 33% at 555 nm, a band for which R RS average and variability are low in a large part of the ocean.

IV. PIXEL-BASED ANALYSIS
Section III is based on annual statistics using level-3 data. While these provide general results, they do not allow a finer analysis that is instead possible with the pixel dataset. The first point that can be addressed with the pixel-based dataset is the importance of the selection of aerosol models. All ancillary data have, more or less directly, an impact on the radiance terms at 765 and 865 nm, and therefore, on the selection of the four aerosol models used in the AC for each pixel. The pixel dataset can be used to quantify how often the perturbations in ancillary data within the CERA ensemble result in a change in aerosol model selection (i.e., when at least one of the four models selected for the pixel is changed) for at least one ensemble member. When only WS, SLP, PW, RH, and [O 3 ] are perturbed using the CERA ensemble, a change in aerosol models happens for 32%, 7%, 15%, 84%, and 6% of the pixels, respectively. In all cases, σ ens found for R RS when a change in aerosol models occurs is on average higher than σ ens without any change in model selection: for instance at 443 nm, σ ens is higher by a factor that varies from 1.6 (for the SLP case) to 3.4 (for the RH case). The examination of individual images reveals that high variations can occur for isolated pixels due to a change in aerosol model selection A change in aerosol model selection has an impact on the whole visible spectrum and will tend to create the largest perturbations in R RS . Not surprisingly, the RH case is associated with the highest incidence of change since RH has a direct action on the model selection, and RH is the variable having the largest weight in the ensemble spread observed for R RS (Fig. 3). On average, WS is the second contributor to σ ens and is also associated with a fairly high incidence of the aerosol model change. Even for the other variables, a change in aerosol model selection is not unusual. When all five ancillary data are perturbed, a change in aerosol models happens for 89% of the pixels, indicating that even the small variations in ancillary data observed in the CERA ensemble have a widespread impact on the functioning of the AC.
The pixel dataset also allows a direct, pixel-based, analysis of the relationships between differences in ancillary data and differences in OC products and tentative interpretations of the observed behavior. For this, each member of the CERA ensemble (for the cases where one ancillary variable is perturbed) is selected and compared with the reference dataset, since this offers a larger dynamic range than when comparing two members of the ensemble. Results are illustrated for R RS at 443 nm (and at 555 nm for the [O 3 ] case) by Fig. 8 in the form of scatter plots colored according to the density of points (complete results for all R RS bands and aerosol products are shown in Figs.

3-8 in the Supplementary Material). The
Pearson correlation computed on all (3.7 million) pairs is always significant ( p < 0.01) for all variables (τ a at 443 and 865 nm, α, and R RS ) for all five cases. However, the simple linear regression running on all points does not always appear representative of the average behavior of the data distributions, likely overly reacting to outliers (see Fig. 8). Therefore, differences in ancillary data first were binned and the median difference in the OC product was computed for each bin. Then, a linear regression was performed between the middle values of the bins associated with differences in ancillary data and the resulting median values of the differences in OC products. Only bins contributing to at least 1% of the dataset were considered (the black crosses in Fig. 8). The slopes of linear regression given in Table III Table III can be seen as representative estimates of the sensitivity coefficients ∂ R RS /∂ x (where x is one ancillary variable) of the relationship between ancillary data and OC products. These estimates should, however, be treated with caution as there are sometimes large variations in this relationship (Fig. 8) and results might depend on the considered dataset.

A. Wind Speed
In general, in the WS case, τ a (865) is barely affected, while τ a (443) and α tend to decrease with WS (negative slope in Table III) and R RS tends to increase in all bands [ Fig. 8(a) and Fig. 3

in the Supplementary Material]. Equation (1) can be written as
by assuming the Sun at its mean distance and neglecting the BRDF effects for simplicity. As discussed in Section II-B, WS has a direct influence on L r , L f , and L g , and an indirect impact on L a through the influence on the former radiance terms in the NIR. WS is also found in the above-to-below irradiance transmittance term acting in the conversion between L w and R RS but with a small impact only seen for high solar zenith angles, typically above 60 • [19]. Neglecting the impact on τ a , and therefore, on the transmittance terms, partial differentiation of (11) with respect to WS (simply noted w) leads to  Colors go from blue to red with an increasing density of points, with light gray associated with the 2-D bins of the scatter plots contributing less than 1% of the maximum density. The gray line represents the linear regression performed over all points. Dots and crosses represent a 2-D binning of the distribution (see text), with dots associated with bins with less than 1% of the data. The black line shows the linear regression using the binned data with only bins accounting for more than 1% of the data (crosses).
by noting t d the product t d,s t d,v . The importance of each term can be quantified by approximating ∂Y/∂w by Y/w, where Y is any of the relevant variable (e.g., R RS ), w is the difference in WS between the CERA and NCEP values and Y the resulting difference for Y . Only pixels for which w ≥ 0.5 m s −1 (66% of the total number of pixels) are considered (to avoid too small differences for numerical reasons). Fig. 9 shows the various terms of (12) with associated median R RS /w. First, it is worth noting that the observed distribution of R RS /w is consistent with the slope reported in Table III [compare the first line of Table III with the median values of R RS /w in Fig. 9(a)]. Let us first consider the term [4]. For 73% of the analyzed pixels (w ≥ 0.5 ms −1 ), L g is equal to 0 at all bands, which is explained by the fact that the glint correction only operates in the algorithm for glint conditions retained as moderate [38]. The remaining values are widely spread (mostly positive, leading to [4] being negative) without showing distinct patterns (so that [4] is not included in Fig. 9). For example at 412 nm, [4] has a median of −0.94 10 −4 sr −1 m −1 s when it is not null, which means that the impact of WS through the glint radiance can be significant for individual pixels. The radiance due to white caps L f increases as a cubic function of WS, which might decrease R RS in the radiance budget of (11), everything else being assumed equal. The term [3] is equal to −1/π∂ρ f /∂w where ρ f is the white cap reflectance. For 55% of the analyzed pixels, the term [3] is null at all bands: L f is equal to 0 when WS is lower than 6.33 ms −1 , while for winds higher than 12 ms −1 , the effects of WS do not operate as the calculation of L f saturates above this value. For medium WS, the term [3] remains anyway very small, tailing from 0 toward the negative values as expected [ Fig. 9(d)].
Even though only slightly and not systematically, L r tends to decrease with WS in the blue (favoring an increase in R RS ), a behavior that tends to be reversed for higher wavelengths [27], leading to the distributions observed for the term [1] in Fig. 9(b). These effects are reinforced with high solar zenith angles, so more frequently in high latitude regions (which is consistent with the geographical patterns seen in Fig. 1(b). It is the term [2] related to L a /w that has the largest weight in defining the distribution of R RS /w [compare Fig. 9(a) and (c)]. This action originates in the NIR in three ways related to L f , L g , and L r . In the NIR, where L w is usually very small, an increase of L f with WS tends to decrease L a , a change that is propagated across the visible domain. In addition, the white cap reflectance decreases with wavelength in the NIR [39] so that a higher L f leads to a lower (ratio of aerosol reflectance), again leading to lower L a values in the visible bands. The pixel dataset was also created removing the correction for white cap reflectance, and the results are comparable with those observed in Fig. 8(a) (or Fig. 3 in the Supplementary Material), with a distribution of points closer to the regression line (not shown), indicating that the impact of L f is small, either through the term [3] or through its action on L a . In the same way, variations in L g affect L a in the NIR, and therefore, across the spectrum. Another pixel dataset was generated removing the corrections for white cap reflectance and for the glint, leading to a distribution of points again comparable to the case with both corrections but still closer to the regression line (see Fig. 4 in the Supplementary Material). Therefore, the term [2] mostly comes from the variations of L r in the NIR, generally leading to a decrease of L a across the spectrum [and a term [2] usually positive, Fig. 9(c)].
In conclusion, R RS /w is largely the result of changes of L r at the considered wavelength and through their actions on L a in the NIR, with a decreasing amplitude for higher wavelengths (see slopes in Table III). This primary effect of WS is perturbed for some pixels by the actions of L f and L g (terms [3] and [4]), which are negative at a given visible wavelength but are indirectly compensated through the NIR by their actions on L a . These indirect actions actually seem to be prevalent as the observed sensitivity factor ∂ R RS /∂w is lower when L f and L g are ignored in the AC process. For instance, the slope for R RS (443) associated with the first ensemble member [ Fig. 8(a)] is 2.9 10 −5 sr −1 m −1 s; removing the impact of white caps, this slope becomes 2.5 10 −5 sr −1 m −1 s, and 2.4 10 −5 sr −1 m −1 s when white caps and glint are neglected.

B. Sea-Level Pressure
In the SLP case, there is a clear negative relationship between the OC products (for both aerosol and R RS ) with fewer outliers compared with the other cases, so that the slopes of linear regression computed on all data and on binned data are similar. The amplitudes of changes are, however, small [see Fig. 1(d) and (b) and Fig. 5 in the Supplementary Material], with variations close to the digitization levels for R RS (670) and τ a (865) (2 10 −6 sr −1 and 0.0001, respectively) so that no sensitivity factors are given for these quantities in Table III.
An increase of L r in the NIR due to SLP might lead to a decrease in L a in the NIR, and thus, across the spectrum, favoring an increase in R RS . In addition, the relative increase of L r with SLP changes very little with wavelength, so that for a given SLP increase, L r will increase more at 765 nm than at 865 nm (L r being higher at the former wavelength), so that L a will be relatively more depressed at 765 nm, decreasing and again L a at all bands. This is consistent with the negative regression slopes found for τ a (443) and α (Table III). However, the most direct impact of a change in SLP is a change of the same sign for L r [27], leading to an opposite change in L w and R RS . For a given wavelength and everything else being equal, a higher SLP also leads to lower t d,v and t d,s , and consequently higher R RS [since their product appears in (11)]. These contrasting impacts can be written at a given wavelength in the following way, neglecting BRDF effects and assuming a mean Sun-Earth distance, with partial differentiation with respect to SLP [noted p in (13) [3] . (13) In that case, there is no dependence of L f and L g on SLP. Fig. 10 shows the distribution of the various terms of (13) in the approximate form Y/p for p ≥ 0.5 hPa. Again, the median values obtained for R RS /p are in excellent agreement with the slope of linear regression obtained on the binned data given in Table III. The distributions of R RS /p are fairly narrow with some tails toward more negative values, which is consistent with the distribution of points along the regression line in Fig. 8(b) (or Fig. 5 in the Supplementary Material). In (13), the term [3] (that can be written as +τ r /2AM/ p 0 R RS with τ r the Rayleigh optical thickness, AM the air mass equal to 1/ cos θ s + 1/ cos θ v , and p 0 the standard SLP of 1013.25 hPa) is approximately an order of magnitude lower than the other terms (Fig. 10). The main term is term [1] associated with the direct effect of L r on R RS , that is partly compensated by the term [2] (indirect effect of L r on L a ), leading to the remarkable relationship seen in Fig. 8(b) (and Fig. 5 in the Supplementary Material).

C. Precipitable Water
There is a larger variability in the PW case [ Fig. 8(c) and Fig. 6 in the Supplementary Material], with a generally positive relationship between differences in PW and in R RS . The slope is also slightly positive for τ a at 865 nm, but becomes negative at 443 nm and for α (Table III). These results are consistent with a slight increase of L a at 865 nm, associated with a decrease of , and therefore, of L a across the visible bands, ultimately leading to an increase in L w or R RS (Table III). This action operates predominantly in the NIR bands for which water vapor absorption becomes more significant (e.g., [40]). As anticipated in Section II-B, PW affects out-of-band corrections performed in multipleto-single and single-to-multiple scattering conversions of the aerosol reflectance, and a deeper understanding of the mechanisms at play would require a dedicated study. It is anyhow recalled that the impact of PW (within the range of values tested) is small among the five cases considered in this study.

D. Relative Humidity
The frequency distributions found for the RH case do not show a linear pattern even though the relationship tends to be slightly negative for R RS and positive for the aerosol products [ Table III, Fig. 8(d) and Fig. 7 in the Supplementary Material]. As described in Section II-B, the calculation of the aerosol reflectance in the visible is performed as a weighted average of contributions from (up to 4) aerosol models, and the position of the observed RH with respect to the eight tabulated values (30%, 50%, 70%, 75%, 80%, 85%, 90%, and 95%) defines the weights given to the models associated with values bracketing RH. If RH moves across a tabulated value, the selected aerosol models are changed, which introduces a divergence in the calculations of L a . This happens more often for high RH as the steps in tabulated values become finer (5%), an aspect which is not really relevant as RH over the oceans is mostly above 70% [8]. The ultimate effect on the AC of a change in RH depends, among other things, on the selected aerosol models and their associated fine-mode fraction. For instance, the Ångström exponent α usually increases with RH for low fractions (below 0.3), and therefore, low α, while it decreases for high fractions (above 0.8) and high α. The results of Table III merely suggests that the former scenario happens more often, which is consistent with the fact that most of the open ocean, at least in low-to-mid latitudes, is characterized by α not exceeding 1 [37], [41], [42].

E. Ozone Concentration
In the [O 3 ] case, there is not much impact of the differences in ancillary data for τ a (865) but a larger, positive, slope for α and τ a (443) is seen (Table III,  Neglecting the impact of NO 2 , gaseous transmittance t g (the product of t g,s and t g,v ) for SeaWiFS bands is e −τ O3 .AM with τ O3 the optical thickness due to ozone, i.e., the product of [O 3 ] and the absorption coefficient k O3 . Again, the effect of variations in [O 3 ] on R RS can be seen as the sum of a direct effect through t g , and an indirect effect through the NIR on L a , leading to the following differentiation: which can be rewritten as .
(15) Fig. 11 shows the distributions of these terms in the approximate form R RS /O 3 , the median values of which are comparable to the slope of linear regression obtained on the binned data given in Table III. For a given L t detected by the sensor, a larger [O 3 ] increases the value of the TOA radiance that is input to the AC after correction for gaseous transmittance.
In the visible bands, everything else being assumed equal, this is translated into an increase in R RS in a manner proportional to ozone absorption [ Fig. 11(b)]. In the NIR (where L w is very small), this same effect instead leads to an increase in L a , a change that is relatively higher at 765 nm (O 3 absorption being 2.44 times higher at 765 nm than at 865 nm), both effects that tend to increase L a at all visible bands (as well as α and τ a ) and decrease L w . Therefore, when the ozone absorption is small such as at 412 or 443 nm, the term [2] dominates and the sensitivity factors ∂ R RS /∂O 3 are negative. When ozone absorption increases, the opposite effect is seen, for instance at 555 nm where term [1] becomes an order of magnitude higher than term [2], leading to a largely positive sensitivity factor, 5.5 10 −6 sr −1 dobs −1 .

V. DISCUSSION AND CONCLUSION
This study aimed at quantifying the impact of uncertainties of ancillary data on the outputs of a standard AC algorithm. Of course, the results are applicable only for the considered algorithm (l2gen) and for SeaWiFS, but they are expected to hold in general terms for other multispectral sensors when processed by l2gen since the bands are similar and the physical relationships underpinning the algorithms are identical. They might have a more general validity at least in terms of orders of magnitude and/or signs of the impact. Indeed other AC codes based on various algorithmic approaches rely on LUTs for Rayleigh radiance calculations (for which a similar dependence on atmospheric pressure and WS is usually included), model a contribution from white caps, and/or include a correction for ozone (e.g., [43]- [46]). On the contrary, RH is not of general use in AC algorithms. Some results are also conditioned by the considered wavelengths so that they might be affected for sensors with different band sets. Moreover, dedicated analyses would be needed to estimate how uncertainties in ancillary data, particularly the concentrations of absorbing gases (e.g., O 3 or NO 2 ), might affect the subtle spectral features provided by hyperspectral missions, such as Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) [47].
As explained in Section II, any perturbation occurring in the NIR bands has an impact across the spectrum for all variables (through the calculation of the aerosol radiance), while any perturbation occurring at one visible band can also have an impact at all bands if the black-pixel assumption is not valid and the iterative process is activated (this is true for 443, 555, and 670 nm through the L w NIR correction, and possibly for other bands if used in the calculation of Chl-a). The algorithm contains LUTs (for aerosol and Rayleigh reflectance) and divergence points (associated with the selection of aerosol models), and some radiative contributions are affected by thresholds and saturation, like the calculation of the white cap and moderate glint radiances. Changes in ancillary data have an impact on the AC process in various, sometimes compensating, ways. All these elements make understanding the effects of variations in ancillary data particularly challenging and, in general, impede an analytic derivation of these effects. Based on simplifying assumptions, the development and analysis of (12), (13), or (15) go some way into explaining some of the observed behaviors but similar developments are not easily done for all cases. The present analysis relies on calculations performed globally across a year and therefore covers a representative sample of the large diversity (in the geometry of illumination/observation, water content, atmospheric, and meteorological conditions) surveyed by the satellite imagery so that the net effect of uncertainties in ancillary data could be documented. Focusing on R RS , it can be observed that an increase in WS and PW is usually translated into an increase at all bands, while an increase in SLP tends to reduce R RS . No clear trend can be ascertained when RH is varied. In the case of an increase in [O 3 ], R RS tends to increase for bands above 490 nm (most noticeably at 555 nm) but to decrease at 412-443 nm where ozone absorption is lowest.
Thus, it appears that R RS might respond differently to variations in ancillary data, in time and space but also spectrally. The focus of this study is on the outputs of the AC but, for the sake of discussion, the impact of the variations observed on R RS on the main OC derived product, Chl-a, is here briefly presented in terms of σ ens . Chl-a is the standard NASA product obtained by a combination of a blueto-green maximum-band-ratio algorithm and a band-difference algorithm (for oligotrophic regions) [48], [49]. The resulting map of σ ens [ Fig. 12(a)] does not obviously compare with the variations observed for R RS [Fig. 2(c) and (d)], while it is more like a Chl-a average map with high values at high latitudes or in upwelling areas. For instance, while there is a local minimum for R RS in the eastern equatorial Pacific cold tongue, there is a local maximum seen for Chl-a extending much further west. Also notable are the patterns of Chl-a's σ ens in the Southern Ocean that are much less intense than for R RS , while it is consistently high in the northern high latitudes. The global median of σ ens is 0.0034 mg m −3 with most values below 0.1 mg m −3 [ Fig. 12(b)], while the median ratio of σ ens and annually averaged Chl-a is 2.2%. The high range of σ ens , largely associated with the high latitude northern hemisphere, is significantly contributed by variations in [O 3 ] [ Fig. 12(b) 3 ] is due to the inverse response observed for R RS in the blue and green bands (Fig. 8), all the more active in high latitudes with larger optical paths. Obviously, these results could be very different for other derived products, such as absorption or backscattering coefficients, and this example on Chl-a suggests that dedicated uncertainty propagation analyses would be needed for each form of bio-optical algorithms.
An important assumption of the study is that the spread within the ensemble data and/or the differences between the CERA and reference (NCEP/TOMS) data (with frequency distributions and median values given in Fig. 1 in the Supplementary Material) are representative of the uncertainties affecting the ancillary products. The former might be an underestimate [32], while the latter is typical of the differences that can exist between ancillary data distributed by various organizations. As the two products cannot be both close to the truth, their differences are related to the uncertainties of at least one of them. It is interesting to notice that the ensemble spread within the CERA data, σ ens , and the RMS difference between the CERA and reference datasets, σ dif , have important patterns in common (e.g., in the Southern Ocean and North Pacific for WS or SLP, in the subtropical region for PW and RH, a latitudinal gradient for [O 3 ] …), which supports the idea that they offer a representative view of the relative geographical distribution of the ancillary data uncertainties, or at least allow a realistic sensitivity analysis on the AC (see [8] for a more complete discussion). These elements, associated with the application of the AC at a global scale for a full year, indeed suggest that this study encompasses a fairly comprehensive description of the influence of the ancillary data uncertainties on the OC output products. This being said, it is stressed that the results need to be considered with caution as the variations observed for the OC products are fully dependent on those associated with the ancillary data (either expressed as σ ens or σ dif ) and that a similar analysis should be repeated with actual uncertainty estimates for the ancillary data. It is also acknowledged that the observed variations do not include contributions associated with the fairly coarse resolution (spatial and temporal) of the ancillary data (typically 1 • and 6-h-to-daily). For instance, WS in coastal regions or in the wake of islands, or even NO 2 in coastal regions [14], show variations that are not well represented at these scales. Therefore, conducting similar sensitivity analyses at regional scales with appropriate datasets would be welcome. It is also recalled that the study focused exclusively on the effect of uncertainties in ancillary data, ignoring uncertainties associated with the way they are used in the algorithms (related to model errors). A typical example is significant uncertainty related to the conversion from WS to white cap coverage and radiance (e.g., [50]) that might be comparable to the uncertainty on WS.
A comment is due to the systematic differences between the CERA and reference datasets. Indeed, a change of ancillary dataset would require an update of the system vicarious calibration (SVC) coefficients [11], [51] that would remove the impact of systematic differences (or bias) between these data. It is, however, stressed that the systematic difference observed between CERA and reference data (δ dif ) accounts for only a part of the observed differences (expressed by σ dif ).
More importantly, the location of the SVC site operating the Marine Optical Buoy (MOBY) close to the Hawaiian archipelago [52] corresponds to fairly low values of δ dif for the meteorological variables and a negative value for [O 3 ] (Fig. 5); using bilinear interpolation on the four grid points surrounding the buoy site, δ dif is +1.1 m s −1 , 0 hPa, −2.3 kg m −2 , 0.9%, and −13 dobs for WS, SLP, PW, RH, and [O 3 ], respectively. Therefore, a correction of the SVC coefficients due to a change in the source of the meteorological data would not correct for much larger differences seen elsewhere, e.g., in the Southern Ocean for WS and SLP, or in the tropics for PW and RH. In addition, these larger differences can either be positive or negative, so an update of the SVC coefficients could compensate for these differences only in some regions. This is typically observed for [O 3 ]: R RS products processed with the CERA ozone data with updated SVC coefficients would be brought closer to the standard R RS products obtained with TOMS data at the latitudes of MOBY, but their differences would further diverge at high latitudes. Actually, what is known about the patterns of uncertainties affecting ancillary data should be added as auxiliary information when documenting the location of SVC sites [53].
The distributions obtained for statistical parameters such as σ ens , σ dif and δ dif for a full-year show distinct spatial patterns that consequently are often found in the OC output products. Considering the range of perturbations applied, WS and RH appear as the ancillary variables having the largest impact on the OC variables such τ a and R RS , with [O 3 ] playing a role for bands with high absorption, such as 555 nm. In the variance space, the sum of variations observed in R RS when only one ancillary variable is perturbed is close to the variations seen when all five variables are perturbed, which suggests that at first order, the effects are fairly uncorrelated. In general, the variations observed for R RS are fairly low, representing a few % of the average signal, but they are by no means negligible; for instance, the median ratio of σ ens and the annual average μ yr is 2.8% at 555 nm, while the median ratio of σ dif and μ yr is 7.6% for the same band. When comparing with the natural variability expressed by σ yr , the median ratio is in the interval 6.3%-12% for σ ens , and 13.6%-32.9% for σ dif .
In terms of recommendations, including the uncertainties associated with ancillary data when building uncertainty budgets for AC algorithms is certainly required, and for this, having uncertainty estimates for each ancillary datum would be useful. The distribution of representative ensemble datasets would already represent a noteworthy step in the right direction.

ACKNOWLEDGMENT
The Ocean Biology Distributed Active Archive Center (OB.DAAC) of the National Aeronautics and Space Administration (NASA), Greenbelt, MD, USA, is acknowledged for the distribution of the SeaWiFS level-1 data and ancillary files for NCEP and TOMS/TOVS data. The European Centre for Medium-Range Weather Forecast (ECMWF) is thanked for the availability of the CERA-20C ensemble data that allowed this research, and Dr. Bill Bell is thanked for his advice and support on using these data.