Fisheries and Oceans Canada maintains a network of scientific buoys in the St. Lawrence estuary and gulf. Among a suite of environmental parameters documented, the in-water upwelling radiance is measured using an ocean color radiometer located underneath the center of the buoy. The shadow effect from the 1.05-m-radius buoy on the measured upwelling radiance is estimated and empirical models to correct for it are proposed. On average, the shading error (i.e., the percent of missing radiance) was and for the 555- and 412-nm channels, respectively. Two analytical models were tested to predict the shading error using measured inherent optical properties, the sun zenith angle, and the fraction of diffuse sky irradiance. Neglecting light scattering led to overestimates of the shading error. In contrast, the bias was removed when the scattering coefficient was accounted for, but the overall error was only barely improved (root-mean-square error ). Empirical relationships based on the uncorrected reflectance ratio measured by the buoy were used to predict both the shading error and the diffuse attenuation of the upwelling radiance, two quantities needed to calculate remote sensing reflectance . Overall, was retrieved with an averaged absolute percent difference ranging from 12% to 20%, which appears adequate for the validation of ocean color data such as the Moderate Resolution Imaging Spectroradiometer (MODIS-Aqua) and Visible Infrared Imaging Radiometer Suite (VIIRS) products in the optically complex waters of the St. Lawrence estuary.
Satellite ocean color missions rely on accurate in situ measurements to ensure data quality through a calibration and validation (CAL/VAL) program (Zibordi et al. 2015). Autonomous platforms, such as the Marine Optical Buoy (MOBY) (Clark et al. 1997) and the Bouée pour l’acquisition de Séries Optiques à Long Terme (BOUSSOLE) (Antoine et al. 2008), have been deployed and maintained to provide high-quality data for CAL/VAL purposes in clear oceanic waters. For optically complex waters, a cost-effective solution has been developed by the ocean color component of the Aerosol Robotic Network (AERONET-OC) (Zibordi et al. 2004, 2009). In this network, the water-leaving radiance is estimated using a single instrument, the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) Photometer Revision for Incident Surface Measurements (SeaPRISM), which performs multiple sky- and sea-radiance measurements at given viewing and azimuth angles and at multiple visible and near-infrared wave bands (for details see Zibordi et al. 2009).
The estuary and Gulf of St. Lawrence (EGSL) is among the most productive coastal zones in Canada and North America, sustaining high planktonic biomass during the summer season (Levasseur et al. 1992; Plourde and Runge 1993; Therriault and Levasseur 1985; Parent et al. 2015). The spatiotemporal variability of phytoplankton in the EGSL has been evaluated using ocean color remote sensing data from the coastal zone color scanner (CZCS) in the mid-1990s and later using SeaWiFS (Fuentes-Yaco et al. 1997b,a; Le Fouest et al. 2006). However, the satellite-based quantitative assessment of phytoplankton biomass and primary production in the EGSL is still lacking. Using a limited set of 39 matchups distributed throughout the EGSL, Yayla (2009) showed that the primary ocean color product— that is, the spectral remote sensing reflectance or normalized water-leaving radiance —was not well estimated using SeaWiFS data. Because of the optical complexity of these waters (Nieke et al. 1997; Le Fouest et al. 2006; Cizmeli 2008), a better assessment of the radiometric products currently distributed by various space agencies is required.
The Department of Fisheries and Oceans Canada (DFO) has implemented monitoring programs of its territorial waters. In the EGSL, a network of four to five moored buoys that carry a suite of hydrographic and meteorologic sensors has been deployed (Fig. 1b). Data are either stored internally or transmitted to shore, and some are posted online in near–real time (NRT). In 2002, the buoys were equipped with multispectral radiometers with seven visible-spectrum channels to measure above-water incident irradiance () and in-water upwelling radiance () in order to provide in situ data for satellite validation purposes. The measurements, however, suffer from the shadow created by the buoy in the underlying water column. In addition, the upwelling radiance is measured at a fixed distance from the sea surface (0.86 m) and a correction scheme is needed to properly estimate the right amount of energy leaving the water, that is, the water-leaving radiance . A validated processing protocol for the radiometric quantities provided by the network is still lacking, which is one of the motivations for the present study. Another motivation for this work was to develop a robust system capable of providing NRT in situ data for the validation of future ocean color satellites, including the recently launched Ocean and Land Color Imager (OLCI) on board the European Space Agency’s Sentinel-3 satellite series.
The objectives of this study were 1) to quantify the shadow effect on measured upwelling radiance under the buoy structure using in-water light profiles obtained at the mooring site over a complete summer season, 2) to test analytical models and to develop an empirical model to predict the shadow error, and 3) to implement an algorithm to process the radiometric data from the buoy in order to derive remote sensing reflectance in an autonomous mode. The two biggest challenges faced were (i) to correct the for the buoy shadow and (ii) to estimate given that is measured at a fixed depth below the sea surface. To achieve these objectives, one of the buoy sites (IML-4, Fig. 1) was visited on a weekly basis during the 2015 summer season during which apparent optical properties (AOPs) and inherent optical properties (IOPs) were measured together with key biogeochemical parameters.
In the following, we first present the theoretical background on instrument self-shading for underwater optical measurements, which is the main issue to derive accurate radiometric products from the buoy radiometric data. Next, the buoy’s characteristics and field measurements are presented in detail. Based on in situ observations, the shadow effect was determined experimentally and modeled using semianalytical solutions. We further present an empirically based scheme to estimate from the buoy’s radiometric measurement and to discuss the main sources of its uncertainty. Finally, we demonstrate the potential of the DFO buoy network to validate the primary satellite products from MODIS-Aqua and VIIRS.
Any in-water measurement of upwelling radiance suffers from instrument self-shadow, which itself reduces the light intensity that could potentially be backscattered toward the sensor. This problem was first described by Gordon and Ding (1992) for an instrument that can be considered as a single two-dimensional disk. Gordon and Ding (1992) defined the measurement shadow error, ε, as
If we consider the instrument as a single two-dimensional disk and ignore scattering processes, then ε can be approximated using the following expression (Gordon and Ding 1992):
where a is the total absorption coefficient of the medium and is the vertical depth affected by the shadow of the instrument below the sensor located at depth z. The value of 2 comes from the following: 1) the light at depth is the radiance at depth z that has been attenuated from that depth to —so, for a sensor located at the sea surface (), ; and 2) because of the absence of direct sunlight in the shadow, there is no added contribution of light between the sensor depth and , so . The two contributions are additive, such that . Assuming that the upwelling radiance attenuation is equal to a, then we obtain in Eq. (2).
The term in Eq. (2) can be substituted by the length determined by the refracted sun zenith angle in water () and the radius of the instrument (R), that is, . The general equation is written as (Gordon and Ding 1992)
where GD92 is Gordon and Ding (1992) and . Gordon and Ding (1992) performed Monte Carlo simulations to account for scattering processes in order to replace the value of 2 by a more realistic estimation of this parameter (their Table 1). From the Gordon and Ding (1992) numerical simulations for point-sensor self-shading, Zibordi and Ferrari (1995) derived a simple equation for k, that is, .
Leathers et al. (2004, hereinafter L04) revisited this model for buoyed instruments that have a different geometry from the simple one considered by Gordon and Ding (1992). In the case of a buoy, the sensor is often located at a certain depth under the buoy (), which creates the shadow (see Fig. 1 in L04). Note that in the case of a buoy, we shall use the term shadow instead of self-shadow to indicate the combined effects of both superstructure perturbation and instrument self-shading. So, the distance between the sensor and the base of the buoy must be subtracted and Eq. (2) becomes
In addition, instead of assuming , L04 used a more robust approximation of . These two differences result in the following general expression:
where is the radius of the buoy and k is now redefined as
Equation (5) is valid until the buoy completely shades the instrument. When the sun elevation is low (or high ), it is possible that the sensor itself becomes exposed to direct sunlight and then the GD92 model would apply,
L04 further proposed an analytical solution that includes the effect of light scattering. This effect may be particularly significant for large buoys when the sun is high over the horizon. Scattering increases the amount of light reaching the sensor and thus decreases the shading error. Equation (5) can be rewritten as
where is the reduced shade depth given by
where is the shaded depth obtained using , b is the total scattering coefficient, and is a system-specific factor that can be obtained from radiative transfer simulations or, as it shall be done here, determined empirically from measurements.
Finally, the shading error must be computed for the contribution of both direct and diffuse components of the downwelling irradiance, such that
3. Materials and methods
a. The Viking buoy
The DFO buoy network is deployed from May to early November due to the presence of sea ice during the winter months. The buoys are equipped with a suite of meteorological (rainfall, pressure, air temperature, humidity, and wind speed and direction—Vaisala WXT520) and oceanographic (temperature and salinity—Sea-Bird Scientific SBE 37-SI; chlorophyll-a and colored organic matter fluorescence—WET Labs WS3S and WSCD; currents—Teledyne RD Instruments 300-kHz ADCP) instruments (Fig. 1b). In addition, the buoys are also equipped with multispectral radiometers (Satlantic OCR-507R and OCR-507I) to measure in-water upwelling radiance for and above-water irradiance at seven spectral bands: 412, 443, 490, 510, 555, 670, and 683 nm. To prevent biofouling, the radiance sensor window is protected by a copper plate that rotates out of the field of view during measurements. The sensors are calibrated each year by the manufacturer prior to their deployment, or just after their recovery. The buoy (Fig. 1a) has a radius of 1.05 m (). The sensor head is installed at the center of the hull within a circular tube of 0.23-m radius () that is flushed with bromine to reduce instrument biofouling. The OCR-507R head is at 0.53 m from the base of the buoy () but adding the buoy draft puts the sensor head at 0.86 m from the water surface. With this configuration, Eq. (9) is valid at all sun zenith angles. For example, at of 70°, the shadow extends 1.06 m below the base of the buoy. At this (), the radius of the shaded surface at 0.53 m below the buoy (i.e., where the instruments are placed) is 0.52 m, which is much larger than the hull.
b. Field observations
We visited the IML-4 buoy site on 18 occasions between 3 June and 5 November 2015 on board DFO small craft boats (i.e., 30-ft length). The mooring is located in the middle of the St. Lawrence estuary (SLE) about 23 km north of Rimouski, Quebec, Canada, in 330-m water depth (inset in Fig. 1b). Light profiles were performed using a free-falling Compact-Optical Profiling System (C-OPS, Biospherical Instruments; 19 wavelengths from 305 to 780 nm). In addition, we measured in-water IOPs using an optical package, which includes a CTD probe (Sea-Bird Scientific SBE 19plus), an a-Sphere (Hobi Labs), and a HydroScat-6 (Hobi Labs). All instruments were calibrated by the manufacturers prior to the field season. Surface waters were sampled using a clean 20-L carboy kept in an iced cooler until it was processed in the laboratory later that same day. The water sample was analyzed for a number of physicochemical [suspended particles matter (SPM), dissolved organic carbon (DOC), salinity, nutrients], biological (fluorometric chlorophyll-a, pigments, flow cytometry), and optical analyses . Here we report only on the optical properties.
In 2015, the ocean color radiometer (OCR) sensors were installed only on the buoy on June 30, yielding 14 dates with concurrent radiometric measurements from the C-OPS and the moored OCR (Table 1).
1) Light profiles and Apparent optical properties (AOPs) derivation
The C-OPS can sample the light field in the upper ocean [up to nearly 0 m for and 0.35 m for ] with very high precision and accuracy at a frequency of 15 Hz. The system has been fully described in Hooker et al. (2013) and Morrow et al. (2010). The C-OPS measures and at 19 wave bands: 305, 320, 330, 340, 380, 412, 443, 465, 490, 510, 532, 555, 589, 625, 665, 683, 694, 710, and 780 nm. Simultaneous above-surface downward irradiance  was measured with the radiometer attached to the top of the boat, making sure that no obstructions were in the field of view.
Three to five light profiles were performed with the C-OPS on each date. The boat was drifting during the measurements and the instrument was kept outside any disturbance or shadow created by the boat. In general, the boat drifted within a distance 1 km from the buoy. The data were processed in R software with the package developed by Dr. B. Gentili at the Laboratoire d’Océanographie de Villefranche (LOV), respecting NASA protocols (Gentili 2016; Mueller et al. 2003). AOPs were derived from vertical profiles of upwelling radiance and downward irradiance . Each profile was carefully inspected and records showing high instrument tilt were discarded, that is, for in-water and for . In-water radiometric quantities in physical units [ and ] are normalized with respect to simultaneous measurements of the global solar irradiance, , with t explicitly expressing the time dependence, according to
where identifies the radiometric quantity that is normalized to account for irradiance fluctuations during the vertical profiles. After the normalization and filtering of , the data are fitted using a linear regression applied to log-transformed radiometric quantity versus depth, allowing extrapolation of the quantities of interest just below the sea surface () (D’Alimonte et al. 2013). In the case of upwelling radiance, the extrapolation method yields the spectral diffuse attenuation coefficient for upwelling radiance, , that is needed to calculate the water-leaving radiance from the radiance measured under the buoy at 0.86 m below the sea surface—that is,
where is the diffuse attenuation coefficient for upwelling radiance for a given surface layer. The thickness of the surface layer considered for the extrapolation varied as a function of wavelength. We begin with a 0.5-m-thick layer extending below the very first C-OPS measurements in the water column (typically around 0.4–0.5-m depth). For example, if C-OPS measurements begin at 0.4 m, then the layer extending from 0.4 to 0.9 m will be used. A linear model is fitted to the log-transformed versus z within this depth span, and the coefficient of determination correlation is calculated [Eq. (13)]. An iterative method then increases the thickness of the layer until reaches 0.99 or the thickness of the surface layer reaches 2.5 m. This method is similar to that adopted by Antoine et al. (2013). In general, the thickness of the layer considered is thinner in the UV, red, and near-infrared.
For comparison, we also computed the diffuse attenuation of downwelling irradiance, , using the same method that was used for described above. Then was corrected for instrument self-shadow following the procedure proposed by Gordon and Ding (1992) and Zibordi and Ferrari (1995) described above, with measured total absorption coefficients (see below) at each C-OPS channel. Water-leaving radiance, was calculated from the extrapolated upwelling radiance just below the sea surface as
where the factor 0.54 accounts for the partial reflection and transmission of the upwelled radiance through the sea surface. The remote sensing reflectance was computed as
For a given date, the AOPs derived from three to five profiles were averaged after eliminating spectra that showed a large discrepancy relative to the mean spectrum (i.e., if the difference between the mean and replicates was >10% in terms of ).
The fraction of diffuse skylight to the total downwelling irradiance was determined using the Biospherical Instruments BioSHADE system fitted to the reference radiometer following a simplified version of the procedure detailed in Morrow et al. (2010). Briefly, the BioSHADE is a motor that moves a black aluminum band (shadowband) that is 1.5 mm thick and 2.5 cm wide back and forth above the sensor. Under clear-sky conditions, when the shadowband completely blocks direct sun at time , the radiometer measures the diffuse skylight (minus a part of the sky that is also blocked by the shadowband), . When the shadowband is horizontal, the sensor measures the global solar irradiance. So, to assess the global solar irradiance at , we interpolate just before and after the shadowband started to shade the sensor. This allows for approximating the fraction of diffuse skylight to the total downwelling irradiance as
2) Vertical profiles of IOPs
Vertical profiles of the volume scattering function at , were measured using a Hobi Labs HydroScat-6 backscattering meter with bands at 394, 420, 470, 532, 620, and 700 nm. The instrument was factory calibrated before the field campaign following the method of Maffione and Dana (1997). Term was corrected for attenuation along the viewing pathlength of the detector using the total absorption coefficient measured in parallel with the a-Sphere (see below) (i.e., the σ correction in Maffione and Dana 1997).
The particulate backscattering coefficient is derived from as follows (Maffione and Dana 1997):
where is equal to 1.081 as provided by the manufacturer and where the contribution of pure seawater to scattering at is computed following the method of Zhang et al. (2009) using the temperature and salinity measured at the same depth with the CTD sensor.
Vertical profiles of the spectral nonwater absorption coefficient, , were measured using a Hobi Labs a-Sphere, which is a submersible Teflon integrating sphere. The instrument was factory calibrated before the field campaign with pure water. The a-Sphere measures the absorption at 1500 wavelengths between 360 and 764 nm, which are binned at 1-nm resolution. Raw data were converted into absorption coefficients using the manufacturer’s software and calibration files. Term was corrected for temperature and salinity as measured with the CTD at the same depth using the coefficients published by Röttgers et al. (2014).
Instrument calibration in air was checked after each profile following the manufacturer’s protocol. In addition, in the middle of the field season the sphere was filled with nano pure water and its absorption was recorded for 1 min. The averaged spectrum of the pure water was considered as a blank offset. In general this offset was negligible () relative to these highly absorbing waters (see results). We observed spectral anomalies in the spectra at nm after 16 September due to low battery voltage; these wavelengths were discarded from further analysis.
The vertical profiles of IOPs were binned at 1-m-depth intervals using a loess smoothing function in R (Cleveland et al. 1992).
3) Determination of shading error from in situ optical measurement
The shading error ε was determined by comparing raw to instrument self-shading-corrected [as for ], where is the mean time of the C-OPS profiles and t is the closest time of the buoy data acquisition to . The time difference between C-OPS and buoy data acquisition ranged from a few seconds to 13 min with a mean time difference of 5 min. To account for changes in downwelling irradiance between the two measurements, we normalized as
Assuming that C-OPS measurements are error free, the shading error was calculated using Eq. (1),
Note that this “measured” ε is not purely a shading error but also includes differences impeded by the radiometric calibration and spectral responses of each instrument. In addition, uncertainties in ε may be due to variations in solar illumination resulting from scattered clouds during measurements. Poor illumination conditions were observed on 6 and 24 October, where the coefficient of variation of the downwelling irradiance at 555 nm during the 2-min recordings of the buoy reached 19% and 7%, respectively. For this reason, those dates were not considered for further assessment and algorithm development.
4) Laboratory analyses
(i) CDOM absorption
Seawater was filtered under low vacuum on 47-mm glass fiber filter (Whatman, GF/F 0.7-mm nominal pore size), which was pre-ashed for 2 h at 450°C and pre-rinsed with 50 mL of MilliQ water. Colored dissolved organic matter (CDOM) absorbance was measured with a Perkin Elmer double-beam Lambda-850 spectrophotometer using a 10-cm quartz cell between 220 and 800 nm against nano pure water. Absorption coefficients were then calculated according to
where is the absorption coefficient of CDOM (m−1) at wavelength λ, is the absorbance measured at wavelength λ, and l is the pathlength of the optical cell (0.1 m).
(ii) Particulate absorption
Measurements of particulate absorption were performed using the filter-pad technique. A known volume of the water sample (in two to five replicates) was filtered through Whatman GF/F glass fiber filters shortly after sampling (3 h). Each filter was then placed in the center of a 150-mm integrating sphere equipped with a spectralon filter holder (see Röttgers and Gehnke 2012 for technical details). The optical depth of the particles retained on the filter was then measured using a PerkinElmer Lambda-850 spectrophotometer, from 300 to 800 nm at 1-nm resolution. Optical depth was converted to the spectral particulate absorption coefficient, , using
where is the optical density of a blank filter, A is the clearance area of the particles on the filter (m2), V is the volume of sample water filtered (m3), and β is the pathlength amplification factor. The relationship between β and OD derived experimentally by Stramski et al. (2015) was used.
After the OD scanning, phytoplankton pigments were extracted for 18–24 h using methanol (Kishino et al. 1985), which removed nearly all pigments (95% of sample). The filter was then placed again in the integrating sphere to measure the absorption coefficient of nonalgal particles, . The absorption coefficient of phytoplankton was obtained by subtracting from the total particulate absorption coefficient.
(iii) SPM and SPIM
Known volumes V in liters of seawater (1 L, depending on turbidity) were filtered in triplicate through pre-ashed (1 h at 450°C) and pre-weighed (mass , mg) 25-mm glass fiber filters (Whatman, GF/F 0.7-μm nominal pore size) at low vacuum (Van Der Linde 1998). Each filter was then rinsed with Milli-Q water, dried for 24 h at 60°C prior to weighing under a dry atmosphere (mass , mg) to obtain the SPM concentration (mg L−1) as
Organic matter lost on ignition (LOI) was determined after baking the filters for 3 h at 500°C, weighing again, giving the concentration of suspended particulate inorganic matter (SPIM). The average coefficient of variation for triplicates was 11% and 14% for SPM and SPIM, respectively.
(iv) Chlorophyll-a (Chl-a)
Duplicate subsamples (100 ml) for Chl-a determination were filtered onto Whatman GF/F. Chl-a concentrations were measured using a Turner Designs 10-AU fluorometer, following a 24-h extraction in 90% acetone at 4°C in the dark without grinding (acidification method: Parsons et al. 1984).
a. Optical and biogeochemical variability at the IML-4 buoy site
Compared to oceanic waters, estuarine waters of the St. Lawrence can be considered as moderately turbid ( mg L−1), CDOM rich ( m−1), and mesotrophic to eutrophic ( mg m−3) (Table 2). On average, CDOM, phytoplankton, and nonalgal particles account for 64.8%, 22.5%, and 12.7%, respectively, of the light absorption at 443 nm.
Spectral absorption measured at the surface using the a-Sphere displayed an exponentially increasing shape toward short wavelengths that is typically found in estuarine waters (Fig. 2a). In these brackish waters, with salinity varying from 25.79 to 29.3, CDOM dominates the total absorption at 500 nm. The presence of phytoplankton pigments is significant, as illustrated by the presence of an absorption peak at 676 nm, which almost reached the level of pure water absorption on a few occasions (15 and 26 June). Nonalgal particles make a significant contribution to , with being often greater than the phytoplankton absorption coefficient at 500 nm (not shown).
Particulate backscattering, , varied between 0.0042 and 0.028 m−1. For comparison, the minimum at IML4 was larger than the maximum value reported by Antoine et al. (2011) at the BOUSSOLE site in the Ligurian Sea. The spectra were relatively flat with values from 0.41 to 0.81, with an averaged value of 0.57 0.12, for the exponent of the power-law function that describes the spectral dependency of (Fig. 2b).
Figure 3 shows the apparent optical properties measured at IML-4. Diffuse attenuation of upwelling radiance spectra is similar to that of total absorption (Fig. 2a vs Fig. 3a), only with larger amplitude due to the variability of . For comparison, the averaged spectrum of (black line) and (gray dashed line) are shown. Both were almost identical except at 600 nm, where was systematically lower than . At 683 nm, for example, and were on average 0.45 and 0.68 m−1, respectively, whereas pure water absorption at this wave band is 0.48 m−1. The relatively low values in the red part of the visible spectrum indicate, as expected, that inelastic scattering processes (Raman and sun-induced chlorophyll-a fluorescence) make an important contribution to upwelling radiance in this spectral region (Li et al. 2016).
The remote sensing reflectance typically peaks in the green wavelengths (532, 555, and 589 nm), where the variability and the dynamic range are greater than at other wavelengths (Fig. 3b).
b. Measured shading error
The measured shading error (%) of the Viking buoy was 63% 13% (average ± standard deviation), but as expected it varies with wavelength following the spectral total absorption coefficient (Fig. 4). Measured ε values were 47% 8% and 82% 6% at 555- and 412-nm channels, respectively. These values are much higher than for other commercially available buoyed instruments, such as the Satlantic tethered spectral radiometer buoy (TSRB) for which shading error correction methods have been proposed (Leathers et al. 2001). For a given channel, however, the dynamic range of ε did not exceed 30%.
c. Modeling shading error
1) Semianalytical models
As mentioned above, L04 formulated semianalytical models to predict the shading error for radiance sensors maintained below a larger buoy frame. The first model neglects the scattering processes [Eq. (5)], such that ε depends only on the sun elevation, the fraction of diffuse sky irradiance (f), and the total absorption coefficient, . Figure 5 shows the predicted ε for both direct (gray curves) and diffuse (dashed curve obtained using of 35°) components of the downwelling irradiance as a function of for the Viking buoy geometry. Most observed ε, which were acquired with varying between 31° and 64° (Table 1), fall within the predicted range. Note that ε at 670 and 683 nm are lower than those at 412 and 443 nm for similar , as expected in the presence of sun-induced fluorescence by phytoplankton pigments.
We applied the L04 models to our dataset using measured IOPs and the modeled fraction of the total downwelling irradiance that is skylight (). The latter was computed using the Gregg and Carder (1990) irradiance model at the buoy data acquisition time. The atmospheric visibility was changed in order to get modeled equal (within 5%) to the measured irradiance at that wavelength. Under stable atmospheric conditions, modeled and measured with the BioSHADE were similar (within 10%, not shown).
Spectral absorption and particulate backscattering from the a-Sphere and HydroScat-6, respectively, were averaged within the first 4 m of the water column. For a few dates, the a-Sphere data were not available (N = 3; Table 1), either due to nondeployment or instrument failure. For these dates, discrete measurements of and were used instead. For dates when HydroScat-6 data were not available (N = 5), we derived from the C-OPS-measured using a tuned version of the quasi-analytical algorithm (QAA) of Lee et al. (2013).
To apply the semianalytical model that includes scattering [Eqs. (9) and (10)], was converted to assuming a backscattering ratio of 1.8% and the pure water scattering coefficient calculated following Mobley (1994) was added. Next we estimated the value of of Eq. (10) for our system by minimizing the root-mean-square error (RMSE) and the bias calculated between predicted and measured ε (here = 0.17):
Figure 6 shows the comparison between modeled and measured ε, while statistics are presented in Table 3. The semianalytical model that neglects scattering [Eq. (5)] tended to overestimate the shading error at all wavelengths (Fig. 6a). When scattering was included in the semianalytical solution [Eqs. (9) and (10)], the positive bias on the predicted ε was removed, but the data points were more scattered (Fig. 6b), as shown by the decreased (from 0.73 to 0.67). Interestingly, the residuals of both models were not spectrally flat (Figs. 6a,b and 8). The ε residuals were significantly (p 0.05) correlated to for both semianalytical models (Fig. 7), with an of 0.40 and 0.33 without and with scattering included, respectively. In the red channels, however, the semianalytical models clearly overestimated ε. Again, this is due to sun-induced fluorescence by phytoplankton that allows more photons to reach the shaded sensor. Similarly, at the shortest wavelengths (i.e., 412 and 443 nm), sun-induced CDOM fluorescence may explain part of the ε overestimation.
2) Empirical model
With the aim of developing a shade correction algorithm to estimate from the IML-4 buoy OCR data, empirical relationships between the parameters measured by the buoy instruments and ε were explored. This empirical approach at the very least will be applicable to our study site. First, the uncorrected upwelling radiance measured by the sensor under the buoy was normalized by the measured downwelling irradiance above the sea surface (), calculated as
Visual comparison between the spectral shape of and C-OPS showed that both spectra shared striking similarities, suggesting that even without any correction, the measured radiometric quantities potentially held useful information on the water column optical properties.
A multiple linear regression (MLR) analysis was thus performed to predict ε at each OCR channel using band ratios of , sun zenith, angle and the fraction of diffuse sky irradiance. Band ratios () included , , , and , and the ratio proposed by Lee et al. (2013) for the most recent QAA version,
The general MLR equation used was
where are MLR coefficients, BR is one of the band ratios mentioned above, and f is the fraction of diffuse to total downwelling irradiance above seawater [see Eq. (16)]. The best set of predictors (i.e., f, BR, ) were selected based on the highest adjusted (Table 4). The variable BR was included for all channels, but f was significant at 443 and 683 nm and at 412, 670, and 683 nm. The ratio was the best predictor of ε at all wavelength except at 412 nm, where provided the best performance. In general, the low values of adjusted were due to the relatively small dynamic range in ε for a given channel (Fig. 4) and the small number of observations (N = 11; after excluding dates where poor illumination conditions were observed and where was ).
As expected, the locally tuned empirical model performed better than the analytical models (Fig. 6c). The RMSE were half of those obtained with the analytical models and the modeled ε residuals were spectrally flat (Fig. 8c).
d. Correction scheme for autonomous remote sensing reflectance retrievals
A correction scheme was developed to estimate the spectral remote sensing reflectance from the buoy radiometric measurements. The scheme includes three steps.
Step 1: The empirical models presented above were employed to estimate from and to correct the measured upwelling radiance,
Step 2: The diffuse attenuation coefficient of upwelling radiance was estimated to calculate the water-leaving radiance. Empirical relationships were established between the mean diffuse attenuation for upwelling radiance within the surface layer, , and BR and/or . The general MLR applied to the data is
where are empirical coefficients presented in Table 5. Figure 9 shows the scatterplot of the measured versus retrieved using the relationships. The best band ratio for 412 and 443 nm was ξ, while the ratio was selected for other wavelengths. The determination coefficient tends to decrease from 0.79 to 0.34 with increasing wavelengths (Table 5). This decrease in toward the red is due to the smaller dynamic range of at longer wavelengths (Fig. 3a), while inelastic scattering may also explain the low obtained at 683 nm. The root-mean-square log error (RMSLE) varied between 0.037 to 0.057 at 555 and 683 nm, respectively. Term was used to estimate the upwelling radiance just below the sea surface as
e. Application of the correction scheme to the IML-4 buoy
The performance of the correction algorithm was evaluated based on six statistical values: the percent bias, the absolute percent relative difference (AD), the RMSE, the , and the slope and intercept of the linear regression (type 2) between C-OPS and the buoy estimates. The results are presented in Table 6 and Fig. 10.
A strong linearity is observed between measured and retrieved with a and a slope close to unity (Table 6). Overall, was retrieved with an AD ranging from 12% to 20%. Negative bias was observed at 412, 443, 670, and 683 nm. This is likely due to the empirical models that underestimated both ε and , which cumulated and led to the larger underestimation of .
We further partitioned the uncertainty in resulting from the ε and corrections. The was estimated using the true and together with the empirically derived from the buoy data. This is equivalent to assuming a perfect shadow correction of . The results (dots in Fig. 10 and numbers in parenthesis in Table 6) clearly show a major improvement in the retrieval, with very strong linearity between measured and retrieved ( at all λ). AD decreases by a factor of 3 at all wavelengths. These results globally suggest that the uncertainty in is driven more by the shadow correction than by the attenuation-based correction.
Fig. 11 shows six examples of spectra retrieved from the buoy as compared with the measured spectra obtained by the C-OPS. In the worst cases (Figs. 11a,f), the spectral shape of the retrieved was still well preserved. The large bias observed on 24 October (Fig. 11f) was expected, since it was obtained under poor illumination conditions (see above). As a result, the reflectance band ratios normalized by were retrieved with better accuracy than the absolute (Table 7; Fig. 12). In fact, the AD and BIAS were 3% for and , and 10% for and .
a. Uncertainty of the buoy’s retrieval
Seventeen years after the launch of the new generation of satellite ocean color sensors (i.e., SeaWiFS), there has been no systematic validation of the ocean color primary products, such as the remote sensing reflectance (or ) in the EGSL. The lack of field observations has been a major limitation to validating satellite-derived and to proposing new or improved algorithms. The network of buoys deployed by DFO every summer, which are all equipped with OCR, represented a great opportunity to provide the needed in situ data to achieve satellite data validation goals. However, we showed that the radiometric data they provide require corrections due to 1) the shading effect caused by the buoy and 2) the diffuse light attenuation between the depth at which is measured and the sea surface.
The correction scheme presented here is a first attempt to retrieve from the buoy data in an autonomous mode. Since it is empirically based, its application to the other mooring sites with different optical properties may not be appropriate. However, since IOPs are unknown a priori at the buoy sites in real time, the application of analytical models of shadow correction (Gordon and Ding 1992; L04) and diffuse attenuation (Lee et al. 2005) in an autonomous mode is not trivial. To apply an analytical or semianalytical approach in an operational acquisition mode, one would need to derive the IOPs from the parameters measured by the buoy instruments in real time. In the current setting, in addition to the radiometric measurements, only in vivo fluorescence instruments for Chl-a and CDOM are available. Poor relationships were found between fluorescence variables and IOPs and therefore were not used to predict the shadow effect. Adding a backscattering meter would provide , but the most important parameter to know is the absorption coefficient. Analytical models would also need to be refined to include sun-induced fluorescence processes, which are particularly evident in the red channels.
Attempts were made to retrieve IOPs using a QAA. After an initial guess for ε and to obtain an initial guess of , IOPs were retrieved using a tuned version of the QAA (Lee et al. 2013). Then an iterative scheme in which and are used to estimate ε with analytical models [i.e., Eqs. (5) and (9)] and with a semianalytical model for (Lee et al. 2005) was implemented. In general, this approach failed to converge toward a satisfactory solution and introduced large artifacts in the final . One reason explaining the failure of this approach is the fact that Lee et al.’s (2005) model relates IOPs to for the surface (0–1-m depth) rather than . Our in situ measurements of and indicate that both quantities were nearly equal in the blue and green wave bands (412–555 nm), but not in the red channels, where due to inelastic scattering (Fig. 3). In addition, retrieved using the QAA was not perfect, despite a relatively good retrieval (i.e., absolute percent relative difference 12% and , not shown). Another approach would be to develop an inversion model based on precomputed lookup tables (LUTs) derived from radiative transfer simulations, including inelastic processes. For example, Huot et al. (2007) developed a method based on LUTs, but their mooring was equipped with sensors at two depths, allowing a direct estimation of . The DFO’s moorings, however, rely on only a single measurement at 0.86 m.
The empirical algorithms developed here are by definition limited to the range of optical (Table 2) and illumination conditions encountered at the buoy site. Fortunately, our field observations cover a full summer season and half of an autumn season with fairly good range of sun elevation (Table 1). One limitation is the fact that the lowest encountered was 31°, which is higher than the actual at the time of the satellite overpass near the summer solstice (i.e., 26°). This is because our field observations were made at the end of the morning [around 1000–1100 local time (LT)], while the satellite overpasses are closer to noon. Under clear sky, ε increases rapidly with decreasing as the sun get higher in the sky (Fig. 5). Therefore, more observations around noon under clear sky would be required to improve the predictive capability of our ε empirical model.
The uncertainties estimated here are also large because the measured ε is not error free. In situ determination of ε is uncertain due to several reasons:
The location of the light profiles is not always exactly beside the buoy. As the wind blows, the boat can drift for several hundred meters away from the buoy while the profile casts are made. Because of its large size, some spatial heterogeneity is expected in the St. Lawrence estuary.
Differences in the sensor spectral response between OCR and C-OPS (neglected here).
Vertical gradients in IOPs were sometime strong in the top 10 m, resulting from thin fresher water layers.
More data would be needed to fully quantify the uncertainty of our approach. In fact, with only 11 days of concomitant observations that pass quality control, it was not possible to split the dataset into two subsets: one to derive the empirical coefficients and one to validate. Therefore, the metrics used to evaluate the model (bias, RMSE, , slope, intercept) may actually be better than in reality.
b. Sample results from the IML-4 buoy
The SLE is characterized by cold, nutrient-rich estuarine waters showing moderately turbid, CDOM- and phytoplankton-rich conditions. It differs markedly from other coastal validation sites of AERONET-OC (Zibordi et al. 2009). The presence of mineral particles (SPIM = 3.77 2.04 mg L−1; Table 2) and high makes these waters relatively bright at long wavelengths. A preliminary assessment of the SeaWiFS data quality reported by Le Fouest et al. (2006) suggested that standard atmospheric correction (AC) algorithms would fail to retrieve at short wavelengths in these waters. Recent improvements in turbid waters AC (Bailey et al. 2010; Wang and Shi 2007; Lavender et al. 2005; Ruddick et al. 2000) could increase ocean color data quality retrieved over the estuary, but this remains to be quantitatively assessed. The DFO’s mooring network may have the potential to fulfill this important task.
The new Viking buoy design was deployed at IML-4 (Fig. 1) for the first time in 2013 and every year since then. Some measured meteorological and oceanographic parameters are depicted in Fig. 13. One of the problems in obtaining accurate radiometric measurements from moored buoy in relatively deep waters (here 330-m depth) is the buoy tilt (Antoine et al. 2008). The median buoy tilt observed over the 3 years of deployment was only 4° with 98% of the observations having a tilt <10°. The higher tilts observed in late summer 2014 are explained by the addition of a small winch on one side of the buoy, changing the center of gravity. In 2015 the buoy was slightly modified to correct this problem. Additional weight near the center of the buoy reduced the tilt with 88% of the observation being <5° and a mean value as low as 3.8° (cf. to 4.9° in 2013–14). Therefore, unlike BOUSSOLE (Antoine et al. 2008) the buoy tilt is not a major problem for the Viking buoy, especially the most recent design that includes a winch.
Figures 13c,d show that there is a very strong short-term variability of temperature and salinity. For example, changes as large as 5 psu in salinity and 10°C in temperature can be observed within less than 24 h. Variability of optical properties at the IML-4 mooring site is thus expected to be significant at short time scales, making satellite ocean color validation challenging using punctual ship-based observations.
Figures 13e and 14 show the in situ Chl-a concentration (in vivo fluorescence sensor) and estimated corresponding to the median value of all buoy measurements taken between 1000 and 1400 LT. The in situ chlorophyll-a concentration determined in water samples (in vitro fluorescence) and C-OPS-, MODIS-Aqua-, and VIIRS-derived reflectances (red, blue, and green symbols, respectively) are also shown. For the satellite data, only those measurements having met the quality criteria are plotted (Bailey and Werdell 2006).
Although a full assessment of the current satellite products of MODIS and VIIRS is out of the scope of this study, a preliminary assessment of the satellite retrieval of both and Chl-a highlighted interesting features. First, it is clear from Figs. 14a,b that the retrieval of the at short wave bands (i.e., 412 and 443 nm) tends to be largely underestimated relative to in situ data with frequent negative values in summer. MODIS-Aqua appears to be more problematic than VIIRS with 53% versus 38% of the (412 or 410 nm) observations being negative, respectively. Second, the retrievals in the green bands (551 or 555 nm) fall within the lower range of the in situ data (Fig. 14c). Third, the in the red channels (670 nm) is better retrieved than that at shorter wave bands (Fig. 14d). Overall, these results indicate either atmospheric correction problems or inaccurate sensors calibration. The former is more likely in the presence of moderately turbid waters, such as those of the St. Lawrence estuary. We hypothesize that the standard atmospheric correction scheme overestimates the spectral dependency of the aerosols, leading to a larger underestimation of in the blue channels relative to the red channels. Another consequence of the spectrally dependent underestimation of is observed in the retrieved blue-to-green reflectance ratio, which is often underestimated (Fig. 14e). As a result, Chl-a will be overestimated using any standard band ratio algorithm used by the space agencies. Therefore, the overestimation of the Chl-a depicted in Fig. 13e probably results from the combined inaccuracy of the blue bands’ retrievals and the inadequacy of the globally tuned empirical Chl-a algorithm.
Fisheries and Oceans Canada has developed a network of buoys moored at strategic sites that are part of an extensive monitoring program. The optical data provided by this network offers great potential to support space-based studies of environmental variability and to track changes in the marine environment under increasing anthropogenic pressures. The results showed that in situ can be retrieved with an averaged absolute percent difference ranging from 12% to 20%. Most of the difference results from the shadow correction (65% of the difference), while the estimation of water-leaving radiance based on a single measurement of at 86 cm from the sea surface accounts for 35% of the difference. This level of uncertainty is larger than the 5% accuracy target of the space agencies (Clark et al. 2003). Considering the systematic bias of current satellite-derived that likely results from the aerosol selection during the atmospheric correction process, a 20% accuracy is already a major improvement that can lead to better validated satellite-derived products over the complex EGSL waters.
New developments are foreseen to improve the radiometric data quality from the Viking buoys. First, installing the OCR sensors at the extremity of extension arms on the buoy, such as those of MOBY or BOUSSOLE, would significantly reduce the shadow error. An alternate solution tested in fall 2016 consists of installing a second sensor on a smaller buoy (diameter of 28 cm) tethered to the main buoy to dramatically reduce the instrument self-shading. Second, the addition of deeper sensors could be used to directly estimate , further reducing the uncertainty of in situ . Considering the strong light attenuation in the EGSL, the second set of sensors would not need to be as deep as those of MOBY or BOUSSOLE (9 m). Based on the attenuation coefficient measured at the IML4 site, a second sensor between 2- and 3-m depth would allow for good estimation of . Considering that inelastic scattering processes may be significant in the red part of the spectrum (Li et al. 2016), resulting in abrupt changes in within the first meters of the water column, we recommend putting the sensors as close to the sea surface as possible. Third, a backscattering meter may improve the empirically based shadow correction. Finally, as pointed out by Zibordi et al. (2015), hyperspectral radiometric measurements should be considered in order to increase the applicability of the in situ for the CAL/VAL purpose of ocean color sensors with different spectral bands’ responses.
Retrieving Chl-a or IOPs from space in the estuarine waters of the EGSL is challenging. The use of semianalytical or quasi-analytical algorithms should be preferred in these optically complex waters (Lee 2006). However, these methods would fail as long as retrievals are imperfect (e.g., negative reflectance in the blue/violet bands). More work must be done to improve atmospheric correction over the EGSL or to develop local vicarious calibration coefficients that would improve the quality of the satellite retrievals. Using the buoy data, a preliminary assessment of the retrieved by the recently launched OLCI sensor [on board the European Space Agency (ESA)’s Sentinel-3 satellites)], which has more spectral bands in the near-infrared for atmospheric correction, is promising.
We are grateful to Gabriel Ladouceur, Frédéric Diotte, and Roger Pigeon for their assistance in the fieldwork, and to Jacques St-Pierre from Multi-Electronique for providing the buoy geometry. We thank the NASA OBPG for providing the satellite ocean color data freely. We thank David Antoine, Giuseppe Zibordi, and one anonymous reviewer for their valuable advice on the manuscript. This work was supported by an NSERC Discovery grant and a Université du Québec à Rimouski grant to Simon Bélanger, and a CSA grant to Pierre Larouche for the buoy development.