1. Introduction

Accurately retrieving the spectral diffuse attenuation coefficient of downwelling irradiance ($K_{d}(\lambda )$) from Ocean Color radiometry is essential to quantify both the amount and spectral features of the underwater solar radiation available as a function of depth. The diffuse attenuation coefficient is an apparent optical property (AOP) that depends on the inherent optical properties (IOPs) of the water body, mostly absorption ($a$) and back-scattering ($b_b$) [1] while being moderately affected by external environmental conditions impacting the light field (e.g., solar angle and cloudiness) [2]. Together with surface irradiance, $K_{d}(\lambda )$ is necessary to constrain bio-physical processes such as heating, carbon fixation, photo-chemistry, and is used as an input in many assimilative physical or biogeochemical models [3] as well as primary production models [4].

Because $K_d$ varies as a function of the IOPs of a water body, algorithms have been developed to retrieve $K_{d}$ from remote sensing reflectance ($R_{rs}$), using 1) explicit-empirical fits to in-situ data (Standard Level-2/Level-3 product from NASA/ESA [5,6]), 2) non-parametric techniques such as neural networks (NNs) [7], or 3) semi-analytical approaches that first retrieve IOPs before using them to compute $K_{d}$ [8,9]. All of these algorithms incorporate coefficients obtained through an optimization process. In this process, the coefficients are determined by fitting the algorithm to data collected in the ocean or simulated through radiative transfer calculations. Unfortunately, existing algorithms were constrained and validated with a relatively small number of in-situ data points that do not represent the full variability of the global ocean, or with radiative-transfer model runs using a range of input data assumed to represent the global ocean, but whose distribution does not match the spatial/temporal distribution of optical properties in the whole ocean. In a previous study [10], building on previous work [11], we compiled a novel global database of radiometry data at three different wavelengths and for PAR obtained with sensors onboard BGC-Argo profiling floats and matched them to coincident observations from six different satellite sensors (MODIS-Aqua, MODIS-Terra, VIIRS-SNPP, VIIRS-JPSS, OLCI-S3A, OLCI-S3B). Results showed a significant positive bias of all algorithms evaluated (empirical, NN, and semi-analytical) for low values of $K_{d}(\lambda )$ with $R_{rs}$-derived $K_{d}(\lambda )$, resulting in an underestimation of light at depth. This bias was attributed to the fact that the datasets used to tune the coefficients of the algorithms lacked observations representing the clearest waters of the open ocean where the lowest $K_{d}(\lambda )$ values are found.

Given that the clearest waters of the world represent a substantial portion of the global ocean’s surface area (17.8 % have a $K_d(490)^{NASA/ESA}$ < 0.026 $m^{-1}$ based on MODIS climatology), the objective of this study is to develop a robust algorithm capable of accurately retrieving $K_d$ across the entire range of its oceanic variability. Here, we first set out to recompute (’tune’) coefficients for the semi-analytical $K_{d}(\lambda )^{L05}$ [9,12] algorithm and for the NASA/ESA $K_d(490)^{NASA/ESA}$ [5,6] algorithm, using a more globally representative database of satellite $R_{rs}$ match-ups with $K_{d}$ derived from concurrent BGC-Argo float measurements. Since the data collected in this database are not uniformly distributed across all oceanic regions, we considered the number of data points and their distribution within and across various biogeochemical regions when recomputing the coefficients. Diving deeper into the float-satellite match-up dataset, we find that it does not represent all observed $R_{rs}$ spectra, even in the Mediterranean where most data were collected. Consequently, we elected to assess several analytical algorithms published in the literature. If any of these algorithms demonstrate satisfactory performance with the available data at hand, we hypothesize that they will perform effectively outside the range of this dataset, including in Case-2 waters, as they are less likely to be biased by being tuned to an incomplete dataset.

2. Methods

In this section, we expand on the method used to compile the float dataset, on the subsampling used to make it more representative of global distributions, and on the optimization, based on this dataset, of the coefficients in existing empirical and semi-analytical algorithms. We also list the independent datasets and metrics used to validate the performance of the algorithms with the new optimized coefficients. Lastly, we formulate the new algorithm we are presenting ($K_d^{GF}$).

2.1 Datasets used in this study

2.1.1 BGC-Argo floats match-up dataset

Building on the work of [10,11], we use a dataset consisting of match-ups between Satellite sensors and BGC-Argo floats to evaluate current existing algorithms. There are many advantages to using BGC-Argo floats data as opposed to more classical in-situ data, including extensive (and increasing) spatial coverage [13] that is unbiased temporally, as well as the fact that the radiometers do not suffer from shading such as occurs when deploying them from research vessels (in the least of some skylight). In this study, the BGC-Argo float dataset is used to evaluate existing algorithms and is used to reparameterize coefficients for the new version of said algorithms (which are then evaluated using the in-situ datasets described in Section 2.1.2).

$K_{d}(\lambda )^{float}$Kd(L) computation

$K_{d}(\lambda )$ was computed from the BGC-Argo float dataset following [11]. In brief, float-retrieved $K_d$ (herein $K_d^{float} (\lambda )$) is calculated from downwelling irradiance ($E_{d}(\lambda )$) measurements that are acquired from 0-250m at varying vertical resolution, with the float surfacing around local noon. Since the radiometry is not available at $z=0$, we extrapolated a linear regression of $ln(E_d(\lambda, z))$ with depth, in order to compute $E_d(\lambda,0^-)$ for each profile at each available wavelength. Then, we computed the penetration depth ($z_{pd}$) as $E_d(\lambda,z_{pd}) = Ed(\lambda,0^-)/e$, and used it to obtain $K_{d}(\lambda )^{float} = 1/z_{pd}$.

We compared different techniques for $K_d$ retrieval in [10] and found no significant bias in the accuracy of the retrieval. The method described in [11] was therefore chosen, as it implies that the relative error linked to the radiometer slope factor of $E_{d}(\lambda )$ dominates the uncertainty, as is expected for large irradiance values near the surface (as opposed to the small offset associated with the dark current which may dominate the uncertainties at very low light). Before retrieving $K_d^{float}(\lambda )$, all $E_{d}(\lambda )$ profiles were quality controlled following [14] which removed effects of passing clouds, wave focusing, and other bias-inducing effects. A minimum of 5 $E_{d}(\lambda )$ values in the upper 10m are also required as part of this quality control.

