Particulate Backscattering in the Global Ocean: A Comparison of Independent Assessments

How well do we know the particulate backscattering coefficient (bbp) in the global ocean? Satellite lidar bbp has never been validated globally and few studies have compared lidar bbp to bbp derived from reflectances (via ocean color) or in situ observations. Here, we validate lidar bbp with autonomous biogeochemical Argo floats using a decorrelation analysis to identify relevant spatiotemporal matchup scales inspired by geographical variability in the Rossby radius of deformation. We compare lidar, float, and ocean color bbp at the same locations and times to assess performance. Lidar bbp outperforms ocean color, with a median percent error of 18% compared to 24% in the best case and a relative bias of −11% compared to −21%, respectively. Phytoplankton carbon calculated from ocean color and lidar exhibits basin-scale differences that can reach ±50%.

There are currently three ways to measure b bp globally: (1) autonomous profiling floats (Bittig et al., 2019), (2) passive (or "ocean color") satellite remote sensing, and (3) active satellite remote sensing (light detection and ranging, lidar). Recent progress in using b bp derived from satellite lidar measurements (hereafter "lidar b bp ") to study ocean biology (Behrenfeld et al., 2013(Behrenfeld et al., , 2017Lu et al., 2016Lu et al., , 2020 has established lidar as a prominent tool, such that we are now entering a "satellite lidar era in oceanography" (Hostetler et al., 2018). There is little doubt that lidar b bp observations will continue to advance our understanding of ocean processes because lidar b bp can observe polar ecosystems in the absence of sunlight and at low sun angles, potentially provide constraints on inversion algorithms for passive remote sensing approaches, and contribute another independent measurement of ecosystem stocks.
For decades, it was not possible to assess passive satellite performance of b bp retrievals on global scales because there were few in situ observations. For example, the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) onboard Aqua was launched in 2002, but its most spatially extensive b bp performance assessment was not realized until 2019 (K. M. Bisson et al., 2019), when a global network of Argo floats equipped with backscattering sensors was used. Here, we conduct a similar analysis of lidar b bp .
Identifying in situ matchup observations with MODIS is straightforward relative to lidar because passive ocean color satellite instruments produce wide swaths of data, often stretching 2,000 km in the cross-track direction. In contrast, defining matchups for lidar and in situ observations is challenging because a single lidar pulse, like an in situ measurement, gives a snapshot in time for pinpricks in space. On regional scales, MODIS and lidar have been compared before in polar regions and in the North Atlantic (Behrenfeld et al., 2017;Lacour et al., 2020). Satellite lidar has not been validated globally. In this study, we ask, "Are ocean color and lidar b bp retrievals consistent on global scales?" Satellite b bp (λ) needs to be contextualized with known biases and assessed errors because the fidelity of past and future modeling efforts relies on the accuracy of b bp as an input product. Here, we introduce a method to globally validate lidar backscattering from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument aboard the NASA Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite, and we compare CALIOP b bp , MODIS b bp , and Argo b bp with the goal of quantifying satellite b bp performance and bias.

Argo b bp
Vertical profiles of b bp (700 nm, m −1 ) were downloaded from the Argo Data Assembly Center (ftp://ftp.ifremer.fr/ifremer/argo/dac/ on May 20, 2020) and processed as in K. M. Bisson et al. (2019). Float b bp profiles were despiked with a three-point moving median and outliers in log-space were removed (given by those b bp values outside the bounds of 1.5 times the interquartile range). After outliers were removed, there were 37,337 data points at independent locations. To make the float profiles comparable to remote sensing products (where CALIOP data represent a fixed 22.5 m layer, and MODIS data are exponentially weighted toward the surface), the mean b bp value is reported within the calculated mixed layer depth (MLD, given by the depth where density is greater than 0.03 kg m −3 relative to the density at 10 m). The median MLD is 18 m for the global Argo data set, with an interquartile range of 3.9 m. Choosing the first light attenuation layer rather than the MLD did not significantly change the values of Argo b bp .

Retrieving Ocean Color b bp
The retrieval of ocean color b bp (λ) is an ill-posed inverse problem that requires spectral remote sensing reflectances (R rs (λ); sr −1 ) as input and is constrained with a set of assumptions about the absorbing and backscattering constituents in the ocean. Our analysis is focused on the MODIS instrument onboard Aqua and the generalized inherent optical properties algorithm in its default configuration (GIOP-DC, Werdell et al., 2013) because MODIS outperformed the other contemporary global ocean color satellites such as Visible Infrared Imaging Radiometer Suite (VIIRS) and Ocean and Land Color Instrument (OLCI) in K. M. Bisson et al. (2019). Likewise, the GIOP-DC outperformed the other inversion algorithms such as Quasi Analytical Algorithm (Lee et al., 2002) and Garver-Siegel-Maritorena algorithm (Maritorena et al., 2002) when confronted with Argo float b bp in K. M. Bisson et al. (2019).
We obtained MODIS Level-3 9-km remote sensing reflectance data (R rs, λ = 412, 443, 488, 531, 547, and 667 nm) to generate global b bp maps using the GIOP-DC algorithm, as well as MODIS Level-2 1-km R rs (same wavelengths, all from the NASA Ocean Biology Processing Group, https://oceancolor.gsfc.nasa.gov) to generate coincident matchups with Argo b bp according to the Bailey and Werdell (2006) quality control criteria. We identified satellite matchups that occur within a ±3-h window in a 5 × 5 satellite pixel box, and also within a ±24-h window in a 9 × 9 pixel box centered on the float observation (where the larger box accounts for assumed advection). All R rs (λ) data were corrected to remove Raman scattering through the empirical algorithm of Lee et al. (2013).
GIOP is a flexible inversion algorithm that allows different formulations within the framework to be modified (full details in Supplementary Material Text S1). We ran the GIOP-DC on MODIS R rs observations and report our derived b bp at 532 nm so that MODIS and CALIOP b bp are compared at the same wavelength. Finally, because MODIS b bp is a function of eigenvector choices for b bp (λ), phytoplankton absorption (a ph (λ); m −1 ), and nonalgal particle plus colored dissolved organic matter absorption (a cdm (λ); m −1 ), we performed a sensitivity analysis to quantify MODIS b bp performance depending on which specific assumptions are used (see Supplementary Material Text S2, Table S1).

Satellite Lidar Retrievals of b bp
The CALIPSO satellite was launched in 2006 with the primary goal of observing the vertical distribution of clouds and aerosols. Like MODIS, CALIPSO flies in the A-train constellation and has a 16-day repeat cycle (Winker et al., 2009). CALIPSO's main instrument is CALIOP, which is a nadir-pointing lidar with two measurement wavelengths, 532 nm and 1,064 nm, and has a footprint diameter at the ocean surface of ~100 m. CALIOP measures the copolarized and cross-polarized channel component of column integrated backscatter. Although CALIOP was not intended for ocean research, its polarization properties have been used to characterize b bp at 532 nm for the first vertical 22.5-m bin in the ocean (Behrenfeld et al., 2013). Since 2013, there have been refinements to the lidar b bp algorithm. In this analysis, we use the daytime lidar product published in Behrenfeld et al. (2019), which is freely available online (data access details are in the acknowledgments, and data processing details are in Supplementary Material Text S3). We made one key modification to the In our study we choose a beta(π)/b bp value of 0.32. Because the Behrenfeld et al. (2019) CALIOP data were processed using a value of 0.16, we multiply the retrieved b bp product by 0.5. Using this factor, the global b bp frequency distributions between CALIOP and MODIS are similar ( Figure S1). Given this, we focus our efforts on point by point comparisons of spatiotemporal matchups common to CALIOP, MODIS, and Argo observations. Argo b bp is used from 2015 to present, and the CALIOP b bp product used in this study spans 2006-2017, so we restrict our analysis to 2015-2017.

Identifying Matchups Across CALIOP, MODIS, and Argo b bp
Observations from CALIOP and Argo are single points separated by distance and time, so we could not use a method that relies on their intersection for comparison. Instead, we adopted a decorrelation approach to quantify near-coincident space (Figures 1(a) and 1(b), black and purple lines) and time windows (Figures 1(c) and 1(d), black and purple lines) that yield sufficient matchups (cyan lines) for analysis. Rather than group all data together, we chose to subset regions by annually averaged sea surface temperature (SST), where an SST of 15°C was used to distinguish warmer waters that are permanently stratified within the euphotic layer from cooler, high latitude waters with deeper active mixing (after Behrenfeld et al., 2006). Regions with annual SST < 15°C (Figures 1(a) and 1(c)) and >15°C (Figures 1(b) and 1(d)) represent different physical environments because the first baroclinic Rossby radius of deformation (defining the length scale of baroclinic variability) is dependent on the Coriolis parameter (and therefore on latitude, Chelton et al., 1998). Higher latitudes are expected to exhibit shorter decorrelation length scales of physical variability, which are expected to influence the decorrelation in b bp .
We calculated distances (in km) between Argo and CALIOP using the haversine formula. Pearson's correlation (r) is used to quantify similarity between CALIOP and Argo b bp on log-10 transformed data. We defined coincidence with MODIS following Bailey and Werdell (2006). Backscattering spectral slopes (γ) calculated as part of the GIOP-DC inversion were applied to the Argo b bp at 700 nm to derive Argo b bp at 532 nm so that all b bp are comparable at the same wavelength (see Equation 4 in Supplementary Material Text S1). We calculated the median percent error (MPE, the median of 100% × |satellite b bp /Argo b bp -1|) and relative bias (the median of 100% × [satellite b bp − Argo b bp ]/Argo b bp ) to compare MODIS, CALIOP, and Argo b bp (all at 532 nm; m −1 , and data are not log-transformed prior to these calculations because the data are normalized by Argo b bp ).
We used the shapes of decorrelation for r and MPE to find cutoff distance values where the slope of MPE increases and the slope of the correlation decreases (Figure 1). The intent of this approach was to maximize the number of paired observations while maintaining a high correlation and low MPE. Based on this analysis, we chose 15 km radius for matchups in regions with annual SST < 15°C and 50 km for regions with annual SST > 15°C (Figure 1). In all cases, the correlation is similar across all hours (up to 24). With these matchup criteria, we take a subset of the Argo observations common to both CALIOP and MODIS b bp (within a ±3-h window, n = 93 as well as ±24-h window, n = 261) so that all three sensor types can be compared. One alternative approach to the paired matchup method as outlined here is to look at general correspondence between distributions of CALIOP and Argo b bp in particular regions, as is done in Lacour et al. (2020) in the North Atlantic. Differences between Lacour et al. (2020) and the current study are discussed further in Supplementary Material Text S4.

Results
The spatial distribution of matchup Argo b bp observations common to both CALIOP and MODIS exhibits good global coverage with representation in the Southern Ocean, Arctic Ocean, South Pacific Gyre, Atlantic basin, and Indian Ocean (Figure 2(a)). Global annually averaged maps of MODIS and CALIOP b bp reveal similar patterns (Figures 2(b) and 2(c)), with elevated b bp in coastal and/or upwelling regions and lower b bp in the oligotrophic gyres.
Evaluation of equivalent matchup data between Argo observations of b bp and retrievals from CALIOP and the optimum parameterization of MODIS GIOP-DC reveals a superior performance of CALIOP (black bars vs. purple bars in Figure 3), having lower MPE (Figures 3(a) and 3(c)) and relative bias (Figures 3(b) and 3(d)). This improved performance is observed both for the 3-h matchup data (where CALIOP has 18% MPE vs. 24% MPE for MODIS [Figures 3(a), S1(a), and S1(c)] and CALIOP has a lower relative bias [−11%] compared to MODIS [−21%] [ Figure 3(b))] and the 24-h matchup data (Figures 3(c), 3(d), S2(b), and S2(d)). CALIOP also exhibits superior performance in b bp retrievals compared to four alternative ocean color inversion (GIOP) parameterizations (MPE and relative bias bars to the right of the vertical red line in Figure 3).

We find clear inconsistencies between CALIOP and MODIS b bp (Figures S2(c) and S2(g)).
For b bp values around 0.001 m −1 , MODIS exhibits a higher dynamic range of b bp compared to CALIOP, spanning nearly an order of magnitude in the ±24-h case ( Figure S2(g)).
In a qualitative sense, MPE and bias are better for both sensors in the ±3-h window compared to the ±24-h window, and both MODIS and CALIOP underestimate Argo b bp in general (Figures 3(b), 3(d), S2, and S3).
As an illustration of the ecological and biogeochemical significance of CALIOP and MODIS b bp differences, we converted these data into estimates of PhytoC concentrations using the linear relationship reported by Graff et al. (2015). While annual global average PhytoC estimates from MODIS (PhytoC M ) and CALIOP (PhytoC C ) are similar (17 mg C m −3 and 18 mg C m −3 , respectively), notable regional differences are observed ( Figure  4). For example, PhytoC C is ~20% higher than PhytoC M in the South Pacific gyre and temperate regions of the North Pacific and South Atlantic. In contrast, PhytoC M exceeds PhytoC C by ~20% in the Equatorial Pacific and the central gyres of the South Atlantic and Indian Oceans. The largest differences between retrievals are found in the North Indian Ocean, the equatorial Atlantic west of Africa, and the Artic/Subarctic, where PhytoC C may exceed PhytoC M by up to 50%, and in the Southern Ocean where PhytoC M exceeds PhytoC C by 50%.

Discussion
The improved performance of lidar b bp retrievals relative to ocean color reported here is a somewhat unexpected finding because spatial coverage of lidar data is so restricted compared to ocean color data. In other words, the average spatiotemporal coincidence between Argo b bp data and lidar is far broader than that for wide-swath, 2-day repeat cycle ocean color measurements, suggesting (a priori) that ocean color b bp should yield better performance when compared to float data, at least in spatiotemporal heterogenous waters. Instead, the relatively low MPE (18%) and relative bias (−11%) for the CALIOP data is a clear improvement over all contemporary ocean color satellites, even when considering their highest performing b bp algorithm (MODIS-24% MPE, bias = −21%, this study, VIIRS-31% MPE, bias = −11% [K. M. Bisson et al., 2019], and OLCI-45% MPE, bias = 2% [K. M. Bisson et al., 2019], with biases recalculated according to our definition here, Figure S1). Although the study of K. M. Bisson et al. (2019) featured more matchups between Argo and ocean color b bp , we note that lidar MPEs are below 25% even for ±24-h matchups at distances >50 km in the SST > 15°C case, which represents ~75% of the ocean by area. We also note that MODIS performance degrades with distance and time (as expected), with an MPE of 27% in the 9 × 9 pixel box, ±24-h matchup case. While the performance of MODIS b bp is indeed affected by choices within the GIOP inversion, no particular configuration can produce the performance metrics of CALIOP b bp . If more exact spatial matchups were possible between lidar and Argo b bp data, we would expect the enhanced performance of lidar compared to ocean color to be even more pronounced.
Although we have found good agreement between CALIOP and Argo b bp , CALIOP b bp is imperfect, particularly at low Argo b bp values ( Figure S4). Future lidar products may especially benefit by optimizing b bp to float values (as is done currently with the MODIS SST algorithm, Kilpatrick et al., 2015), especially because there are sufficient (~750) matchups between CALIOP and Argo observations. CALIOP b bp is also sensitive to the scattering phase function used (which might vary regionally/seasonally) and data are only available along the orbit track as opposed to the large ocean color swaths.
There are necessarily limitations of ocean color b bp . Lidar is a more direct measurement of b bp compared to MODIS, as the latter retrieval uses the remote sensing reflectances, R rs (λ), with assumptions about the absorption and backscattering spectral shapes of the ocean components and specific relationship between R rs (λ) and inherent optical properties. R rs (λ) is retrieved following atmospheric correction, which removes radiometric effects from ocean surface glint and white-caps, as well as molecular and aerosol absorption and scattering. The chemical composition and size distribution of aerosols are assumed (Gordon, 1997;Gordon & Wang, 1994) and inferring the aerosol signal from satellite observations can be challenging since the atmospheric signal contribution is typically 90% at 440 nm at the top of the atmosphere while the residual signal is from the ocean. Even worse, the contribution of the ocean signal quickly decreases at longer wavelengths (i.e., >500 nm), making it more challenging to accurately estimate R rs (λ). While useful for sensor-to-sensor comparisons, the bidirectional reflectance distribution function correction, as part of the atmospheric correction, can impart additional uncertainty in R rs (λ) retrievals, as it cannot ubiquitously represent all conditions at all times (Mobley et al., 2016). Small uncertainties in the aerosol correction lead to large uncertainties in R rs at green and red bands due to two inherent limitations: (1) the ocean signal is small relative to the aerosol signal and (2) the dynamic range of R rs in the green and red wavelengths is small compared to the more dynamic aerosol signal. A future assessment is required to quantify the impacts of atmospheric correction on b bp retrievals. Despite the issues outlined above, we find a good overall correspondence between ocean color and Argo b bp. CALIOP, Argo, and MODIS observe b bp in different areas of the volume scattering function. CALIOP is nadir viewing (scattering angle of 180°), typical backscattering sensors used on Argo float have a nominal scattering angle of 142° (but some have 124° and 149°, Poteau et al., 2017), while MODIS has viewing angles relative to nadir spanning ±49.5° (corresponding to scattering angles between 131° and 180°; https://aqua.nasa.gov/modis). As R rs is known to be influenced by viewing angle and particle phase function (with variations up to 65% in some cases, Xiong et al., 2017), the viewing angle differences between sensors are a potential source of error for the retrieved b bp products. Another source of discrepancy between sensors is the water column depth used to generate b bp observations. MODIS, CALIOP, and Argo consider slightly different depths of the water column, which may be potentially important for instances when the sensor depth exceeds the mixing depth. A further source of error arises from the different sensor wavelengths used in this study. Given that ocean color b bp performance is affected by assumed backscattering and absorption shapes, it would be a meaningful improvement for future floats to be equipped with a backscattering sensor including green wavelengths. Having Argo b bp observations in the green would eliminate the influence of the backscattering spectral power-law-fit slope assumption and also the influence of absorption assumptions because the green bands are minimally influenced by phytoplankton and water absorption.
In this study, we validated global lidar b bp and compared it to the best case ocean color sensor and algorithm pairing. Regional differences in derived PhytoC between CALIOP and MODIS quantify the consequences of b bp product choice. PhytoC is essential for calculating phytoplankton physiology through the chlorophyll:PhytoC ratio and it is a central term in state-of-the-art NPP algorithms (i.e., Silsbe et al., 2016;Westberry et al., 2008) and carbon export models (where differences in data products have wide effects on model outcomes, K. M. Bisson et al., 2018). Although there are clear spatial differences between CALIOP and MODIS b bp , we choose not to focus on regional differences within our analysis because there are too few observations (93 globally, with only 5 observations poleward of 50 degrees) to make rigorous statements about performance on regional scales.
Because CALIOP has limited spatial coverage compared to MODIS, an optimal approach may come from generating products that combine data from the two sensors. Continued efforts are needed to improve CALIOP lidar retrievals in low b bp areas. Nevertheless, the CALIOP record provides a less uncertain and an independent global data set of b bp that presents an opportunity for evaluating and improving satellite ocean color retrievals of this fundamental optical property related to plankton ecosystem structure and biogeochemistry.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.    Comparison of CALIOP performance metrics and those for variants of MODIS inversions. V1 changes the b bp slope used, V2 changes the a cdm shape, V3 changes the assumed a ph shape, and V4 is the GIOP result using R rs data that were not corrected for Raman scattering.