Approach for identifying optically shallow pixels when processing ocean-color imagery.

Ocean-color remote sensing is routinely used to derive marine geophysical parameters from sensor-observed spectral water-leaving radiances. However, in clear geometrically shallow regions, traditional ocean-color algorithms can be confounded by light reflected from the seafloor. Such regions are typically referred to as "optically shallow". When performing spatiotemporal analyses of ocean color datasets, optically shallow features such as submerged sand banks and coral reefs can lead to unexpected regional biases. Most contemporary approaches mask or flag suspected optically shallow pixels based on ancillary bathymetric data. However, the extent to which seafloor reflectance contaminates the water-leaving radiance is dependent on bathymetry, water clarity and seafloor albedo. In this paper, an approach for flagging optically shallow pixels has been developed that considers all three of these variables. In the method, the optical depth of the water column at 547 nm, ζ(547), is predicted from bathymetric data and estimated water-column optical properties. If ζ(547) is less then the pre-defined threshold, a pixel is flagged as potentially optically shallow. Radiative transfer modeling was used to identify a conservative threshold value of ζ(547) = 20 for a bright sand seafloor. In addition, pixels in waters shallower than 5 m are also flagged. We also examined how varying bathymetric datasets may affect the optically shallow flag using MODIS data. It is anticipated that the optically shallow flag will benefit end-users when quality controlling derived ocean color products. Further, the flag may prove useful as a mechanism for switching between optically deep and shallow algorithms during ocean color processing.