Match-ups between floats and Ocean Color sensors

Match-ups between BGC-Argo float radiometry and satellite $R_{rs}$ are computed as in [10]. We have updated the dataset [15] to include match-ups as recent as March 2023 (compared to April 2022 in the previous work available online [16]). It contains match-ups of $K_d$ (derived from $R_{rs}$) between six different contemporary ocean color satellite sensors and BGC-Argo floats measuring radiometry. Match-up criteria are based on [17]. Match-ups include 2,610 with MODIS-Aqua, 2,757 for MODIS-Terra, 3,426 for VIIRS-SNPP, 1,998 for VIIRS-JPSS, 317 for OLCI-S3A and 232 for OLCI-S3B, for at least one wavelength of $K_{d}(\lambda )$.

Match-ups between floats and MERRA-2

In addition to computing match-ups of BGC-Argo floats with ocean color sensors, match-ups were performed with the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) to retrieve the GMAO MERRA2 meteorological products necessary for the analytical algorithm described below. MERRA-2, an atmospheric reanalysis, is formed by assimilating observations (satellite data) into forecast models, particularly the Goddard Earth Observing System (GEOS) model [18]. This integration enables the creation of global atmospheric datasets that are spatially and temporally continuous, with a notable high temporal resolution of 1 hour. The compiled dataset accessible through OBPG integrates multiple MERRA-2 datasets to offer optimal ancillary data for algorithms, such as atmospheric correction, utilized in the retrieval of ocean color products. This includes essential meteorological and aerosol data. For match-ups with the BGC-Argo float dataset, given the high temporal resolution of MERRA-2, we chose a 1-hour criterion. If more than one image matched up with a float position, the closest image in time was kept. Since MERRA-2 has a coarser resolution than Ocean Color sensors ($0.5^{\circ }{\times }0.625^{\circ }$) match-up data consists only of the data of the pixel closest to the float location (as opposed to a 25$km^2$ box for Ocean Color sensors).

Regional biome-based weighting to compute a representative BGC-Argo dataset

Float coverage and satellite match-ups in the global ocean are unevenly distributed with a substantial portion of the float radiometry data ($\approx 70{\% }$) in the Western Mediterranean, Eastern Mediterranean, and the Southern Ocean (Supplement 1, Table S1). Given that one of the goals of this study is to derive an algorithm that will represent the global ocean, it is important not to create an additional bias by over-fitting empirical coefficients of current algorithms towards specific regions. Given that all of the BGC-Argo floats in this study are in Case-1 open-ocean waters, we elect to use open-ocean biomes. In this study, biomes are defined as regions of ocean basin characterized by a consistent seasonal evolution of sea-surface temperature and chlorophyll-$a$ concentration, based on the reasoning that bio-optical variability is, by design, constrained within those two variables, and thus they represent statistically distinct regions that can be used to fractionate an unevenly distributed dataset such as ours [19].

To perform a biome-weighting necessary to compute the cost function (Section 2.3), two parameters need to be taken into account 1) the surface area of each biome ($Area_i$) - available from [19] for biomes 1:17 to which we added two Mediterranean biomes and 2) the number of match-ups in each biome ($N(Area_i)$) - calculable from the float match-up database. An individual weight ($W_i$) for each specific match-up can then be computed from $W_i = Area_i / N(Area_i)$.

Weights for each region for the full (all sensor combined) float match-up database are listed in Table S1. Biomes 1 (North Pacific Ice), 2 (North Pacific Subtropical Seasonally Stratified), 5 (West Pacific Equatorial), and 17 (Southern Ocean Ice) have less than 15 match-ups (0, 3, 0, and 11 respectively) across all sensors and are therefore not considered in the rest of this study as their very high individual weight would likely bias the new coefficients and because these few match-ups likely have high uncertainties relative to the central tendency in their own region.

2.1.2 Additional independent datasets for validation

We evaluated all algorithms with existing in-situ and synthetic datasets that were used in previous works for global $K_d$ algorithm development. Given that the recomputation of the coefficients for $K_{d}(\lambda )_{New}^{L05}$ and $K_{d}(\lambda )_{New}^{NASA/ESA}$ was performed using the BGC-Argo float dataset, testing the performance of those algorithms on different datasets is crucial. The independent datasets also have statistical distributions of $K_d$ which are significantly different from the float database’s [10] and contain data points in Case-2 waters, something not present in the BGC-Argo float database. Therefore, they provide indications regarding the robustness of an algorithm, such as $K_d^{GF}$ across a wider dynamic range (they are comprised of a higher percentage of high values, and have a strong bias towards specific biomes).

In-situ datasets

The NOMAD (NASA bio-Optical Marine Algorithm Data set) dataset [20] is an in-situ dataset spanning oligotrophic to eutrophic waters and has been used to develop the operational empirical, derived from $R_{rs}$, $K_d(490)$ Level-2/Level-3 product distributed by NASA and contains about 800 simultaneous $E_{d}(\lambda )$ and $R_{rs}(\lambda )$ measurements. The COASTLOOC dataset [21] consists of oligotrophic to eutrophic measurements in European waters and contains 195 pairs of irradiance reflectance below the surface ($R(0^-,\lambda )$) and $K_{d}(\lambda )$. $R(0^-,\lambda )$ was converted to $R_{rs}$ with $R_{rs} = 0.133 \times R(0^-,\lambda )$ following [22].

Synthetic datasets

The International Ocean Color Coordinating Group (IOCCG)’s synthetic dataset, originally designed to develop and validate inversion algorithms [2,23] was also used, where data points simulate natural variability over a wide range of Case-1 and Case-2 waters. One thousand paired $R_{rs}(\lambda )$, $K_{d}(\lambda )$ and associated IOPs such as $a(\lambda )$ and $b_b(\lambda )$ are available from this dataset. Note that the simulations did not include inelastic scattering [9,23] and thus their $R_{rs}(\lambda )$ is not corrected for Raman Scattering for algorithms requiring it. The default sky radiance distribution built into Hydrolight (based on the model of [24]) was used for the IOCCG simulations with a wind speed of 5 $m s^{-1}$. Atmospheric parameters needed to compute $f$ (Eq. (8), associated with the $K_d^{GF}$ algorithm introduced below) were not provided in the original IOCCG dataset, therefore a Hydrolight simulation with the same inputs as the IOCCG was run, and atmospheric parameters were retrieved from it to compute $f$. The computed $f$ for the two different solar angles (30$^{\circ }$ and 60$^{\circ }$) were coupled with the retrieved $a(\lambda )$ and $b_b(\lambda )$ to compute a synthetic $K_{d}(\lambda )^{GF}$ (See Section 2.4).

