Evaluation of ocean color remote sensing algorithms for diffuse Evaluation of ocean color remote sensing algorithms for diffuse attenuation coefficients and optical depths with data collected on attenuation coefficients and optical depths with data collected on BGC-Argo floats BGC-Argo floats

: The vertical distribution of irradiance in the ocean is a key input to quantify processes spanning from radiative warming, photosynthesis to photo-oxidation. Here we use a novel dataset of thousands local-noon downwelling irradiance at 490 nm (E d (490)) and photosynthetically available radiation (PAR) proﬁles captured by 103 BGC-Argo ﬂoats spanning three years (from October 2012 to January 2016) in the world’s ocean, to evaluate several published algorithms and satellite products related to di ﬀ use attenuation coe ﬃ cient (K d ). Our results show: (1) MODIS-Aqua K d (490) products derived from a blue-to-green algorithm and two semi-analytical algorithms show good consistency with the ﬂoat-observed values, but the Chla-based one has overestimation in oligotrophic waters; (2) The K d (PAR) model based on the Inherent Optical Properties (IOPs) performs well not only at sea-surface but also at depth, except for the oligotrophic waters where K d (PAR) is underestimated below two penetration depth (2z pd ), due to the model’s assumption of a homogeneous distribution of IOPs in the water column which is not true in most oligotrophic waters with deep chlorophyll-a maxima; (3) In addition, published algorithms for the 1% euphotic-layer depth and the depth of 0.415 mol photons m − 2 d − 1 isolume are evaluated. Algorithms based on Chla generally work well while IOPs-based ones exhibit an overestimation issue in stratiﬁed and oligotrophic waters, due to the underestimation of K d (PAR) at depth.


Introduction
Light from the sun fuels oceanic primary production, heats the upper ocean, and oxidizes chemical compounds such as organic molecules.In order to accurately model these processes, the subsurface light distribution is needed.As light attenuates near-exponentially in water, to describe the subsurface light field, the exponent describing this attenuation needs to be known (referred to as the diffuse attenuation coefficient [1]).This exponent is a function of in-water components themselves and, to a lesser degree, a function of the illumination conditions [2].Additionally, in order to constrain the layer in which certain processes take place, specific light horizons are of interest.For example, for photosynthesis, a euphotic or isolume depth is defined based on a relative light level (e.g., 1% [3]) or an absolute intensity (e.g., 0.415 mol photons m −2 d −1 [4]) below which photosynthesis is assumed as zero.Moreover, accurately quantifying the attenuation of sunlight within the upper ocean is essential for physical and biogeochemical models, affecting the modeled upper-ocean temperature (e.g., [5,6]) and ecosystem dynamics (e.g., [7]).Sea-surface sunlight is globally available from space agencies ( [8]), and a variety of models to describe its attenuation through the water have been devised (e.g., [9,10]).However, lack of data has limited their validation on global scales.
In the past decade, with the rapid progress of the Biogeochemical-Argo (BGC-Argo) float technology, the BGC-Argo dataset has become the largest data source for optical and biogeochemical observations in the global ocean [11,12].Downwelling radiometry is one of the six core BGC variables in the International Argo program [13,14], and has been used to study global bio-optical relationships and anomalies [15,16], used to determine the depths of the euphotic zone and a specific isolume [17,18], to improve the parameterization scheme of sunlight attenuation in a biological model [19], as well as to correct and derive other bio-optical measurements performed on the same float [20][21][22].An additional important application of BGC-Argo is to validate ocean color remote sensing products (e.g., [23][24][25]) over the whole year and on the global scale.
In this study, based on a recently compiled global BGC-Argo dataset (described in [15,26,27]), first, we evaluate several ocean color remote sensing algorithms for the diffuse attenuation coefficient at 490 nm (K d (490)), including two empirical [28,29] and two semi-analytical algorithms [9,30]; then, we evaluate the performance of the algorithm for the diffuse attenuation coefficient (K d (PAR)) of photosynthetically available radiation (PAR) [10], at different optical depths; finally, we evaluate several empirical [29,31] and semi-analytical algorithms [32] for the euphotic depth (z eu ) and the 0.415 mol photons m −2 d −1 isolume depth (z 0.415 ).Given our dataset, we focus only on open ocean algorithms.The need for validation stems from the fact that these algorithms have not been validated globally and over the whole year.

