Approach for Propagating Radiometric Data Uncertainties Through NASA Ocean Color Algorithms

Spectroradiometric satellite observations of the ocean are commonly referred to as “ocean color” remote sensing. NASA has continuously collected, processed, and distributed ocean color datasets since the launch of the Sea-viewing Wide-field-of-view Sensor (SeaWiFS) in 1997. While numerous ocean color algorithms have been developed in the past two decades that derive geophysical data products from sensor-observed radiometry, few papers have clearly demonstrated how to estimate measurement uncertainty in derived data products. As the uptake of ocean color data products continues to grow with the launch of new and advanced sensors, it is critical that pixel-by-pixel data product uncertainties are estimated during routine data processing. Knowledge of uncertainties can be used when studying long-term climate records, or to assist in the development and performance appraisal of bio-optical algorithms. In this methods paper we provide a comprehensive overview of how to formulate first-order first-moment (FOFM) calculus for propagating radiometric uncertainties through a selection of bio-optical models. We demonstrate FOFM uncertainty formulations for the following NASA ocean color data products: chlorophyll-a pigment concentration (Chl), the diffuse attenuation coefficient at 490 nm (Kd,490), particulate organic carbon (POC), normalized fluorescent line height (nflh), and inherent optical properties (IOPs). Using a quality-controlled in situ hyperspectral remote sensing reflectance (Rrs,i) dataset, we show how computationally inexpensive, yet algebraically complex, FOFM calculations may be evaluated for correctness using the more computationally expensive Monte Carlo approach. We compare bio-optical product uncertainties derived using our test Rrs dataset assuming spectrally-flat, uncorrelated relative uncertainties of 1, 5, and 10%. We also consider spectrally dependent, uncorrelated relative uncertainties in Rrs. The importance of considering spectral covariances in Rrs, where practicable, in the FOFM methodology is highlighted with an example SeaWiFS image. We also present a brief case study of two POC algorithms to illustrate how FOFM formulations may be used to construct measurement uncertainty budgets for ecologically-relevant data products. Such knowledge, even if rudimentary, may provide useful information to end-users when selecting data products or when developing their own algorithms.


A.1 Band ratio Chl model
ChlBR is returned as the solution when ChlBR > 0.2 mg m -3 and is computed as follows: [A1] which has the derivative: where α is a polynomial function. The order of the polynomial, N=4, and the coefficients ai are sensor dependent. Specifically, α is expressed as: [A3] and has the derivative (assuming the coefficients, ai, have no uncertainties): . [A4] The log-ratio, LR, term is: where Rrs,b and Rrs,g are remote sensing reflectances centered on blue and green sensor bands, respectively.

Appendix B: Diffuse attenuation coefficient and uncertainty
NASA's standard algorithm for the deriving diffuse attenuation coefficient at 490 nm, Kd,490 (m -1 ), is based on blue-to-green reflectance ratios (Mueller, 2000). The algorithm was empirically developed using a high quality in situ dataset of coincident Kd(490) and Rrs(λ) data (Mueller, 2000;Werdell & Bailey, 2005) and is computed as follows: which has the derivative with respect to c of: where χ is a polynomial function. The order of the polynomial, N=4, and the coefficients bi are sensor dependent. Specifically, χ is expressed as: and has the derivative (assuming the coefficients, bi, have no uncertainties): . [b4] In this study, we consider the Kd,490 algorithm tuned for SeaWiFS such that the blue and green remote sensing reflectances, Rrs,b and Rrs,g, that are centered on 490 nm and 555 nm, respectively. Also, the fourth order polynomial coefficients b0, b1, b2, b3, and b4 are -0.8515, -1.8263, 1.8714, -2.4414, and -1.0690, respectively. Thus, the log-ratio, LR, term is: The partial derivatives of Eq. B5 are: and . [B6b] The variance in Kd,490 can thus be estimated as: where u(Rrs,490) and u(Rrs,555) are the variances of Rrs,490 and Rrs,555,respectively and u(Rrs,490,Rrs,555) is the error covariance of Rrs,490 and Rrs,555. The partial derivatives in Eq. B8 are computed as: and . [B9b]