The IOCCG dataset provides $K_{d}(\lambda )$ at five different depths (0, 5, 10, 20, and 40 meters). In order to compute $K_d$ only in the layer in which there is light, we used a depth weighted $K_{d}(\lambda )$ following the equations found in [9]. This weighting ensures that for wavelengths where light is attenuated quickly, a small weight is given to $K_d$ values where there is a very small amount of light. Those weighted $K_ds$ are used herein as the "Synthetic $K_d$" to which we compare "$R_{rs}$-derived $K_d$".

2.2 Recomputation of algorithms’ coefficients for $K_{d}(\lambda )$Kd(L) retrieval from $R_{rs}$Rrs

Algorithms deriving $K_{d}(\lambda )$ from $R_{rs}$ will be referred to as $K_{d}(\lambda )^{R_{rs}}$. Two of the algorithms evaluated in [10] and their associated coefficients are here re-parameterized.

2.2.1 Empirical algorithms currently used for the operational products

The first algorithm we set out to recompute is the operational Level-2 and Level-3 product for $K_d(490)$ distributed by both NASA (Eq. (1)) and ESA (Eq. (2)). Although NASA and ESA present them as two distinct products, both are based on the same empirically derived formulation between in-situ $K_d(490)$ measurements and blue-to-green $R_{rs}$ band ratio [5], originally developed for CZCS with "blue" defined as the wavelength closest to $490$nm and "green" as the wavelength between $547$nm and $565$nm. It is important to note that these algorithms were designed for open ocean Case-1 waters and are known to not perform adequately in Case-2 waters [7].

NASA’s version is computed as follows, with the $A_i$ coefficients tuned to each specific satellite sensor, and $K_w(490) = 0.0166\: m^{-1}$ being the value used for the diffuse attenuation at $490$nm due to seawater.

For the sake of comparison, NASA’s and ESA’s empirical products (applied to their respective sensors, with their coefficients being sensor dependent) will be grouped in the statistical metrics and referred to as $K_d(490)^{NASA/ESA}$.

2.2.2 Lee et al. semi-analytical algorithm based on IOP retrieval: $K_{d}(\lambda )^{L05}$

$K_{d}(\lambda )$ retrieval from the semi-analytical algorithm (herein $K_{d}(\lambda )^{L05}$) was first published in 2005 [8] and further refined in 2013 [9]. It is based on the relationship between $K_{d}(\lambda )$, IOPs, and the solar zenith angle ($\theta$) and is designed to retrieve $K_{d}(\lambda )$ from satellite $R_{rs}$ for both Case-1 and Case-2 waters. It was constrained using simulations using the radiative transfer code Hydrolight (grouped together to form the synthetic IOCCG dataset [23], see Section 2.1.2). The IOPs for absorption ($a(\lambda )$) and backscattering ($b_b(\lambda )$) at any wavelength are retrieved from Satellite $R_{rs}(\lambda )$ (corrected for Raman Scattering) using the Quasi-analytical algorithm (QAA) [25] version 6 with $\eta _w(\lambda ) = {b_{b_w}(\lambda )}/ {b_b(\lambda })$ and the pure water absorption values from [26]:

We recompute the value of the $[Y, m_1, m_2, m_3]$ coefficients using the float/satellite match-up dataset for $K_d(490)$ and keep constant the coefficient multiplying the sun angle ($m_0 = 0.005$ in Eq. (3)). The coefficients are constant as a function of wavelength [9]. We also compute exclusively the $m_2$ coefficient using the full match-up dataset.

2.3 Cost functions used for the algorithms’ coefficients re-computation

New coefficients were computed for the $K_{d}(\lambda )^{R_{rs}}$ algorithms by minimizing the following cost function:

Uncertainty computation

The uncertainty ($U_i$) in the cost function (Eq. (4)) depends on two parameters for a given match-up: The uncertainty in $K_{d}(\lambda )^{float}$ and the uncertainty in $K_{d}(\lambda )^{R_{rs}}$. We assume the uncertainty of the floats’ estimate of $K_d$ ($U_{K_{d}(\lambda )^{float}}$) to be due to both sensor-related uncertainty (a constant) and a relative uncertainty due to environmental variability (e.g., even in an upper ocean with constant IOPs we expect $K_d$ to vary near the surface). We linearly regress the difference between $K_{d}(\lambda )^{float}$ computed from the linear fit method (extrapolating to $E_d(0^-)$) and $K_{d}(\lambda )^{float}$ computed from a least square regression of an exponential fit (as described in [10]) against $K_{d}(\lambda )^{float}$. The intercept of this regression is our constant error whereas the slope multiplied by $K_{d}(\lambda )^{float}$ is the relative error. We obtain $U_{K_{d}(\lambda )^{float}}=max(2 \times 10^{-3}, 5 \times 10^{-2} \times K_{d}(\lambda )^{float})\: m^{-1}$.

Concerning $U_{K_{d}(\lambda )^{R_{rs}}}$, we used the average absolute percentage differences (AAPDs) between $U_{K_{d}(\lambda )^{R_{rs}}}$ and a validation dataset detailed in [9]. They are wavelength and algorithm-dependant and range from 15.3% ($K_d(490)^{NASA/ESA}$) to 26.4% ($K_d(490)^{Lee05}$).

The total uncertainty for a given match-up for the purpose of the computation of the cost function is computed as the root mean square of the uncertainties, $U_i = \sqrt {U_{K_{d}(\lambda )^{float}}^2 +U_{K_{d}(\lambda )^{R_{rs}}}^2} \: m^{-1}$.

2.4 New algorithm design

After realizing the limitation of existing datasets used to optimize the coefficients of existing algorithms (see Section 3.3), we elected to use a simple published analytical algorithm derived from radiative transfer theory that performed well across various existing datasets. Other semi-analytical algorithms with fewer coefficients than $K_d^{L05}$ and with assumptions about the state of the sky are listed and evaluated in the Supplement 1.

The analytical algorithm (whose output is denoted by $K_d^{GF}$, with GF standing for Gordon and Frouin) is based on [27] and [28] and characterized as follows:

Note the distinction between downwelling irradiance in water ($E_d(0^-)$) and downwelling irradiance in air just above the air/water interface ($E_d(0^+)$). Assuming a cloud-less environment, there are two sources of incident downwelling irradiance ($E_d(0^+)$) to the ocean’s surface: direct sunlight ($E_d(0^+, sun)$) and diffuse skylight ($E_d(0^+, sky)$). The radiance distribution is influenced by the magnitude of the two sources and by the solar zenith angle. In the case of a uniform radiance distribution, Gordon suggests the following parametrization of $D_0$ [27]:

An in-depth discussion regarding the accuracy of estimating the $E_d(0^+, sun)/ E_d(0^+, sky)$ ratio using the $T_{dir}/T_{tot}$ ratio, including simulations at the MODIS wavelength can be found in the Supplementary. Computing Transmittances requires several inputs ($\tau _r$, $\tau _{aer}$, $\omega _a$, $a_{aer}$, Table 1) that are all calculated during the atmospheric correction procedure required to derive Level-2 from Level-1 ocean color data. For a given ocean color sensor and a chosen atmospheric correction model, all those inputs are therefore retrievable directly from the aerosol processing pipeline.

Table 1. Inputs needed for the calculation of $D_0$ following Gordon’s characterization [27] with Frouin’s estimation of $f$.

View Table | View all tables in this article

In this paper, we are working with Level-2 Ocean Color data. Thus, we elected to retrieve those inputs from the OBPG ancillary product (GMAO MERRA2) consisting of aggregated datasets from the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2). Definitions, formulas, availabilities, and references for each input are listed in the table below (with $\alpha$ the Angstrom coefficient [30]).

Using Eq. (8–11) to compute $f$, $K_d^{GF}$ is then computed from:

2.5 Statistical performance metrics

As in [10], the following statistical metrics were used to evaluate the performance of $K_{d}(\lambda )$ algorithms. These specific metrics are widely used to validate $K_d^{R_{rs}}$ performance across the literature [7,8,34].

1. Absolute percentage difference (APD) of the log-transformed $K_d^{R_{rs}}$ and $K_d^{float}$ [8]. The APD gives equal weight to over and under estimations of $K_d^{float}$ while covering a range of value larger than one order of magnitude.

2. Root mean squared difference (RMSD) is calculated as follows on the natural log of the $K_d$ values, with $N$ the total number of match-ups:

3. The bias is defined as the median of the ratio of the natural log of the float and remote observations:

4. $R^2$ is the square of Pearson’s correlation coefficient on the natural log of the $K_d$ values.

5. The slope is defined as the slope of a robust (bi-square weighting function) linear fit using the Matlab integrated function $fitlm$ between the natural log of $K_d^{R_{rs}}$ and the natural log of $K_d^{float}$, not allowing for an intercept.

Statistical metrics and regression coefficients are computed on a subset of our dataset representative of the global ocean, using a Monte-Carlo-type subsampling. For each sensor, the total number of match-ups in each biome is computed. Within each biome, a certain number of match-ups are extracted to obtain a subset that is representative of the percentage of the ocean covered by this biome. The number of match-ups selected in each biome is based on the largest possible amount of match-ups in the biome that had the largest discrepancy between its relative surface area and the number of actual match-ups. For example, in the case of MODIS-Terra and $K_d(490)$, biome 16 has the largest discrepancy with only 16 match-ups and a coverage of 10.76% of the ocean. Therefore we selected all 16 match-ups and sub-sampled the other biomes by the proportional surface area listed in Table S1, for a total number of 172 match-ups. This process is repeated 200 times (each time picking random match-ups from each biome other than the one for which we use all the match-ups), and sub-samples are added together to compile an overall "proportional dataset" for each sensor that is proportionally representative of the global ocean. This overall proportional dataset is used to compute the statistical metrics and not for the cost function minimization. When performing Statistical metrics on the full dataset, this subsampling is performed again on the full dataset to account for the fact that different satellite sensors might have different match-up distributions.

3. Results

3.1 Kd(490) retrievals from reparametrized algorithms

For all satellite sensors, the new coefficients resulted in a significant improvement for both types of algorithm ($K_d(490)^{NASA/ESA}_{New}$ and $K_d(490)^{L05}_{New}$), with a bias closer to one and a lower ADP (Table 2). When deciding on the parameters (iterations, threshold, etc $\ldots$) for the cost function, we noticed that in some cases, for MODIS-Terra and VIIRS-JPSS, the cost function never fully converged within the set number of iterations, resulting in the function not reaching its tolerance threshold. This meant that different sets of coefficients independently had very similar statistical performance when used in computing $K_d(490)^{L05}_{New}$ and comparing it to $K_d(490)^{float}$ for those two satellite sensors. We elected to list the coefficients (Table 3) obtained with 10,000 iterations as the cost function did converge with those values. Note that a cost function not reaching convergence is a sign of over-fitting.

Table 2. Statistical metrics for $K_d(490)^{R_{rs}}$ vs $K_d(490)^{float}$ for each of the six studied satellite sensors with the New (from the reparametrization) and Old (from the original publications) coefficients.

View Table | View all tables in this article

$K_d(490)^{L05}_{New}$ performs similarly to $K_d(490)^{NASA/ESA}_{New}$ for all measured statistical metrics. The bias for the low values of $K_d(490) (< 0.026 \: m^{-1}$) is no longer present in $K_d(490)^{L05}_{New}$ although there appears to be a larger error ratio for some of the higher $K_d(490)$ values than with $K_d(490)^{R_{rs}}_{Original}$ (Fig. 1). Following the reparametrization, the number of match-ups with an error ratio (i.e., the normalized difference between $K_d(490)^{R_{rs}}$ and $K_d(490)^{float}$) larger than $25{\% }$ decreased from 23% to 16% for $K_d(490)^{NASA/ESA}$ and from 32% to 16% for $K_d(490)^{L05}$.

 