BGC-Argo Data
The globally distributed BGC-Argo dataset used here was published in the SEA scieNtific Open data Edition (SEANOE, www.seanoe.org/data/00383/49388/)and is comprised of more than 5000 local-noon profiles of downwelling irradiance (different numbers for different wavelengths due to quality control filtering [26]) obtained using 103 Satlantic OCR504 radiometers deployed on Argo floats, spanning 3 years from 20 October 2012 to 26 January 2016.Each radiometry data depth profile includes downward irradiance (E d ) at three wavelengths (380, 412, 490 nm, [µW cm −2 nm −1 ]), and the instantaneous photosynthetically available radiation (iPAR, [µmol photons m −2 s −1 ]), which is the downwelling photon flux integrated over a wavelength range spanning from 400 to 700 nm.All radiometry data has been quality controlled [26], removing the points and profiles that are significantly affected by clouds and wave focusing.The ocean basin and trophic environment of each profile have been identified in the dataset (www.seanoe.org/data/00383/49388/data/49825.pdf),however, in this study, we reorganized them into 10 regions (Table 1 and Figure 1).All symbols used in this study are listed in Table 2.We note that no correction for dark offset has been applied in this dataset beyond that obtained from the manufacturer.Small (O(0.03 µW cm −2 nm −1 )) nonzero dark values are observed in data from profiling floats but these values are variable in sign and hence will not introduce a bias in our current work.The attenuation of E d (490) and of PAR with depth is approximated as an exponential decrease, and the layer-averaged diffuse attenuation coefficients, K d (490) z and K d (PAR) z , from surface (0 − ) to any given depth, z, are defined as: The BGC-Argo radiometry data processing flow chart is shown in Figure 2, where the surface E d (490,0 − ) and iPAR(0 − ) are determined first.Since BGC-Argo has no radiometry exactly at z = 0 − , they are computed by extrapolating a linear regression for ln(E d (490,z)) and a second-degree polynomial regression for ln(iPAR(z)) with depth, z, respectively, at the top 10 m of each profile.To ensure the accuracy of extrapolation, we only process the profiles with at least 5 valid values within the top 10 m (which include 4882 E d (490) and 2548 iPAR profiles, as shown in Table 1).This extrapolation method is validated with in-situ cruise data in Appendix A. Flow chart of BGC-Argo radiometry data processing.Note that all derived K d values are layer-averaged ones (the layer from sea surface to depth z, thus, z as the subscript in the symbols, e.g., K d (490) z ), rather than K d at a specific depth.The same procedure could be used to derive the layer-averaged attenuation coefficient to any depth of interest.
We then determine the penetration depth at 490 nm (z pd ) and euphotic depth (z eu ), as the depth where E d (490,z) reaches E d (490,0 − ) × e −1 ~0.37 × E d (490,0 − ) [36] and the depth (denoted by z eu ) where iPAR(z) reaches iPAR(0 − ) × 1% [3], respectively.z pd and z eu , respectively, are determined through an interpolation of E d (490) and iPAR profiles to depths where E d (490) and iPAR reach the calculated E d (490,z pd ) = E d (490, 0 − )/e and iPAR(z eu ) = 0.01 iPAR(0 − ) (Figure 2).Note that, the definition of penetration depth is spectrally dependent (i.e., z pd varies with wavelength), in this study, the wavelength is omitted for simplicity, and all z pd without designated wavelength represents z pd (490).By definition, the remotely-sensed diffuse attenuation coefficient K d (490) is the layer-averaged one from surface to z pd , and it is determined in situ as: Note that Organelli et al. [15] used an approximation to estimate the penetration depth, to derive K d (490) zpd , namely z pd = z eu /4.6.In Appendix B we show that this approximation results in a significant bias in K d (490) zpd in oligotrophic regions.[4], the isolume depth z 0.415 , defined as the depth where the daily PAR reaches 0.415 mol photons m −2 d −1 , has been used as a threshold depth below which light is insufficient to support photosynthesis, and has been regarded as a useful estimate of the euphotic depth [31, [37][38][39].Phytoplankton living at depth are expected to react to absolute light level rather than a relative level [40], and thus, the isolume is more applicable for phytoplankton growth than the depth of 1% of sea surface light intensity (z eu ).Note that recently, Behrenfeld and Boss [41] argued for a lower threshold irradiance value.
Generally, it is difficult to directly measure z 0.415 , as well as any other daily-isolume depth, as it needs the radiometry at both continuous vertical scale (at least a resolution of meters) and continuous temporal scale (at least hourly).It can be estimated using a single observation during a day (on Argo and ships), with the assumption that the layer-averaged K d (PAR) at z 0.415 is constant during the day (not a constant K d (PAR) profile), i.e.,: Equation (4) allows us to determine iPAR(z 0.415 ) based on a single profile of radiometry and PAR sat , the surface daily photosynthetically available radiation product (see Section 2.2.1).α is the transmission of sun light through the air-sea interface, dependent on the latitude and day of year [33].z 0.415 is determined as the depth where iPAR(z) reaches iPAR(z 0.415 ) which is calculated from Equation (4).Using a similar procedure, one can derive the depth of any other isolume.

Satellite-Based K d (490)
We download three available satellite K d (490) products and derive three other literature-based K d (490) for evaluation using the floats' data.They encompass different approaches (empirical vs. semi-analytical) and different sources (MODIS-Aqua vs. GlobColour): (i) K d (490) M-KD2M is an operational MODIS product derived from an empirical algorithm "KD2M" that uses the blue-green reflectance ratio ("B/G") of satellite-measured remote sensing reflectance [28] and is available from NASA (https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/Daily/4km/Kd_490/); (ii) K d (490) GC-M07 is an operational product retrieved from an empirical relationship between K d (490) and Chla (Equation (5), [29]), provided by the GlobColour project [43], using products merged between MODIS and a visible infrared imaging radiometer (VIIRS) for our data period (ftp://ftp.hermes.acri.fr/GLOB/merged/day/); (5) (iii) K d (490) M-M07 is derived from MODIS-retrieved Chla sat and Equation (5) (the same empirical algorithm as K d (490) GC-M07 ); (iv) K d (490) GC-L05a is an operational product retrieved from a semi-analytical IOPs-K d (490) algorithm [9] available from GlobColour (ftp://ftp.hermes.acri.fr/GLOB/merged/day/); (v) K d (490) M-L05a is derived from MODIS-Aqua R rs , based on the same algorithm as K d (490) GC-L05a ; and (vi) K d (490) M-L13 is derived from MODIS-Aqua R rs , based on a similar but updated semi-analytical algorithm [30].Compared to L05a, the L13 algorithm added a new input: the backscattering coefficient of pure sea water (b bw (490)).Here we take it as a constant, 0.001387 m −1 , which corresponds to the sea-surface water with temperature at 20 • C and salinity at 35 psu, based on [44].We have tried to use the Argo-observed temperature and salinity to compute a more accurate b bw (490) following [44], but found no obvious improvement for K d (490) retrieval, when compared to using a constant value.

Satellite-Based K d (PAR) z
Based on radiative-transfer computations with Hydrolight [2], Lee et al. [10] proposed an IOPs-K d (PAR) algorithm, which can be used to estimate the layer-averaged K d (PAR) at any depth, with the inputs of remotely-sensed surface a(490) and b b (490).In this study, we evaluate K d (PAR) M-L05b at six optical depths (from z pd to 6z pd ).
(ii) z 0.415-L07 , is calculated following [32], whose method for z eu , can also be used for z 0.415 .

Satellite-Float Matchup Criteria
For each float profile, the median value of the corresponding satellite data within a 3 × 3 pixel box centered on the profile's surface position was used if at least 5 values in the box were available and that the profile was within 3 hrs of satellite overpass (Figure 3

Statistical Metrics
Three statistical metrics are used to evaluate the differences between products derived from the float database and satellite products: (1) Mean absolute difference (MAD), which represents the absolute differences between the measured (observation) values, and satellite-derived or model-estimated values; (2) Mean absolute percentage difference (MAPD), which represents the relative differences; (3) and Mean percentage difference (MPD), which represents the relative system bias.They are defined as: Here, M i is the in-situ measured value, E i represents the satellite-retrieved or model-estimated value, and n is the number of observations.

Results and Discussion
3.1.Distribution of K d (490) zpd , z eu , and z 0.415 for the BGC-Argo Dataset In the BGC-Argo dataset, K d (490) zpd varies by over an order of magnitude, from 0.02 m −1 to 0.3 m −1 (Figure 4a).Regionally, the median K d (490) zpd is observed to be lowest in subtropical gyres (0.024 m −1 ) and highest in the Black Sea (0.112 m −1 ), with the highest dynamic ranges in the Western Mediterranean (Med.)Sea, Southern Ocean, and subpolar regions (Figure 4b).z eu spans about an order of magnitude in this dataset, from ~25 m in the (North Atlantic) subpolar gyre to ~250 m in subtropical gyres (Figure 4c), with a median value being shallowest (34.0 m) in the Black Sea (Figure 4d) where high CDOM concentration contributes significantly to attenuation [46].z 0.415 has a similar dynamic range and regional distribution as z eu (Figure 4e,f) with median values varying from 23.5 m in the Black Sea to 118.5 m in the subtropical gyres.It exhibits a few very shallow values (<10 m, see Figure 4e) corresponding to very cloudy days (daily PAR above sea surface are very low, e.g., <2 mol photons m −2 d −1 ).In fact, z 0.415 is often shallower than z eu , because the highest daily PAR at sea surface is about 70-80 mol photons m −2 d −1 for a clear summer day [33], and then the corresponding PAR(z eu ) is about 0.7-0.8mol photons m −2 m −1 , slightly higher than 0.415 mol photons m −2 m −1 (as light intensity in water decreases with depth exponentially, about 5.2 m (6.6 m) are needed to decrease PAR from 0.7 (0.8) mol photons m −2 m −1 to 0.415 mol photons m −2 m −1 for a moderate K d (PAR) as 0.1 m −1 ).Thus, for a clear summer day, z 0.415 may be deeper than z eu by about 5 m; but when clouds cause highly depressed daily PAR values, or in other seasons, z eu is deeper than z 0.415 .In the dataset analyzed here, the median z eu is always deeper than the median z 0.415 in all regions (Figure 4d,f), by 3.64 m in subtropical gyres and by 14.7 m in the Southern Ocean.

Assessment of Satellite Algorithms for K d (490)
The performance of classical B/G algorithm (K d (490) M-KD2M ) and semi-analytical algorithm (K d (490) M-L13 ) outperform the others with this dataset (Figure 5a,b, Table 3).The former has the lowest relative differences (MAPD = 14.1%) and no significant mean bias (MPD = −0.3%),but it has a few points overestimated in the high-value range (>0.2 m −1 ); the latter has the lowest absolute differences (MAD = 0.009 m −1 ) and it does not exhibit the overestimation issue at high-values as the B/G algorithm.Note, however, that both algorithms underestimate K d (490) in the Black Sea.It is likely due to high colored dissolved organic matter (CDOM) concentration in this region [16,46], which affects the performance of both band-ratio and semi-analytical algorithms that have been devised using open-ocean data.K d (490) M-L05a (Figure 4c) also performs well in these highly attenuating waters with similar but slightly worst results than K d (490) M-L13 .Chla-based K d (490) M-M07 has the highest MAPD (19.1%) among all 4 MODIS-Aqua products, with an overestimation bias (MPD = 8.7%) mainly in oligotrophic waters (Figure 4d).As for the GlobColour products, K d (490) GC-L05a exhibits an obvious overestimation in all waters (Figure 4e), having high relative and absolute differences (MAPD = 52.4%,MPD = 51.5%).Comparing MODIS-Aqua and GlobColour with the same Chla-based empirical algorithm (M07), we find little difference (Table 3); however, for the semi-analytical algorithm (L05a), K d (490) M-L05a , is much closer to float-observed K d (490) zpd , than K d (490) GC-L05a , suggesting that the overestimation of K d (490) GC-L05a is associated with GlobColour, rather than the model of Lee et al. itself.In addition, we find that the relationship of float-observed K d (490) zpd and Chla sat (Figure 5f) is different than that in [29], having a larger scale factor and a higher exponent (Figure 5f and Equation ( 11)): We applied this relationship to retrieve K d (490) zpd from Chla sat and found it to have similar statistical differences to K d (490) M-KD2M (MAD = 0.013 m −1 , MAPD = 15.1%,MPD = 2.1%).Although our results show that the updated Chla-K d (490) equation has limited improvement on the K d (490) retrieval from satellite, it suggests the need to re-consider the bio-optical relationship with more in-situ data.
In summary, the GlobColour-provided semi-analytical K d (490) product seems to have a retrieval problem while the classical B/G-based KD2M algorithm performs best except in highly turbid waters.The semi-analytical algorithm [30] also performs well in the open ocean, and even better in the high-value range as it is designed for solving both clear and turbid waters.

Assessment of the Satellite Algorithm for K d (PAR)
The IOPs-K d (PAR) algorithm [10] was the only one to date to estimate the layer-averaged K d (PAR) at any depth.It was developed on the basis of Hydrolight simulation to provide a layer-averaged K d (PAR) at any depth from remote sensing.However, the Hydrolight simulation conducted by Lee et al. [10] assumed a homogeneous distribution of IOPs in the water column.Since the surface a(490) and b b (490) only represent the IOPs within the upper mixed layer, changes of IOPs below the mixed layer depth (MLD) are not accounted for in the IOPs-K d (PAR) algorithm.Satellite-retrieved K d (PAR) z-L05b is evaluated with the BGC-Argo dataset from z pd to 6z pd (Figure 6 and Table 4).It performs well at all six optical depths in the high latitude (e.g., Southern Ocean, Subpolar Gyre, and Arctic Ocean), where most water-columns are well-mixed.However, in the stratified waters (including almost all profiles in subtropical gyres, most in the Eastern Mediterranean Sea, and many in the Western Mediterranean Sea), the algorithm increasingly and gradually underestimates K d (PAR) from 2z pd to 4z pd (from 5.7% to 0.8% of relative light intensity for PAR).In subtropical gyres, MPD decreases from 4.4% at z pd to −14.6% at 4z pd (Table 4).Since Chla increases with depth from MLD to the deep chlorophyll maximum (DCM) in oligotrophic waters [47], both a(490) and b b (490) are expected to increase with depth as well.As a consequence, the IOPs-K d (PAR) algorithm inevitably underestimates K d (PAR), even if the QAA algorithm could retrieve sea-surface a(490) and b b (490) correctly.z eu corresponds roughly to 3-4 times z pd (see Appendix B), and is near the DCM in oligotrophic waters [4,47].It follows that the largest underestimation of K d (PAR) is at 4z pd .At 5z pd and 6z pd (correspondingly 0.43% and 0.28% of relative light intensity for PAR), there are fewer valid BGC-Argo observations and the underestimation issue is not clear.In addition, given that below the DCM Chla, a(490) and b b (490) are observed to decrease with depth, the underestimation of K d (PAR) at such depths will weaken.IOPs-based 374(49) 0.011 m −1 (0.004 m −1 ) 12.5%(10.9%)−0.4%(−1.3%)Next, we evaluate the satellite products of two important light-level horizons, the euphotic layer depth z eu , and isolume depth z 0.415 , which represent a relative and an absolute light horizon, respectively.Both z eu and z 0.415 have optical and biological significance, and have been used to estimate primary production [48], retrieve the vertical distributions of Chla [47], and understand the mechanisms on formation and variability of the DCM [17,18].Based on the BGC-Argo dataset, the Chla-based z eu algorithm (z eu-M07 ) performs better, with the lower absolute and percentage differences of 8.0 m and 12.5%, respectively (Figure 7a), than the IOPs-based one (Figure 7b).Regionally, z eu-M07 exhibits a slight underestimation in the subtropical gyres (MPD = −8.4%)and the eastern Mediterranean Sea, and an overestimation in the Black Sea where CDOM deviates significantly from its global relationship to Chla [16,49].The CDOM index (relative CDOM concentration to the global mean CDOM-Chla relationship) is extremely high in the Black Sea and low in oligotrophic gyres [50].Therefore, the Black Sea (Subtropics) is expected to have a shallower (deeper) z eu than the empirical Chla-z eu relationship would provide (High CDOM shoals z eu and low CDOM deepens z eu ).The overall performance of z eu-L07 is still satisfactory (Figure 7b), with MAD of 14.0 m, but with a bias (MPD = 8.4%) which appears in the clearest waters (z eu > 80 m).In subtropical gyres, the overestimation is most remarkable, with MPD reaching 12.0% (Table 4).This issue is likely due to the constant IOPs with depth assumed in this model design which does not include a DCM, consistent with the underestimation of K d (PAR) -L05b at 3z pd and 4z pd .
We assess two satellite z 0.415 algorithms (Figure 7c,d), which to our knowledge has never been done, although the equation (Equation 7) proposed by Boss and Behrenfeld [31] has been used in several subsequent studies (e.g., [51][52][53][54]).Overall, the patterns of two scatter plots of z 0.415 (Figure 7c,d) are similar to those for z eu (Figure 7a,b), with similar statistics (Table 4).z 0.415-B10 displays a slight but obvious underestimation in Subtropical Gyres and Eastern Mediterranean Sea as it is based on z eu-M07 (Figure 7c) and z 0.415-L07 exhibits the same overestimation in subtropical gyres (Figure 7d) as z eu-L07 (Figure 7b).
In summary, the evaluation of satellite-retrieved z eu and z 0.415 suggests that the Chla-based algorithms outperform the IOPs-based ones, with lower scatter and lower relative system bias.Both IOPs-based z eu and z 0.415 exhibit the overestimation in subtropical gyres due to the underestimation of K d (PAR) z-L05b at depth.

Final Remarks and Conclusions
The synergy and joint use of BGC-Argo and satellite remote sensing data contribute to studies using both observing assets [55].BGC-Argo floats provide the largest dataset for validation and evaluation of satellite products in the global ocean, extend the satellite ocean color observations from surface to depth, and fill missing data in satellite coverage due to low sun angle, high latitude winters and clouds; remote sensing is helpful in guiding the deployment of BGC-Argo floats (e.g., in the subtropical gyres or seasonal bloom regions), identifying the spatial scale of float-observed phenomena (basin-, meso-or submeso-scale), extending the BGC-Argo observation from discrete locations to continuous temporal and spatial distributions, and even as a method for calibration of chlorophyll fluorometers deployed on floats (e.g., [23]).
In this study, we present the use of BGC-Argo data to assess existing satellite products.First, for the MODIS-Aqua K d (490) products, both the B/G algorithm and the semi-analytical algorithm [30] perform well, in the open ocean.Moreover, our dataset shows, the B/G algorithm, while having the lowest bias in the open ocean, exhibits an overestimation at large-values region (K d (490) > 0.2 m −1 ), not exhibited by the semi-analytical one.All K d (490) algorithms underestimate in the Black Sea due to extremely high CDOM.The GlobColour K d (490) GC-L05a product has an abnormal overestimation for all the data, which is likely not related to Lee et al.'s model embedded in it.The agreement between float and satellite platforms can be regarded as a "consistency check": on one hand, it provides the validation of existing products, and on the other hand, it suggests that the quality-controlled float-observed E d (490) [26] is of high quality.
K d (PAR) is a critical variable for the retrieval of z eu , and z 0.415 , and is also useful for marine ecosystem modeling (e.g., [56]) and the estimation of biological heat effect on the upper-layer oceans [5].However, K d (PAR) is affected by differences in diffuse attenuation coefficient at different wavelengths.Near the surface, the decrease of PAR is mainly dominated in the open ocean by the losses at the red band due to water absorption [57].With the depth increasing, as red and near-UV light nearly disappear, the green and blue bands make increasing contributions to K d (PAR).Finally, K d (PAR) at deep waters will be close to K d of the blue/green band (440-510 nm), which is most penetrative in the clear waters [5].Owing to its spectral sensitivity, the layer-averaged K d (PAR) z at a certain depth (z) is not a linear average of K d (λ) z spectra, but a weighted average, and the weight varies as function of both downwelling irradiance spectra (E d (z)) and attenuation spectra K d (λ) z [57].This characteristic makes devising an algorithm for it challenging close to the sea surface [58].While, in the deeper waters, the vertical change of IOPs becomes the main error sources of retrieval algorithm of K d (PAR) z .Subsurface structure of IOPs can vary for the same surface conditions, therefore, it is difficult to predict the whole PAR profile accurately, when only relying on the sea-surface IOPs information.Even so, overall, the IOPs-based model [10] performs well near the sea surface until ~2z pd , and also estimates the vertical change of K d (PAR) well in well-mixed waters, but underestimate remarkably below 2z pd in stratified waters (including almost all profiles in Subtropical Gyres, most in the Eastern Mediterranean Sea, and quite a few in the Western Mediterranean Sea), due to the presence of a DCM not accounted for in its design.In turn, affected by the underestimation of K d (PAR) below 2z pd , both IOPs-derived z eu and z 0.415 exhibit a bias in the oligotrophic waters.Chla-based algorithms to estimate z eu and z 0.415 [29,31] perform better although both exhibit some underestimation in subtropical gyres and eastern Mediterranean Sea.
Besides the validation of remote sensing products performed here (and the identification of the significant bias associated with some regions), the statistics associated with our validation can be used to assign errors to the remote sensing products, for example, to propagate them when used in models.For example, many primary-production algorithms parametrize light attenuation using the euphotic depth or include the diffuse attenuation of light attenuation directly (e.g., [48,59]).Simple error propagation (e.g., [60]) can provide the uncertainty in primary production resulting from their uncertainties.

Appendix A. Independent Evaluation of Processing Methods of BGC-Argo Data
In Figure A1, the processing methods of BGC-Argo radiometry data are validated with traditional ship-borne measurements in the Biogeochemistry and Optics South Pacific Experiment (BIOSOPE) cruise [61], which was conducted in the Southeast Pacific during October-December, 2004.Its radiometry dataset had 39 spectral-irradiance profiles recorded by a Satlantic profiler and corresponding surface irradiance measurements recorded by the Satlantic TSRB (Tethered Spectral Radiometer Buoy), spanning various trophic environments, from eutrophic (west of Marquesas Island, and the upwelling conditions off Chile) to ultra-oligotrophic (center of South Pacific subtropical gyre).We first conducted the same quality-control procedure as Organelli et al. [26] on all BIOSOPE radiometry data, removing the noisy profiles and points, then the measured irradiance values above the sea surface (E s ) are converted to the ones just below sea surface E d (0 − ), by multiplying the transmission coefficient α [33], and using a similar procedure to obtain iPAR(0 − ).Following the same steps as Figure 2, z pd , z eu , K d (490), and K d (PAR) are calculated.All these values are denoted as "measured", as they are obtained based on measured E s .We follow the same extrapolation methods mentioned above to estimate E d (490,0 − ) and iPAR(0 − ) from below water measurements.Then, extrapolation-based derived z pd , z eu , K d (490), and K d (PAR) are denoted as "estimated".
After quality control procedures similar to those used for radiometry measured with floats, there remain 23 E d (490) and 21 iPAR profiles from BIOSOPE for validation (due to the requirement of at least 5 samples within the top 10 m).We find no obvious differences between the linear and non-linear extrapolation methods for E d (490), although both have a slight overestimation (MPD = 6.1% and 4.9%) (Figure A1a).Since the attenuation coefficient at 490 nm is mainly determined by IOPs, K d (490) is expected to vary little within the top 10 m (it will vary some due to the adjustment in the mean-cosine at that wavelength and Raman scattering).Thus, ln(E d (490)) is expected to be a nearly linear function of depth.As for iPAR, a second-degree polynomial extrapolation is necessary due to the strong attenuation of red wavelengths especially in the top 10 m just below the sea surface due to strong attenuation of red wavelengths [57], causing the linear extrapolation to yield a significant underestimation (MPD = −15.6%, Figure A1b).Derived z pd , z eu , K d (490) zpd , and K d (PAR) zpd are computed to evaluate the error due to propagation of estimated sea surface E d (490,0 − ) and iPAR(0 − ) (Figure A1c-f).The influence of extrapolation is very limited, although K d (490) zpd and K d (PAR) zpd exhibit a slight bias.
We conclude from this validation exercise that our extrapolation methods to obtain iPAR(0 − ) and E d (490,0 − ) and computing z pd, z eu , K d (490) zpd , and K d (PAR) zpd from radiometers on profiling floats to have well constrained and small uncertainties.Organelli et al. [15] proposed a method to derived K d (490) zpd from BGC-Argo data by first determining z eu from iPAR profile as we do here (Section 2.1), then determining z pd based on the proposed equation: z pd = z eu /4.6 (derived based on the fact that -ln(0.01)≈ 4.6).Finally, they applied a linear regression between ln(E d (490)) and depth to retrieve K d (490) zpd .In contrast, here we use Equation (3) to compute z pd .
When we follow their method to derive K d (490) zpd-O17 , it displays a slight underestimation in oligotrophic waters in comparison to ours (Figure A2a).The relationship z pd = z eu /4.6 is likely biased in most oligotrophic waters.On one hand, the diffuse attenuation of visible light is lowest in the blue band for such waters [62], and hence z eu (the 1% light depth for PAR) must be shallower than z 1%490 (the 1% light depth for E d (490), i.e., z eu < z 1%490 in most waters [57], as shown in Figure A2b).On the other hand, from Equation ( 1 which implies that the ratio of z 1%490 to z pd could reach -ln(0.01)(approximately 4.6) only when K d ( 490) is uniform from surface to z 1%490 (i.e., K d (490) z1%490 = K d (490) zpd ).However, due to the deep chlorophyll-a maximum (DCM) observed in highly clear and stratified waters [47], one would expect K d (490) to increase with depth to the DCM depth, which means z 1%490 /4.6 < z pd (Figure A2c).It follows that z eu /4.6 < z pd in highly stratified (oligotrophic) waters, as shown in Figure A2d, consistent with the observed bias in [15] obtained in such waters (Figure A2a).Statistically, for all valid 2272 samples in our BGC-Argo dataset (Figure A2d), z eu = (3.53± 0.83) × z pd .Regionally, the ratio of (z eu /z pd ) reaches its highest values in the North Atlantic subpolar gyre (4.05 ± 0.83), and its lowest values in the subtropical gyres (3.05 ± 0.46) and New Caledonia (2.98 ± 0.60).We note that Lee et al. [30] also showed that, z 1% of blue light (the arithmetic average of z 1%412 , z 1%443 , z 1%488 , and z 1%531 ) was generally 30%-40% deeper than z eu in open oceans.
Figure 2.Flow chart of BGC-Argo radiometry data processing.Note that all derived K d values are layer-averaged ones (the layer from sea surface to depth z, thus, z as the subscript in the symbols, e.g., K d (490) z ), rather than K d at a specific depth.The same procedure could be used to derive the layer-averaged attenuation coefficient to any depth of interest.

Figure 3 .
Figure 3. Diagram of a good match-up between satellite and float data.The 9 squares represent 9 pixels closest to the location of float surfacing, with 5 valid pixels (blue) and 4 invalid ones (white).The spatial resolution is 4 km.In such a condition or with more than 5 valid-value pixels out of 9, a pair of valid satellite-float matchup data is obtained.

Figure 4 .Figure 5 .
Figure 4. Histogram and boxplot of the distribution of the (a,b) near surface diffuse attenuation coefficient at 490 nm (K d (490) zpd ), (c,d) euphotic layer depth (z eu ), (e,f) and isolume depth (z 0.415 ), for each of the 10 regions (SO: Southern Ocean; AS: Arctic Sea; SPG: Subpolar Gyre; BS: Black Sea; TZ: Transition Zone; WMS: West Med.Sea; RS: Red Sea; EMS: East Med.Sea; NC: New Caledonia; STG: Subtropical Gyre).In the histogram, black and red columns (as well as numbers) represent the full dataset and the satellite-matched ones, respectively.In the boxplot, red points beyond the end of the whiskers represent outliers beyond the 1.5 × IQR (IQR = interquartile range) threshold.

#
The numbers in the parentheses represent the statistics in subtropical gyres.

Table 1 .
Regions and Basins Classified in This Dtudy and Corresponding Profile Numbers of E d (490) and iPAR.

Table 2 .
Symbols used in this study.

Table 3 .
Statistical results of evaluation on the satellite-based K d (490) products *.

Table 4 .
Statistical results of evaluation on the satellite-based K d (PAR) and related products (all products compared here are calculated for the MODIS-Aqua platform).
3.4.Assessment of Satellite Algorithms for z eu and z 0.415