Introduction
Ocean-color inversion algorithms are routinely used to derive inherent optical properties, IOPs, from the sensor-observed spectral remote-sensing reflectances, R rs (λ) (sr −1 ). Indeed, a range of ocean-color inversion algorithms have been developed that are robust and able to retrieve IOPs with high precision and low bias [1]. Many of these approaches [2][3][4][5][6][7][8] have been implemented within the NASA Ocean Biology Processing Group's L2GEN code, which is distributed as part of the NASA SeaWiFS Data Analysis System (SeaDAS) software suite (https://seadas.gsfc.nasa.gov). SeaDAS offers a user-friendly environment for processing, visualizing, and analyzing ocean-color data.
Most ocean-color algorithms have been developed for optically deep, oceanic waters where the IOPs, namely the total absorption coefficient, a(λ) (m −1 ), and total backscattering coefficient, b b (λ) (m −1 ), are deemed to co-vary with phytoplankton abundance. However, efforts continue on developing ocean-color inversion algorithms for optically complex and optically shallow waters [1,9]. Optically complex waters are characterized as regions where the IOPs are influenced by phytoplankton, as well as non-algal particulate matter (NAP), and colored dissolved organic matter (CDOM), whose concentrations do not necessarily co-vary with phytoplankton abundance [9,10]. Optically shallow waters are those in which light reflected off the seafloor contributes significantly to the water-leaving signal [11]. This is known to effect geophysical variables derived by ocean-color algorithms, often leading to biased values [7,[12][13][14]. For example, Fig. 1 shows how the Pedro Bank, a highly reflective shallow region southwest of Jamaica, and other shallow features are associated with anomalously high retrievals of chlorophyll-a pigment concentrations (Chla; mg m −3 ). We note that the Chla image in Fig. 1 was derived using the standard NASA algorithm as described in the online Algorithm Theoretical Basis Document (ATBD; https:// oceancolor.gsfc.nasa.gov/atbd/chlor_a/). When using L2GEN to process from top-of-atmosphere radiances (level-1) to geophysical products (level-2), binary quality control flags (l2_flags) are assigned to each pixel [15]. The l2_flags are most often used to identify pixels of suspected poor quality, thereby allowing the end-user to exclude them from further analysis. Pixel exclusion is particularly important when performing algorithm validation [16] and/or spatiotemporal binning and mapping. Presently, a coastal zone l2_flag (COASTZ) exists to flag geometrically shallow continental shelf pixels, with the general assumption that these pixels are potentially optically shallow. COASTZ uses a 1 arc-minute resolution global bathymetric dataset and flags pixels with geometric watercolumn depths less than 30 m.
A number of radiative transfer-based studies have discussed how light reflected off the seafloor affects R rs (λ) [7,[17][18][19][20]. The extent to which R rs (λ) is influenced by light reflected off the seafloor is indeed dependent on the geometric depth of the water column, but also the benthic albedo (ρ; unitless) and the IOPs of the overlying water column. For example, in very clear waters up to 20 m in depth, a bright sandy seafloor may contribute to R rs (λ).
However, if the water column is highly attenuating, the same sandy seafloor at 20 m may have little-to-no influence on R rs (λ) [7]. While bathymetry alone has most often been used to flag suspected optically shallow waters, IOPs are generally considered highly variable in time and space for most coastal locations. Following, it is likely that in a geometrically shallow region the water may be optically shallow at some times and optically deep at others as the water column opacity evolves (through its variable IOPs). The COASTZ flag that is based solely on bathymetry is therefore insufficient for identifying optically shallow pixels; information regarding IOPs is also necessary.
In the spectral range 470 -570 nm, often referred to as a 'transparency window', light reflected off the seafloor is expected to have the most influence on the water-leaving radiance signal [21]. We therefore expect that knowledge of water column clarity near 550 nm (547 nm for NASA's Moderate Resolution Imaging Spectroradiometer aboard the Aqua spacecraft; MODISA) will inform of whether it is likely that R rs (λ) is radiometrically MCKINNA  influenced by seafloor reflectance. An IOP-based estimate of the optical depth at 547 nm, ζ(547) (unitless) would provide an appropriate metric to quantify water-column opacity, however, a direct estimate of the IOPs at 547 nm is not viable due to potential radiometric contamination of R rs (547) by seafloor reflectance. Instead, we propose to follow the approach of Barnes, et al. [22] and estimate IOPs at 667 nm, a spectral region less susceptible to shallow water effects, followed by extrapolation and estimation of IOPs at 547 nm using a bio-optical model. With a known geometric water-column depth, z, it is then possible to approximate ζ(547). If the water is suitably opaque, we can assert that the benthic contribution is negligible, however, a suitable threshold value of ζ(547) needs to be established. In this study, the threshold depth for ζ(547) will be determined by analyzing the results of radiative transfer modeling.
Within this paper, we use radiative transfer modeling to examine the effect of optically shallow water upon R rs (λ). We then use the outcomes of this study to develop a level-2 quality control flag for identifying potential optically shallow pixels. The flag is demonstrated using sample imagery of the northern Great Barrier Reef (GBR), Australia. We also discuss the effect of varying the input bathymetric dataset used for calculating the flag.

Optical depth
At a given geometric depth, z, a quantitative measure of the water column's vertical opacity is the optical depth, ζ(λ), which is a non-dimensional spectral parameter with values ranging from 0 to ∞ [23]. The more transparent the water column, the smaller ζ(λ) will be. Mathematically, ζ(λ) is related to z and the beam attenuation coefficient, c(λ), via the following relationship [23] ζ(λ) = ∫ 0 z c λ, z′ dz′ . (1) If we consider a shallow water column as optically homogenous, then ζ(λ) over a geometric depth range, Δ, can be expressed as Should the value of ζ(λ) be suitably small, it is likely that light reflected of the seafloor will contaminate R rs (λ). However, the extent to which R rs (λ) is affected is also determined by the reflectivity of the seafloor. For this study, we considered three benthic substrate types: seagrass, coral, and sand.

Radiative transfer modeling
Within this study we used Hydrolight-Ecolight 5.1 (HE5) radiative transfer (RT) software [24] to simulate R rs (λ) at fourteen geometric depths ranging from 3 to 50 m, using 250 IOP combinations. The synthesized IOP dataset established by the International Ocean Color Coordinating Group (IOCCG) [25] was used as optical inputs for HE5 RT modeling. The first 250 IOP combinations from the IOCCG dataset were used as these were deemed representative of oligotrophic and mesotrophic optical conditions encountered in shallow shelf waters (0.03 < Chla < 1.0 mg m −3 ; 0.2 < c(547) < 1.3 m −1 ). RT modeling followed the method described in IOCCG Report 5 [25], however, the water column depth was varied and three different benthic albedos were used. Details of the HE5 modeling are given in Table 1. Note, inelastic scattering processes were not included in these HE5 simulations.
During simulations, three reflective benthic albedos of varying brightness were used (i) sand, (ii) coral, and (iii) seagrass. Simulations were also conducted for a "black" nonreflective benthic substrate. Thus, for each depth and IOP combination a reflective albedo remote sensing reflectance, R rs,Ref (λ), and a black-seafloor remote-sensing reflectance, R rs,Blk (λ), were simulated. By modeling pairs of R rs,Ref (λ) and R rs,Blk (λ) spectra, it was then possible to quantify the relative contribution, RC(λ), of light reflected off the seafloor to the total remote sensing reflectance signal as By computing RC(λ) values for each depth and IOP combination, it was thus possible to see how RC(λ) for a given seafloor type varied with ζ(λ). In addition, the remote sensing reflectance for infinitely deep water, R rs,Deep (λ), was computed for all IOP combinations. Using these data it was possible to calculate the relative difference, RD(λ), between R rs,Ref (λ) and R rs,Deep (λ) for each IOP combination, geometric depth and benthic albedo type, where From considering RC(λ) and RD(λ) values, a suitable ζ(547) threshold could be determined for use in the optically shallow flag.

Satellite-estimates of the optical depth
An approach was developed to estimate ζ(547) using ancillary bathymetry data and IOPs (a(547) and b b (547)) derived using the Quasi-Analytical Algorithm (QAA) [4]. The QAA derives IOPs by algebraically solving the semianalytical model developed by Gordon et al. [26] that simulates R rs (λ) as a function of IOPs. In this study, we have used a modified form of the QAA deemed suitable for optically shallow waters described by Barnes, et al. [22]. For brevity, we have not included a full description of the QAA here, please see Lee, et al. [4] and Barnes, et al. [22] for comprehensive detail.
Once a(547) and b bp (547) are known, the beam attenuation coefficient is estimated as where, the value of the backscatter ratio of pure water, η w , was set to 0.5, and the particulate backscatter ratio, η p , was set to 0.015; halfway between global average oceanic value of 0.01 [27] and well known Petzold average particle value of 0.0183 [28]. Finally, the estimated optical depth at 547 nm, ζ E (547), is calculated as follows where the water column depth, Δz, is taken from an ancillary bathymetric dataset. The final step in developing the flag was to determine an ζ E (547) threshold value below which the water column is considered optically shallow. This threshold value was determined to be 20 by analyzing HE5 radiative transfer modeling outputs, as described in section 3.1.

Test scenes and data processing
All ocean color data used in this study were sourced from the NASA Ocean Biology Processing Group's most recent version of the MODISA dataset (v2018.0) [29]. Data were processed using the default configuration of L2GEN. For testing and evaluation purposes, the optically shallow flag was implemented as a sub-module of the Shallow Water Inversion Model (SWIM) [7]. We initially evaluated the impact of varying the input bathymetric dataset used by the flag. Specifically, we used three freely-available bathymetric datasets: (i) the global ETOPO1 1-arc-minute combined ocean bathymetry and land topography dataset (ETOPO1) [30], (ii) the global General Bathymetric Chart for the Oceans 30 arc-second dataset (GEBCO) [31], and (iii) the region-specific GBR 100 m dataset (3DGBR) [32].
It was also determined how a flagged area varies temporally based on changes in water clarity. To do so, we applied the optically shallow flag to four MODISA test images of the northern GBR, Australia near Lizard Island (14.67°S, 145.46°E). During the austral winter/ spring, when river outflows are least, water clarity in parts of the GBR is governed by periodic wind-, wave-, and tide-driven sediment resuspension events [33]. Accordingly, we processed four test scenes from July-August 2015 during which time wind speed, as measured by an Australian Institute of Marine Science (AIMS) weather station on Lizard Island, varied in magnitude and direction.

Radiative transfer modeling results
The two-dimensional sub-plots in Figs. 2 and 3 show how RC(λ) and RD(λ), respectively, vary with c(λ) and geometric depth. Figure 2 clearly demonstrates the dependence of benthic albedo contribution to R rs (λ) on both water clarity and water-column depth, as well as on seafloor type. RC(λ) for sand is greater than that for seagrass, for example, particularly when ζ(λ) < 5. Figure 3 reveals the relative difference between optically shallow and optically deep values of R rs (λ), namely how RD(λ) varies with water clarity, water-column depth, and benthic albedo. Both figures indicate that sand increases the magnitude of R rs (λ) (> 25 %) at all wavelengths shorter than 667 nm when ζ(λ) < 5.
For the darker benthic albedos, however, the relationship is somewhat different, particularly for bands centered on 412 and 443 nm. For both coral and seagrass, for example, ζ(412) and ζ(443) decreases are accompanied by decreases in RD(λ). We hypothesize this effect is due to pigments (e.g. chlorophyll-a) in the benthic substrates strongly absorbing photons at shorter wavelengths thereby resulting in the R rs,Ref (λ) signal being smaller than the R rs,Deep (λ) signal.
From these plots of RC(λ) and RD(λ), we see that the three MODISA bands in the 'transparency window' centered on 488, 531 and 547 nm were most susceptible to radiometric contribution from the seafloor albedo. In addition, bands centered on the shorter wavelengths 412 and 443 nm were also significantly affected. However, the band centered on 667 nm was least affected, which justifies the use of R rs (667) as a reference wavelength for the QAA as per Barnes, et al. [22]. However, for geometric depths of 5 m or less, R rs (667) is contaminated by benthic albedo. Accordingly, we recommend that any pixel with a water-column depth of 5 m or less be automatically flagged as optically shallow. We note that previous studies [7,22] have suggested caution when processing pixels with associated geometric depths less than 5 m in order to avoid unwanted optically shallow artifacts manifesting in derived ocean color data products.  Table 3.
Generally speaking, a bright sandy bottom represents the 'worst case' scenario in our simulated ocean color data per all figures showing RC(547) and RD(547) and the results in Table 3. Collectively, the results of our RT modeling exercise for sand albedo suggest that benthic albedo effect upon R rs (547) becomes almost completely diminished once ζ(547) > 19. Thus, assigning a threshold value of ζ(547)=20 will conservatively flag any pixel with potential to be optically shallow. We also note that RC(547) and RD(547) diminish almost completely when the geometric depth exceeds 40 m (results not shown).
The RT dataset also facilitated an evaluation of the QAA-based method for estimating ζ(547). A scatter plot of estimated (QAA) versus actual (RT) ζ(547) values is shown Fig. 6.
The results indicate that, to a first-order, the estimated values of ζ(547) were reasonably good with a coefficient of determination (R 2 ) of 0.94 and root mean squared error (RMSE) of 1.45.
We note that in Fig. 6 there is a distinct out-of-family cluster of points where ζ(547) < 5. This cluster of data points is associated with geometric depths less than 5 m. From Figs. 2 and 3, we see that for depths less than 5 m, the sand albedo significantly affects R rs (667). This suggests that the QAA-based approach for estimating ζ(547), which assumes R rs (667) is not radiometrically contaminated, becomes problematic at geometric depths less than 5 m.
We acknowledge that η p , which depends on particle size and type, is variable and may deviate from the fixed value of 0.015 in coastal waters [34][35][36][37]. Indeed, we found that modulating η p by ±0.005 caused ζ(547) estimates to vary in response by approximately ±30% (results not shown). Specifically, using an η p value of 0.010 (the accepted nearsurface global average [27]) estimates of ζ(547) were up to 30% lower and using an η p value of 0.020 estimates of ζ(547) were up to 30% larger. In the context of the flag threshold value of ζ(547)=20, this equates to a potential low-end bias of ζ(547)=14 or a potential high-end bias of ζ(547)=26. A likely impact is that some turbid water pixels (i.e. larger η p ) may be incorrectly flagged. However, it is important to remember the purpose of the flag is to provide a conservative firstorder method for identifying potential optically shallow pixels and the method for selecting η p , either as a fixed or modeled value, can be refined in future as necessary. Alternatively, OPSHAL could be parameterized such that an end-user could input values for η p and also the threshold value for ζ(547) at runtime.
Considering the results of the RT processing, we defined our proposed optically shallow water flag (which we hereafter designate as OPSHAL) by the following binary decisions:

Effect of varying bathymetric datasets
The optically shallow flag is dependent upon the fidelity of the input bathymetry data. To assess how different horizontal and vertical resolutions may affect OPSHAL, we processed a test scene of the northern GBR three times using differing bathymetry datasets with varied spatial resolutions: ETOPO1 (1 arc-minute), GEBCO (30 arc-second), and 3DGBR (100 m). Figure 7 shows the differences in bathymetric detail captured by each dataset and differences in the areas flagged as optically shallow and Table 2 describes the number of shelf pixels flagged as optically shallow.
From Fig. 7, it is clear that ETOPO1 poorly resolves fine-scale submerged features with nearly the entire shelf having geometric depth of approximately 10 m. Further, ETOPO1 results in 90% of all shelf pixels being flagged as optically shallow. By contrast, both the GEBCO and 3DGBR datasets appear to resolve bathymetric features well, resulting in 64 and 58% of all shelf pixels being flagged, respectively. We expect that the region-specific 3DGBR dataset resolved regional bathymetric features best, thus providing the preferred input to the flagging method. We concede, however, that region-specific datasets, while useful for local studies, are not adequate for global data processing. Encouragingly, the spatial area flagged when using the GEBCO dataset in Fig. 7 was very similar to the area flagged when 3DGBR dataset is used. Accordingly, we propose using the GEBCO bathymetry as the ancillary bathymetry dataset for global application of the optically shallow water flag.
We also acknowledge that tidal variation would also affect the geometric depth and hence the estimate of ζ(547). This is particularly pertinent to intertidal waters and lends further merit to automatically flagging all pixels in waters shallower than 5 m. At present in an effort to reduce complexity and maintain computation efficiency we have not included tide as a component of OPSHAL.

Variability in flagged area
We also processed four MODISA images of the northern GBR captured on 31 July, 2 August, 4 August, and 11 August 2015 using 3DGBR bathymetry. Figure 8 shows the area flagged as OPSHAL (orange pixels) on each date. This imagery encompasses a section of the GBR bounded to the west by the northeastern Australian coastline and to the east by the Coral Sea. A broad region of shallow shelf waters extends northward from Cape Flattery (14.946°S, 145.347°E) to Cape Melville (14.178°S, 144.522°E).
The imagery shows that a ribbon of non-flagged water approximately 5 km wide lay parallel to the coastline and the mid-shelf between Cape Bedford and Cape Melville on 31 July (see area inside black ellipses in Fig. 8). During this time the seasonal southeasterly trade winds were intense, with a 24-hour averaged speed of 22 knots, likely to promote sediment resuspension and increase turbidity in the shallow shelf waters. As the wind speed decreased to 18 knots and 15 knots on 2 August and 4 August, respectively, the ribbon of non-flagged water became narrower, indicative of reduced turbidity. Finally, on 11 August the wind speed and direction dropped to 9 knots from an east-southeast direction. Under these conditions, it is likely that wind-induced sediment resuspension and mid-shelf turbidity was reduced and water clarity increased. Indeed, on 11 August when the water seemingly cleared, most of the shelf was flagged as optically shallow.
For the same region shown in Fig. 8, we computed the annual frequency of MODISA pixels flagged as OPSHAL for the year 2015. Figure 9 compares the annual frequency map of OPSHAL-flagged pixels with pixels flagged as COASTZ. The COASTZ flag, which is temporally invariant and based on coarse bathymetry, flags almost all shelf pixels. Conversely, the frequency of pixels flagged as OPSHAL is shown to be temporally variable; aside from geometrically shallow (<5 m) near-shore pixels that are flagged 100% of the time.
Collectively, these results demonstrate how the distribution of optically shallow pixels is dynamic and can vary spatiotemporally. OPSHAL can thus be considered an improvement upon the existing COASTZ quality control flag for shelf waters. Interestingly, these examples also demonstrate how the OPSHAL flag may direct where and when optically deep bio-optical algorithms can be trusted in geometrically shallow shelf waters.

Conclusions
We presented a method for identifying optically shallow pixels in ocean color imagery. Unlike its predecessor, OPSHAL does not rely on bathymetry alone, but rather, calculates marine IOPs on-the-fly to generate time-and space-varying estimates of water opacity. The method is computationally efficient and will be highly useful for multi-level ocean color data processing software such as NASA's L2GEN. The OPSHAL flag will provide an improved level of quality control for end-users, particularly when creating spatiotemporal maps of ocean color data products. We note that the flag may also facilitate the dynamic switching and/or blending of bio-optical algorithms developed for optically deep and optically shallow waters.
It is important to recognize that the approach is limited by the fidelity of the input bathymetry dataset as demonstrated in section 3.2. The GEBCO 30 arc second dataset provides a spatially adequate input dataset for global data processing. However, we concede that GEBCO may not resolve complex marine topography in all areas, thereby limiting the utility of OPSHAL. However, this may be alleviated by end-users utilizing region-specific high quality bathymetric datasets. Indeed, recent work by Barnes,et al. [38] demonstrates how bathymetry can be derived from multi-band ocean color imagery.
Although demonstrated here using MODISA data, the OPSHAL flag can conceptually be applied to past and present ocean color datasets collected by sensors with bands centered on or near 670 nm and 550 nm, including but not limited to: the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), and the Visible Infrared Imaging Radiometer Suite (VIIRS), the Ocean and Land Colour Instrument (OLCI), and the Second Generation Global Imagery (SGLI). We note the approach may also be applicable to future datasets collected by NASA's Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission, which is currently under development and expected to launch in 2022.         Left-hand column shows MODISA true color imagery of the northern Great Barrier Reef, Australia captured on four different dates spanning two weeks. Central column: MODISA imagery with optically shallow pixels colored orange. Clouds over water or pixels with invalid values of R rs (667) (e.g. negative or atmospheric correction failure) values are colored red. Black ellipses indicate non-flagged pixels that lie parallel to the coast. Righthand column: daily-averaged wind speed and direction measured by the AIMS weather station located on Lizard Island.