Fig. 1. Comparison of $K_d^{R_{rs}}$ retrieval for the full match-up dataset (with all sensors grouped) from two different algorithms with the newly computed vs. the original coefficients. (a) $K_d(490)^{R_{rs}}_{NASA/ESA}$ compared to $K_d(490)^{float}$. (b) Error ratio of $K_d(490)^{R_{rs}}_{NASA/ESA}$ to $K_d(490)^{float}$. (c) $K_d(490)^{R_{rs}}_{L05}$ compared to $K_d(490)^{float}$. (d) Error ratio of $K_d(490)^{R_{rs}}_{L05}$ to $K_d(490)^{float}$. New coefficients are listed in Supplement 1, Table S1. The blue vertical line in (b) and (d) is 0.026, the minimum value of $K_d(490)$ in the NOMAD dataset. Note that the coefficients are derived independently for each satellite sensor.

Download Full Size | PDF

Table 3. New coefficients computed for the $K_{d}(\lambda )^{NASA/ESA}$ equation described in Eq. (1) and new coefficients for the $K_{d}(\lambda )^{L05}$ equation described in Eq. (3).

View Table | View all tables in this article

New coefficients for the whole float match-up dataset (treating all satellite $R_{rs}$ as if coming from one single sensor) were also computed. We elected to only recompute the $m_2$ coefficient from Eq. (3), as we observed that several different sets of coefficients could result in a local minimum of the cost function, i.e., an inability to converge to a single independent solution when computing all 4 coefficients. This was also observed with some of the individual satellite sensors. We elected to recompute $m_2$ as opposed to the other coefficients, as recomputing the other coefficients resulted in lower performance on the independent datasets and/or resulted in negative coefficients. The best fit was with $m_2 = 1.2541$ (Table 3). This coefficient is the one used when assessing algorithm performance on the independent datasets and on the full BGC-Argo float dataset.

3.2 Independent validation of the reparameterization

The new $m_2$ coefficient for $K_{d}(\lambda )^{L05}_{New}$ (Table 3) was also validated in two additional ways. First, we used the same float dataset but assessed the performance of the QAA-based algorithm (L05) at a different wavelength: 412 nm. Since the coefficients were re-computed and fitted only using the 490 nm wavelength, $K_d(412)^{L05}_{New}$ retrieval (following Eq. (3)) using those same coefficients is independent. We observe a significant improvement in retrieval for $K_d(412)$, with a smaller bias, a smaller ADP (26%), a smaller RMSD, and a slope closer to 1 (Fig. 2). Indeed, 30% of $K_d(412)^{L05}_{New}$ data points have an error ratio of $> 25{\% }$ compared to $\approx 46{\% }$ for the original parametrization ($K_d(412)^{L05}_{Old}$) coefficients. Since $K_d^{NASA/ESA}$ was exclusively designed for at 490nm it was not computed at 412nm.

 

Fig. 2. Comparison between the reparametrized global $K_{d}(\lambda )^{L05}$ and $K_{d}(\lambda )^{InSitu}$ for the NOMAD, COASTLOOC, and the synthetic IOCCG datasets. (a) $K_d(490)^{L05}$ versus $K_d(490)^{InSitu}$ with the Original (Old) and the New coefficients for the full datasets. (b) $K_d(490)^{L05}$ versus $K_d(490)^{InSitu}$ with the Original (Old) and the New coefficients on Case-2 waters. (c) $Kd(412)^{L05}$ versus $Kd(412)^{float}$ for the Original (orange) and the New (blue) coefficients. (d) Error ratio of $Kd(412)^{L05}$/ $Kd(412)^{float}$ for the Original (orange) and the New (blue) coefficients.

Download Full Size | PDF

Secondly, the global $K_{d}(\lambda )^{L05}$ was evaluated using additional datasets. The NOMAD, COASTLOOC, and the synthetic IOCCG data span conditions from Case-1 to Case-2 waters and have a different statistical distribution than the float match-up database [10]. We elected not to evaluate $K_d(490)_{New}^{NASA/ESA}$ on the independent datasets as it has sensor-specific coefficients that could not be applied. Additionally, it was not designed to work on Case-2 waters, and poor performance is expected regardless of the parameterization.

Performance in $K_d(490)$-retrieval on those three databases both improved and degraded to some extent for the global $K_d(490)^{L05}_{New}$: it resulted in a slightly higher bias and a higher ADP (Fig. 2). It tends to underestimate $K_d(490)$ from the IOCCG dataset, and converge at higher $K_d(490)$ values. However, we see that the overestimation visible for the low values of the NOMAD dataset is resolved by $K_d(490)^{L05}_{New}$. There is little difference in retrieval performance for the reparametrized algorithms when considering the full in-situ datasets or when selecting only data points corresponding to a Case-2 water type from the IOCCG, NOMAD, and COASTLOOC datasets; in both cases, performance is slightly degraded for the ADP. However, in Case-2 waters, the Bias improves. The reparametrization performed better on Case-2 waters, with very minimal change from $K_d(490)^{L05}_{Old}$. Note that in this comparison we have not weighted the match-ups by the area of the oceans they are representing.

3.3 Evaluation of how representative are the $R_{rs}$ match-ups obtained with the floats compared to a standard Level-2 MODIS image of $R_{rs}$

A Level-2 image captured on 2016-07-28 by MODIS-Aqua (part of the float match-up database), was used to evaluate the representativeness of the $R_{rs}$ domain compiled in the float match-up database compared to the variability in $R_{rs}$ within the image (Fig. 3). To do so, we compared the blue/green band ratio (488/547nm) for $R_{rs}$ obtained from MODIS in this Level-2 image with both the blue-green $R_{rs}$ band ratio in the complete float match-up dataset and the blue-green $R_{rs}$ band ratio specific to the MODIS-only match-up dataset.

 

Fig. 3. Level-2 image from 2016-07-28 from MODIS-Aqua of the Mediterranean region. (a) Case-1 (light blue) and Case-2 (dark blue) $R_{rs}(488)/R_{rs}(547)$ from the whole image compared to the $R_{rs}(488)/R_{rs}(547)$ ratios from the whole match-up dataset (yellow) and the match-ups from both MODIS sensors (orange). (b) Map of the area of interest. Light blue points are values of the $R_{rs}(488)/R_{rs}(547)$ ratio identified as Case-1 waters and dark blue are values associated with Case-2 waters. Yellow points represent data points with values of this ratio present in the float/satellite match-up dataset.

Download Full Size | PDF