References
Mueller, J. L. (2000). SeaWiFS algorithm for the diffuse attenuation coefficient, K (490) Where, apoc and bpoc are constants with values of 203.2 and -1.034, respectively. The BR term is a blue-green reflectance ratio with the numerator and denominator being remote sensing reflectances centered on 443 and 555 nm, respectively. [C2] The derivative of Equation C1 with respect to BR is: , and Equation C2 had the following partial derivatives: , and . [C5] The variance of POC is estimated as: where u 2 (Rrs,443) and u 2 (Rrs,555) are variances of Rrs,443 and Rrs,555,respectively. u(Rrs,443,Rrs,555) is the error covariance of Rrs,443 and Rrs,555. The partial derivatives in Eq. C5 are computed as: and . [

Appendix D: Normalized fluorescent line height and uncertainty
NASA's algorithm for normalized fluorescence line height, nflh (mW cm -2 µm -1 sr -1 ), is a measurement of chlorophyll fluorescence emission under natural sunlight (Behrenfeld et al., 2009). The algorithm uses spectral values of normalized water leaving radiances, nLw. Values of nflh are calculated as the difference between the observed nLw678 and a linearly interpolated nLw678 from two adjacent bands (nLw667 and nLw748). Currently, the algorithm is implemented for MODIS only as: . [D1] We note that nLwi is related to Rrs,i as follows: , where, F0,i is the spectral extraterrestrial solar irradiance (Thuillier et al., 2003).
The variance in nflh is estimated as (assuming that F0 has no uncertainties): [D3] where and .
[D4c] The Generalized Inherent Optical Properties (GIOP) is a semianalytical algorithm used to derive standard IOP data products as distributed by NASA's OB.DAAC. Comprehensive discussion of the GIOP can be found elsewhere (Franz & Werdell, 2010;McKinna et al., 2016;Werdell et al., 2013), however, below we briefly overview the algorithm.

E.1 The forward model
At the core of the GIOP is a forward reflectance model that simulates the spectral sub-surface remote-sensing reflectance, rrs,i, as a function of the water-column's inherent optical properties (IOPs). The default configuration of the GIOP uses the quasi-single scattering approximation of Gordon et al. (1988) to model the subsurface spectral remote-sensing reflectance, , as a function of IOPs: , where, ai is the total spectral absorption coefficient, bb,i is the total spectral backscattering coefficient, and g0 and g1 are constants with default values of 0.0949 and 0.0794, respectively. The coefficient ai can be expressed as the sum of absorbing constituent matter present: where, the aw,i is the spectral absorption coefficient of pure water. The two remaining spectral absorption coefficient terms on the right-hand side of Equation E3 are expressed as a product of a normalized spectral absorption coefficient ( ) and its magnitude (x). The subscripts f, and dg denote the constituents phytoplankton and colored dissolved and detrital matter, respectively. Similarly, bb,i can be expressed as: where the subcomponents of water and particulate matter are denoted by the subscripts w and p. Because pure water IOPs and the spectral shapes of other constituent matter can be parameterized at runtime, becomes a function of three free variables: . [E5] A mathematical solution method (default: non-linear least squares optimization) is then employed to find the optimal set of xφ, xdg, and xp such that best matches the sensorobserved sub-surface remote-sensing reflectance, . A "best match" is achieved once some distance metric (e.g. chi-squared) falls below a predefined threshold. We note that is computed from above-water remote sensing reflectance, according to Lee et al. (Lee et al., 2002): Importantly, the GIOP's structure can be varied at runtime to assign the forward reflectance model, the normalized shapes of the IOP subcomponents ( , and ), and the mathematical solution method. We note that the spectral shape coefficients are normalized at 443 nm.

E.2 Bio-optical models
In the GIOP, the normalized shape components , and are parameterized on a perpixel basis using bio-optical models. Below we briefly describe the bio-optical models used in the default configuration of the GIOP.
The spectral shape is modeled per-pixel using the methodology of Bricaud et al. (1998).

Specifically, is a function of Chl as derived in Appendix A and the spectral vectors Ai and
Bi: , where, the scaling coefficient is: . [E8] The subscript 440 denotes that the scaling coefficient is computed at or near 440 nm. The resulting is chlorophyll-specific, hence the scaling factor xϕ has the physical value of chlorophyll concentration (Chlgiop). The spectral shape is modeled using an exponential function of the form: where, Sdg is treated as a constant with a default value of 0.0183 sr -1 (Werdell et al., 2013). The normalized spectral shape has a value of 1.0 at or near 440 nm. Accordingly, the scaling factor xdg is equivalent to adg,440.
The normalized particulate backscattering coefficient, , is modeled per-pixel using a power law: . [E10] The normalized spectral shape has a value of 1.0 at or near 440 nm. The power law exponent, , is calculated following Lee et al. (2002): where r rs,i obs = R rs,i 0.52 + 1.7R rs,i and rrs,550 is centered on or near 550 nm.

E.3 Inverse solution method
The GIOP has a number of built-in inverse solution (spectral matching) methods that an enduser can select at runtime, these including: Levenberg-Marquardt (LM) optimization, Nelder-Mead (amoeba) optimization, and linear matrix inversion (LMI). In this study, we have chosen to use a LM solution method as is it is the current default in NASA's implementation of GIOP. The cost function is a Chi-squared sum of squares metric. When computing the Chi-squared metric, the following Rrs,i bands were included: 412,425,443,460,475,490,510,532,555,583,617,640,655,665 nm.

E.4 Uncertainty propagation
For non-linear least squares, the variance-covariance matrix of derived best-fit parameters, Vx, can be estimated using the Jacobian matrix, J, and the variance-covariance matrix of model inputs, Vrrs as: . [E13] For the GIOP, the uncertainties of xf, xdg, and xp can thus be estimated as the square root of the diagonal elements of Vx. The matrix is the variance-covariance matrix of rrs. The matrix J is computed as: where the i th wavelength element of each column can be expressed as: and .
[E15c] ν = −0.9 r rs,440 r rs,550 The diagonal elements of are equal to the square of uncertainties in sensor-observed subsurface remote sensing reflectances, , and the off-diagonal elements of , , are the covariances between and . The elements are computed as: , and the off-diagonal can be computed as: where, u (Rrs,i, Rrs,j) is the covariance (if known) of above-water remote sensing reflectances Rrs,i and Rrs,j. The partial derivative term in E16 is: . [E18] We note that for the approximation of u(xf), u(xdg) and u(xp), computed from Eq. E13, the spectral shapes coefficients , , and are treated as uncertainty free. In practice, however, we note that and have inherent uncertainties due to their dependence of Chl and Rrs,i, respectively. We remind the reader that in this study we have focused on propagation of data (radiometric) uncertainties and not on the impact of model (spectral shape) uncertainties. We have nonetheless still included relevant FOFM uncertainty formulations for spectral shape models in the following section.

E.5 Uncertainty in spectral shapes
For , variance of the spectral shape, , is driven by variance in the power law exponent, u 2 (g), which can be computed as: . [E19] The partial derivatives in Eq. E19 have the following form: where, V r rs u 2 (r rs,i ) V r rs u(r rs,i ,r rs, j ) r rs,i r rs, j u 2 (r rs,i ) and . [E23] Finally, is computed as: and spectral uncertainty in the derived backscattering coefficient, , is: where the term is the covariance between xp and . When estimating u (bbp,443) the first and third term in the right-hand side of E25 reduces to zero as has a constant value of 1.0 and thus no variance, hence is zero. We have not estimated/parameterized in this study.
For , uncertainty of the spectral shape, , is driven by u 2 (Chl) and is computed as: where, . [E27] The term in Eq. E27 are expressed as: [E28] and ∂γ ∂ν = −2.4e ν ∂ν ∂r rs,440 = − 0.9 r rs,440 ∂ν ∂r rs,550 = 0.9 r rs,440 (r rs,555 ) 2 [E29] Spectral variance in the derived phytoplankton absorption coefficient, , is: where the term is the covariance between xf and . When estimating u(af,443) the first and third terms in the right-hand side of E30 reduces to zero as has a constant value of 0.055 m 2 mg -1 and thus no variance, hence is zero. We have not estimated/parameterized in this study. We also note the spectral coefficients A, B, and the scaling constant of 0.055 each would have asscociated uncertainties, however, quanitifying these model coefficient ucertainties was beyond the scope of this work.
The spectral variance in the absorption coefficient of colored dissolved and detrital matter, , can be estimated as: . [E31] We note that in the default parametrization of the GIOP, Sdg is treated as a constant, thus there is no variance in , and the first and last terms therefore reduce to zero.
For this study, we have statistically estimated the spectral covariance matrix for SeaWiFS remote sensing reflectances VRrs due to radiometric uncertainty in Lt. The objective was to appraise the impact of spectral covariance terms in the analytical FOFM uncertainty framework, not to exactly quantify the covariances for any given sensor. We accept that it does not encompass a wide variety of viewing geometries, scan angles, water types, and aerosol conditions. In addition, we have not considered the uncertainties in ancillary datasets such ozone concentration, and wind speed. This methodology followed two steps: (i) statistically derive VRrs from a number of SeaWiFS images captured over the South Pacific Gyre, and (ii) numerically estimate the Jacobian matrix, JLt, describing how small changes in Lt affect derived Rrs. The SPG was selected as atmospheric and oceanic gradients in this region can be considered quasi-homogenous at local horizontal scales, thus local variability in Lt,i can, to a first order, be attributed to sensor noise. We have used SeaWiFS MLAC data distributed by NASA's Ocean Biology Distributed Active Archive and Center (OB.DAAC) and processed these using NASA's l2gen code which is distributed as part of the SeaDAS data processing and visualization software package (https://seadas.gsfc.nasa.gov/). Following JCGM (2008) we have applied correction factors, including vicarious calibration and temporal degradation gains, in an attempt to reduce systematic effects.
First, we followed a similar approach to Lamquin et al. (2013) to statistically estimate the covariance matrix of top-of-atmosphere radiances, VLt. A selection of 1,928 SeaWiFS MLAC level-1 files spanning the years 1999 -2010 were processed using l2gen. Each L1 MLAC file was encompassed, in part or whole, a 1° x 1° region centered on 26°S, 122°E in the SPG. The derived level-2 (L2) data products were: Lt,i , Rrs,i and Chl. The quality of each level-2 file was then assessed by examining the proportion of a file flagged as: CLDICE, HIGLINT, or PRODFAIL (probable cloud or ice contamination, sunglint: reflectance exceeds threshold, and failure in any product, respectively). This reduced the number of L2 files to 188. Next for each of the 188 L2 SPG extract files were manually inspected in SeaDAS v7.5, and depending on their size, one-to-three 5x5 pixel subsets were extracted where Chl and Rrs,510 appeared to be quasi-homogenous, resulting in a total of 212 L2 5x5-pixel subset files. Finally, the covariance matrix of Lt was computed for each 5x5 pixel subset from which a median average VLt was generated. Using VLt the correlation matrix RLt was then computed. Figure F1 shows estimated VLt and RLt.

Figure F1: Left-hand side: estimated variance-covariance matrix of SeaWiFS top-ofatmosphere radiances, VLt. Right-hand side: estimated correlation matrix of SeaWiFS topof-atmosphere radiances, RLt.
In this study, we approximated the per band top-of-atmosphere radiance uncertainty u(Lt,i) as 0.5% of total Lt,i: . [F2] Next, we estimate a scaled top-of-atmosphere covariance matrix, V'Lt, on a pixel-by-pixel basis which has diagonal elements of and can be computed as: where the matrix S is a diagonal matrix with elements equal to .
A comprehensive method to propagate radiometric uncertainties through NASA's standard atmospheric correction algorithm was beyond the scope of this study. We instead used a numerical approach to approximate the covariance matrix of sensor-derived remote sensing reflectances, VRrs. We first numerically estimated the Jacobian matrix, J Lt, or partial derivatives of Rrs,i (412 -670 nm) with respect to Lt,i (412 -865 nm) . [F4] To generate JLt, for the SeaWiFS MLAC image captured on 14 March 2007 and centered on the SPG, a subset of the image was extracted (25°S-27°S; 121.3°W -123.3°) with mean solar zenith angle (ssolz) of 24.5° (1-ssolz=0.69°) and mean sensor zenith angle (ssenz) of 23.9° (1-ssenz=1.9°). The subset was then processed using l2gen to derive Rrs,i Next, the scene was reprocessed, however, this time Lt,412 was perturbed by a 0.1% of the average scene-wide value of Lt,412. The perturbation process was performed eight times, once for each spectral band. Finally, the partial derivatives in JLt were approximated numerically using a finite difference method, for example: where, is the derived remote sensing reflectance at 443 derived when Lt,443 was perturbed by a small value DLt,443. Finally, VRrs can be estimated as: .
[F6] u n (L t ,i ) = 0.005 × L t ,i u 2 (L t ,i ) An example VRrs matrix and example of u(Rr,is) spectra are shown in Figure F2. These data were derived for the test SeaWiFS MLAC scene captured on 14 March 2007 over the SPG. All elements of VRrs are positive. We note that the estimated relative uncertainties ( Figure F2) are of similar in shape and magnitude to those reported by Hu et al. (2013) for low-Chl waters.