We specifically chose to do a comparison in the Mediterranean Sea (Fig. 3) where we have the most match-ups between floats and satellites (Supplement 1, Table S1). We find that although the majority of the Level-2 $R_{rs}$ band ratio values are present in the match-up dataset, the whole dataset does not account for the full variability in this single image. If the spectral $R_{rs}$ variability in the Mediterranean Sea is not encompassed in the float dataset (of which >5,000 are match-ups from the two Mediterranean biomes), the spectral variability in $R_{rs}$ is even less likely to be well represented in areas of the global ocean with less or no match-ups.

3.4 Result of the new $K_d^{GF}$ algorithmResult of the new KdGF algorithm

Following the realization that existing datasets do not span the full ocean distribution of $R_{rs}$, we elected to evaluate a series of other published algorithms and settled on one that, without tuning of coefficient, exhibited good performance. Other published semi-analytical algorithms with fewer coefficients than $K_{d}(\lambda )^{L05}$ were evaluated, with the coefficients determined from a fit to the BGC-Argo floats dataset. They performed well on the float dataset but showed a significant decrease in performance when evaluated on the independent IOCCG dataset. We concluded that despite fewer coefficients, all empirical and semi-analytical algorithms will be affected by the limitations in the datasets they were constrained with. We, therefore, decided to evaluate a fully analytical algorithm that does not show the best performance but is likely to be robust in all conditions: $K_d(490)^{GF}$, and whose uncertainties are well constrained. An in-depth discussion of the different semi-analytical algorithms evaluated can be found in the Supplement 1.

Performance of the $K_d(490)^{GF}$ algorithm was compared to $K_d(490)^{NASA/ESA}$ and $K_d(490)^{L05}$ for the original, published coefficients, as they are the ones currently used by the scientific community as well as for the new recomputed coefficients. $K_d(490)^{GF}$ performs similarly to $K_d(490)^{L05}_{Old}$ on the BGC-Argo database (Table 4), albeit with a slightly smaller $R^2$. Note that $K_d(490)^{NASA/ESA}$ is not computed on the Hydrolight simulations or in-situ databases since this algorithm has sensor-specific coefficients. When compared to the re-tuned algorithms, $K_d(490)^{GF}$ performs slightly less well, but does not show a strong bias in performance for low values of $K_d(490)$ (Table 4). Visual analysis shows this algorithm to slightly underestimate large $K_d^{float}$ values (Fig. 4, Fig. 5). $K_d(490)^{GF}$ performs well on Hydrolight simulations (when compared to the weighted synthetic $K_d(490)$) with a small bias that appears to be proportional to $K_d(490)$. Note that while the statistics on the floats were derived on a sub-sample of the dataset to account for the uneven distribution of our data points across the globe (See Section 2.1.1.4 for details about the "proportional dataset"), the statistics on the historical in-situ dataset and the hydrolight simulations were not.

 

Fig. 4. Comparisons between $K_{d}(\lambda )^{R_{rs}}$ and $K_{d}(\lambda )^{inSitu}$ for several independent databases. (a) $K_d(490)^{R_{rs}}$ versus $K_d(490)^{float}$ using the L05 (orange circles), the NASA/ESA (grey diamonds) and the GF algorithms (blue triangles). (b) $K_d(412)^{R_{rs}}$ versus $K_d(412)^{float}$ using the L05 algorithm (orange circles) and the GF algorithm (blue triangles). (c) $K_d(490)^{R_{rs}}$ versus $K_d(490)^{InSitu}$ using the L05 (orange) and the GF algorithm (blue). In-situ databases are the NOMAD (triangles), and COASTLOOC (circles) datasets. (d) $K_d(490)^{R_{rs}}$ versus $K_d(490)^{IOCCG}$ using the L05 (orange circles) and the GF algorithm (blue squares).

Download Full Size | PDF

 

Fig. 5. Scatter plot of $K_d(490)^{GF}$ versus in-situ $K_d(490)^{float}$ colored by spatial density of nearby points. Spatial density is estimated using Kernel density estimate. (a) Spatial density for the whole full dataset with all match-up points for each sensor. (b) Spatial density for the "proportional dataset" composed of 200 iterations of subsampling according to the biome weight used for statistical analysis as described in Section 2.5.

Download Full Size | PDF

Table 4. Statistics of $K_d(490)^{R{rs}}$ vs $K_d(490)^{float}$ for each of the three algorithms (including the new and old parametrization of the NASA/ESA and L05 algorithms) for the combined satellite sensors. The statistics for the BGC-Argo dataset were derived on the "proportional dataset" derived from 200 sub-samplings of the full dataset (proportionally to the biomes’ weight) as described in Section 2.5.

View Table | View all tables in this article

To understand the algorithm performance statistics the underlying dataset needs to be understood. $K_d(490)^{float}$ are centered around the [0.02-0.25 $m^{-1}$] and the [0.025-0.03 $m^{-1}$] bins, and 75% of values are below 0.053 $m^{-1}$ (Fig. 5). Performance in the area of high data point density (80% below 0.056) follows the 1:1 accurately, while some bias is observed for the values outside this range, such as high $K_d(490)^{float}$ values.

4. Discussion

In this work, we derived new coefficients for existing algorithms used by space agencies to compute $K_{d}(\lambda )$ from $R_{rs}$ so that they are less biased globally. We used weights in the fitting so that the size of each oceanic biome and the number of match-ups within said biomes were accounted for. Additionally, we assumed a specific uncertainty for the fit (proportional to $K_d$ except at low values) and minimized a cost function based on absolute difference. This, of course, means that the fit is different from that which would be derived by weighing all match-ups equally, assuming the same uncertainty for all points, and minimizing a cost function based on root-mean-square (which is more affected by outliers).

A bias for small values of $K_d(490)$ had been identified in previous studies, albeit with different cutoff values ($0.03\: m^{-1}$ in [7] and $0.026 \:m^{-1}$ in [10]) and had not been addressed to date. However, since the area where $K_d(490)^{NASA/ESA}_{Old} < 0.026 \: m^{-1}$ represents 18% of the ocean (based on monthly MODIS climatology) it is likely to have a significant impact on the quantification of physical and biological processes such as heating and primary production. While the reparameterized algorithms and $K_d^{GF}$ are less biased at low values, some small bias remains, particularly for $K_d(490)^{GF}$ and $K_d(412)^{L05}_{New}$.

4.1 Improvement in statistical performance of $K_d$Kd from the BGC-Argo dataset

Using the $K_d$ algorithms with new coefficients resulted in significant improvement in the statistical metrics of match-ups between floats and satellites at the global scale, but the improvement varied depending on the algorithm, sensor, and dataset. $K_d(490)^{NASA/ESA}_{New}$ showed the smallest but still significant improvement at the individual sensor level. $K_d(490)^{NASA/ESA}$ was designed to work on clear Case-1 open-ocean waters, it is thus not surprising that it performs relatively well on the float dataset that is exclusively comprised of open-ocean measurements. The fact that there remains a bias (although of decreased amplitude) for low values when using the newly computed best-fit coefficients (Fig. 1) is likely due to two factors: 1) These smallest values have a relatively low weight in the overall cost function and 2) The $NASA/ESA$ algorithm uses only the blue/green band ratio of $R_{rs}$. This ratio might not encompass the variability in $K_d$ that is due to other parameters than the blue/green ratio (which is directly linked to the chlorophyll-$a$ content of a water body [35]) such as variable concentrations of colored dissolved organic matter (CDOM) and non-algal particles, and the effect of solar angle among others. Therefore, there are limitations in the algorithm design itself, rather than in the coefficient recalculation per se. Additionally, the blue/green spectral ratio follows an asymptotic shape, meaning that once it reaches a certain value of $K_d(490)$, the ratio will no longer be influenced by changes in $K_d(490)$ [8]. The limitations of the NASA/ESA empirical algorithms are well documented and it is widely accepted to be only suitable for open ocean waters [7,8,36].

On the other hand, the semi-analytical algorithm generating $K_{d}(\lambda )^{L05}$ was developed to work on a wide variety of waters and at all visible wavelengths, which means that the original range for which it was designed is much larger than for $K_d(490)^{NASA/ESA}$. However, no data used in either the computation of its original coefficients or the validation had values of $K_d(490)$ below $0.026 m^{-1}$ [8–10] which likely caused the bias for small $K_d$s. This bias was mostly resolved here by fitting new coefficients. Although the small value bias appears to be mostly resolved by the use of the new coefficients, the larger values in our dataset ($K_d(490) > 0.1\: m^{-1}$) are more under-estimated than previously. This appears to be due to the fact that they have a lower weight in the final computation (See Fig. 5 to assess the spatial density of points in the proportional dataset).

The limited geographical area covered by the data in NOMAD (Gulf of Mexico, Eastern Coastal US, Pacific gyre, and Mediterranean Sea) and COASTLOOC (Mediterranean Sea, Eastern, and Coastal North Atlantic) means that they can almost be considered as regional datasets. Given that we aim for a global algorithm, it is not surprising that a "regional" algorithm performs better than a global one in the specific region with whose data it was developed. The fact that there is little difference in the retrieval for the full datasets versus for the Case-2 water types using the new coefficients suggests that $K_{d}(\lambda )^{L05}_{New}$ is able to accurately retrieve Case-2 waters even though it was not derived with such data (Fig. 2). However, we found a large variability in the values of the new coefficients computed for each sensor (Table 3). Additionally, depending on the parameters used (such as tolerance and number of iterations), the cost function minimization was sometimes unable to converge and find four coefficients resulting in an optimal minimization of the $K_d(490)^{float}-K_d(490)^{R_{rs}}$ function, meaning that several different sets of coefficients resulted in similar performance. Additionally, small changes in cost function parameters resulted in widely different coefficients for most sensors. We elected to present one such set of coefficients with 10,000 iterations and 5,000 function evaluations for each sensor in Table 3, but the reader should be aware that there is high variability in possible coefficients depending on which parameters are used as input in the cost function.

This leads to the conclusion that there may be too many coefficients tuned in this algorithm, resulting in several possible "best" solutions when optimizing the coefficients. Additionally, several of the coefficients ($Y,m_1,$ and $m_3$) were found to be negative following the tuning. This, however, is not physically sound, as we do not expect $K_d$ to depend negatively on either $a$ or $b_b$, another sign of the inter-dependence of the coefficients of this algorithm. We thus elected, when trying to compute global new coefficients using the full dataset, to only tune the $m_2$ coefficient (keeping the others as listed in [9]), to obtain a stable physically possible solution.

4.2 Performance of the new algorithm

We find $K_d^{GF}$ to be a robust algorithm to obtain $K_d(490)$ and $K_d(412)$ from $R_{rs}$. Not only does it have better performance than current algorithms designed for Case-1 water such as the NASA/ESA algorithm, but it also performs better than current (non-reparameterized) global algorithms such as $K_d^{L05}$ (Table 4, Fig. 4). There is a small bias in performance for high values (underestimation) of $K_{d}(\lambda )$, but this bias is of the same magnitude as the one for small values currently identified in existing algorithms. This bias is smaller at 412nm than at 490nm and could be a result of multiple-scattering, an effect not captured in Gordon’s algorithm (multiple-scattering will tend to increase diffuse light resulting in a larger $K_d$, effectively resulting in a smaller $K_d^{Rrs}$ than the measured $K_d^{float}$). Future work should be conducted to try and parameterize the effect of multiple scattering. In this study, we used MERRA-2 data to compute the atmospheric parameters (Table 1) needed as an input in the $D_0$ computation. Given the coarser resolution of MERRA-2 compared to Ocean color sensors (0.5$^{\circ }$ as opposed to 1km), this probably introduced some noise and reduced the statistical performance of $K_d^{GF}$. In the future, inputs for $D_0$ could be extracted directly from the pipeline for the atmospheric correction of ocean color sensors, used when going from Level-1 to Level-2 data. $K_d(490)^{L05}_{New}$ has slightly better performance than $K_d^{GF}$ on the float dataset, which is not surprising given that they were tuned to achieve the optimal performance for this specific dataset. However, algorithms using empirically tuned coefficients may not perform well outside the range of data with which they were constrained, as we saw was the case for $K_d(490)^{L05}_{Old}$, where the clearest waters were not part of the original dataset and therefore its performance was biased in such waters. On the other hand, $K_d^{GF}$ is not tuned to a dataset. The good overall performance on all datasets it was evaluated with suggests $K_d^{GF}$ is a robust global algorithm that will likely provide unbiased estimates of $K_d$ from space. Given the good performance of $K_d(490)^{GF}$ on the synthetic $K_d$ (Fig. 4) with maximum values of $K_d(490) = 2.47\: m^{-1}$, we expect that $K_d^{GF}$ will perform well also in Case-2 waters.

Wind speed has the potential to impact $D_0$, by affecting surface roughness, with a rougher sea surface increasing $D_0$. In the original characterization of $D_0$, wind speed of up to $20 m s^{-1}$ resulted in $<2{\% }$ error for $\theta < 50^{\circ }$ [27]. For $\theta > 50^{\circ }$, wind speed should be taken into account and $D_0$ adjusted accordingly. However, in this study, we did not notice any decrease in $K_d^{GF}$ performance for matchups with $\theta > 50 ^{\circ }$, potentially because the impact on $D_0$ was smaller than the noise introduced by the coarseness of MERRA-2 products. Future work should be conducted to assess and account for the impact of wind speed on $D_0$.

An important note is that in the float match-up dataset, $R_{rs}$ was measured from satellite sensors. On the other hand, in the IOCCG dataset, $R_{rs}$ was obtained from radiative transfer simulations, whereas for the COASTLOOC, and NOMAD datasets $R_{rs}$ was measured in situ. It is well known that uncertainties and noise differ between satellite measurements, in-situ measurements, and simulations of $R_{rs}$ [37], hence they will likely result in differences in performance observed between different datasets.

A reasonable criticism of $K_d^{GF}$ is that QAA itself, the algorithm we use to obtain IOPs for $K_d^{GF}$, includes empirical tuning that is based on Hydrolight simulations. Its success with local datasets [23] and the similarity of its output to other IOPs inversions [23] suggests that its performance may not be as biased as the $K_d$ algorithm derived from it ($K_d^{L05}$). More work should be done with other IOP retrieval algorithms to test the sensitivity of $K_d^{GF}$ to the algorithm used to obtain the IOPs.

4.3 Spectral performance of algorithms

Given the upcoming NASA PACE mission there is a necessity to develop a robust algorithm for retrieving $K_d$ from $R_{rs}$ from the UV to red wavelengths. Given the availability of in-situ data and current satellite bands, we were restricted to evaluating the algorithm only at 412 and 490nm. We, therefore, use the IOCCG Hydrolight simulation dataset to evaluate the performance of $K_{d}(\lambda )^{GF}$ and $K_{d}(\lambda )^{L05}_{New}$ at three wavelengths (380, 550, 670 nm) spanning from the UV to the red (note that these simulations do not include Raman Scattering).

We observe an offset between $K_{d}(\lambda )^{R_{rs}}$ and the Synthetic $K_{d}(\lambda )^{IOCCG}$ (Fig. 6). It appears to be $\approx$10% for $K_d(380)$, $\approx$10% for $K_d(550,670)^{GF}$, and $\approx$20% for $K_d(550,670)^{L05}$.

 

Fig. 6. $K_{d}(\lambda )^{L05}$ (Orange) and $K_{d}(\lambda )^{GF}$ (Blue) versus simulated $K_{d}(\lambda )^{IOCCG}$ at 3 different wavelengths. (a) 380nm. (b) 550nm. (c) 670nm.

Download Full Size | PDF

The evaluation of $K_{d}(\lambda )^{L05}$ in the original publications [8,9,12] did not extend beyond 550nm, primarily because high-wavelength $K_d$ in the ocean is dominated by the absorption of pure water. Note that the IOCCG dataset provides $K_d$ values only at 0, 5, 10, 20, and 40 meters. Consequently, only the surface value contributes to computing "Synthetic $K_d$" for large wavelengths. For more details on the weighting of $K_{d}(\lambda )^{IOCCG}$ with depth, refer to Section 2.1.2.

An explanation for the differences in the algorithms’ performances is that the spectral behavior of $K_{d}(\lambda )^{L05}$ is only based on IOPs while in the case of $K_d^{GF}$ both IOPs and the change of contribution of sun and skylight contribute to the spectral dependency.

Further research is needed to validate the performance of these two algorithms across the entire spectrum, ranging from UV to near-infrared wavelengths. While data in the UV (380nm) has been collected on Argo floats (that could be matched to PACE), there is no large dataset with green to red $K_d$s. Deployment of hyperspectral in-situ radiometer sensors on a large scale is necessary to validate and further refine UV-red spanning hyperspectral algorithms. Note that Raman scattering is not expected to affect directly $K_d$ significantly within the first optical depth [38].

5. Conclusions

Previous algorithms for $K_{d}(\lambda )$ have significantly overestimated $K_d$ in the clear waters sampled in the float dataset, resulting in an underestimation of the depth to which light penetrates. It is therefore likely that they are biased in oligotrophic waters, particularly in the subtropical gyres. The underlying reason for the bias in attenuation was the lack of training data for algorithms from clear ocean regions and thus efforts should be made to continue to gather data representative of all areas of the oceans [10]. Here we find that ’retuning’ coefficients using the more comprehensive float dataset removes this bias and the new version performs well on the initial validation datasets that have a different distribution than the float’s. Since empirical or semi-analytical algorithms are, by design, tuned to datasets, not having a comprehensive dataset accounting for all the ocean’s diversity limits the waters where they are likely to work well and results in estimates of uncertainties unreliable globally. We further demonstrated that current datasets, including the float dataset we have assembled, are still not representative of the full $R_{rs}$ spectral variability in the global ocean, and therefore a future user might want to exercise caution when applying it to $R_{rs}$ spectra and biomes not comprised within the float dataset. $K_d^{GF}$, an analytical algorithm based on IOPs retrieval from space, exhibits improved performance relative to existing algorithms in all water types for which we have data. Given its robust performance at the 412 and 490nm wavelengths with current satellite-derived $R_{rs}$ and over a large fraction of the ocean’s biome, we expect $K_d^{GF}$ to perform well in a wide variability of water ranging from Case-1 to Case-2 waters. However, further work should be conducted to assess the impact of surface effects (such as wind-driven waves). Priority should be given to evaluating the robustness of both $K_d^{GF}$ and $K_d^{L05}$ algorithms at different wavelengths using in-situ data.

BGC-Argo floats have been proven to be a very valuable tool for the validation of global satellite products. Their maintenance, improvement, and further deployments should be encouraged, as many regions of the globe are still under-sampled (e.g., the equatorial Pacific Ocean, [10]), with the ultimate goal of creating a dataset representative of the whole ocean. Additionally, recent efforts to augment them with hyperspectral radiometers [39] should be supported as these will allow the evaluation of spectral algorithms at all ocean penetrating wavelengths, with particular relevance to PACE.