1. Introduction

Analytical solutions to determining water quality from optical remote sensing commonly face three unknown parameters. These are (1) the concentration-specific optical properties (absorption and scattering) of each optically active substance in the water (e.g. phytoplankton pigments, cell material and other organic as well as inorganic suspended solids, coloured dissolved matter), (2) the concentrations of said substances and (3) their individual distribution with depth. Reciprocally, to diagnostically estimate the concentrations of optically active water constituents, several criteria must be met: the absorption and scattering properties of the target substances should be a priori characterized, these properties should not vary over the depth layer observed from remote sensing, and the depth distribution of the substances should be uniform (or at least known) over the depth layer observed by the remote sensor.

Continuous improvements of remote sensing technologies and computing over the last decades have supported a wide range of algorithmic approaches to estimate optically active substance concentrations. Improved multispectral waveband configurations and sensitivity, particularly in the near infra-red (NIR), have advanced water quality remote sensing in optically complex waters. In all but the most turbid water systems, water itself increasingly dominates absorption with increasing wavelength into the NIR. Inversion of subsurface reflectance into spectral absorption ($a(\lambda )$) and backscattering ($b\mathrm {_b(\lambda )}$) coefficients can exploit this behaviour, under the assumption that (1) absorption by other substances can be assumed negligible, (2) scattering by hydrosols is without significant spectral features, and (3) vertical mixing is uniform over the visible water column. This set of conditions allows $b\mathrm {_b}$ to be isolated and extrapolated to the visible light spectrum (with increasing uncertainty towards shorter wavebands). Applications of this inversion mechanism include widely used examples such as the Quasi-Analytical Algorithm [1] and Garver-Siegel-Maritorena semi-analytical algorithm [2,3].

The assumption of vertical homogeneity is not without problems. It can break down in cases of shallow thermal or density stratification or when particle distributions form non-uniform layers in the ocean [4]. These non-uniform distributions can only be corrected for in global biogeochemical models if an appropriate distribution profiles can be assumed [5]. Shallow thermal or density stratification is more likely to occur in coastal and inland waterbodies, which tend to have higher optical complexity, such that generalisations on vertical distributions are less likely to hold. In addition, motile or buoyant phytoplankton may form discrete layers, either near or away from the surface [6]. Under any of these conditions, estimates of optically active substance concentrations derived from remote sensing are likely to deviate from their true values. Observation uncertainty is associated with the optical ambiguity of the remotely sensed reflectance: a thin optical layer with high absorption and scattering is readily confused with a thicker homogeneous layer having smaller absorption and scattering coefficients. Thus, in the absence of a priori information on vertical structure, the quality of the estimate will be strongly influenced by the assumption of either a specific vertical gradient or uniform mixing over the observed depth.

It has been shown [7] that depth distribution and optical properties are unambiguously related in a non-homogeneous water column, such that bulk inherent optical properties retrieved through reflectance inversion are equivalent to the average optical properties over the light penetration depth [7,8]. Knowing the depth distribution of the inherent optical properties and relating these to the depth from where the reflected signal originates, further improves optical closure [9]. While the spectrally variant bulk absorption properties of the water column predominantly determine the curvature of the reflectance spectrum, the amplitude of reflectance is related to light penetration depth [10] as it is modulated by the efficiency of light scattering relative to absorption. While still under the assumption of uniform mixing, the light penetration depth itself may be estimated from reflectance [11].

When light penetration depth exceeds vertical mixing depth, some spectral regions are likely to show larger inconsistencies between reflectance amplitude and the depth-averaged inherent optical than others, based on how much of the light beam was already attenuated in the overlying layer. Observing such spectral inconsistencies in bulk absorption or scattering should be cause for caution to further attempt interpreting the reflectance spectrum, at least in terms of substance concentrations. These inconsistencies can only be identified when their depth-weighted contribution to reflectance is already known, which will not be the case in typical remote sensing scenarios. The well-described absorption and scattering properties of water itself form the only exception: they are predictable across the visible to near-infrared light spectrum and, for current purposes, unlikely to vary greatly with depth. The shape of the reflectance spectrum (under vertically homogeneous conditions) is therefore predictable in parts of the spectrum where water dominates the absorption, and in the absence of significant curvature introduced by scatterers. This has been demonstrated for clear to relatively turbid waters in the NIR part of the spectrum [12]. In the visible domain, increasing overlap between the absorption features of water constituents and those of water itself present a different challenge. Small variations in water absorption are, however, detectable throughout the visible and near infrared, associated with different vibrations and overtones of the water molecule. These features are narrow and not well captured by the most common wavebands on multispectral ocean colour instruments on satellites.

Hyperspectral reflectance measurements extend the state-of-the-art in aquatic remote sensing by offering more degrees of freedom in inverse optimization problems designed to derive optically active substance concentrations. Hyperspectral bio-optical model-based inversion using a priori assumed absorption and scattering models have been used widely for optically complex waters to resolve the individual optical components (e.g. [13–15]). Here, we hypothesize that it is possible to identify specific regions of the reflectance spectrum that are primarily shaped by light absorption by water. The inversion of reflectance in these wavebands then yields a value distribution of $b\mathrm {_b}$, and the shape of this distribution may be analysed to assess consistency in retrieval across the spectrum, potentially bringing conditions of vertical non-uniformity to light. Ideally, practical thresholds for this consistency can then be established so that suspect observations, possibly related to the vertical distribution of optical components in the water column, can be identified.

This work presented here details a waveband selection concept from hyperspectral reflectance, aimed to improve confidence with respect to the assumption of uniform vertical mixing in optical remote sensing. The procedure retrieves the bulk inherent optical properties and exploits the widest possible range of the reflectance spectrum to determine a distribution of $b\mathrm {_b}$ values, compared to conventional NIR band based approaches. We present this concept using a combination of simulated reflectance for water columns with varying vertical uniformity and a wide range of in situ observations from optically complex inland water systems covering oligotrophic to hypereutrophic systems.

2. Methods

2.1 Band selection and inversion procedure

The proposed procedure consists of six steps which are illustrated in a flowchart (Fig. 1) and in four examples of in situ data included in the Results section (Figs. 3–6). These steps are detailed in this section.

 

Fig. 1. Steps taken during (1-2) pre-processing of reflectance and $a\mathrm {_w}$ spectra, (3) spectral matching of features in reflectance and $a\mathrm {_w}(\lambda )$, (3-4) selecting wavebands where $a\mathrm {_w}(\lambda )$ strongly influences the shape of the reflectance spectrum, (5) reflectance inversion, and (6) analysing the distribution of retrieved $b\mathrm {_b(\lambda )}$ values. Each of the numbered steps are detailed in Section 2.1

Download Full Size | PDF

Step 1 rescales the amplitude of reflectance ($R$) and water absorption (${a\mathrm {_w}}$) spectra linearly between 0 and 1 (further noted as $\hat {R}$ and $\hat {a}\mathrm {_w}$). Because of this normalization between minimum and maximum amplitude, any form of reflectance may be used (i.e. subsurface or above-surface, irradiance or radiance reflectance) as long as any effects of bidirectional reflectance at the air-water interface, sun-glint, white caps, floating vegetation or scums masking the water absorption signature, or other disturbances can be assumed negligible and without spectral dependencies. We also note that the case of optically shallow water (bottom reflectance) is not considered here. Reflectance must further be sufficiently spectrally resolved and free of noise to evaluate their shape over multiple relatively narrow (10-nm) intervals. In Step 2, both $\hat {R}$ and $\hat {a}\mathrm {_w}$ spectra are locally smoothed (LOWESS algorithm from statsmodels v.0.9.0 for Python 2.7) over a 10 nm window to reduce noise while maintaining general spectral features. Wherever spectral intervals are indicated in these processing steps, any wavebands with a centre wavelength within the stated interval are included. The $a\mathrm {_w}$ spectrum used here was obtained from [16] adjusted to 15$^{\circ }$C and zero salinity which is assumed representative of sampling conditions in our in situ data set (described in Section 2.2), and interpolated to the spectral resolution of the $R$ of each specific sensor following the above procedures. The strongest effects of temperature on water absorption, when observed, may be avoided by omitting the 750 nm spectral region; this was not done here to allow for any prominent effects to show in the results.

Step 3 locates wavebands where the shape of $R$ is strongly influenced by the shape of $a\mathrm {_w}(\lambda )$, which is done in two parts (Steps 3A-B in Fig. 1). We first determine the ratio of reflectance wavebands sampled along the spectrum over set intervals, which characterizes the changing slope along the spectrum. The along-spectrum ratio of $\hat {R}$ and inverse along-spectrum ratio of $\hat {a}\mathrm {_w}$ at a given waveband $\lambda$ (in nm) are defined, respectively, as the dimensionless $\hat {R}^\prime (\lambda )$ and $\hat {a}^\prime \mathrm {_w}(\lambda )$:

Wavebands, where the difference between $\hat {R}^\prime$ and $\hat {a}^\prime \mathrm {_w}$ are within a specified threshold, are taken forward. A tolerance of $\pm$5 % between the along-spectrum ratios was empirically determined to distinguish between similar and dissimilar spectral regions in our test data sets (see Section 2.2). To determine tolerance thresholds, $\hat {R}^\prime$ and $\hat {a}^\prime \mathrm {_w}$ are assumed to be compared on similar scales and units of measurement, which is valid as long as the amplitude of $R$ is approximately proportional to the absorption coefficient ($R\approx {a}^{-1}$). In the simplest case, for clear natural waters, subsurface irradiance reflectance is proportional to the ratio $b_\mathrm {b}$ / $a$, and because $b\mathrm {_b}$ « $a$ this assumption holds. For increasingly turbid waters, the relation between $R$ and the inherent optical properties is better approximated as proportional to $b\mathrm {_b(\lambda )}$ / ($b\mathrm {_b(\lambda )}$ + $a(\lambda )$), while $b\mathrm {_b(\lambda )}$ « $a(\lambda )$ may still be valid in parts of the spectrum, particularly at longer wavelengths where water absorption dominates the inherent optical properties [17]. As an optional quality control step, omitting wavebands where the reflectance spectrum is noisy can be helpful to remove outlier results including those which may fall within the stated threshold. For results shown here, noise was removed where the coefficient of variation within directly adjacent wavebands of $\hat {R}$ exceeded 1, although instrument dependent thresholds and wavelength intervals may be considered.

The ratio values $\hat {R}^\prime (\lambda )$ and $\hat {a}^\prime \mathrm {_w}(\lambda )$ may intersect in areas where they diverge rather than align, which is an indication that the shape of $R$ is influenced by another optical constituent than $a\mathrm {_w}$ in the given spectral region. These divergent areas should be excluded and may be removed by comparing the derivative spectra of $\hat {R}^\prime (\lambda )$ and $\hat {a}^\prime \mathrm {_w}(\lambda )$ (Step 3B). The first derivative spectra were computed by convolving the curves over a window of $\pm {}$ 5 spectral channels (noting that this implies broader windows for sensors with coarser resolution, but narrow windows increase noise), using the signal convolution function of scipy v1.5.4. Wavebands, where the derivative spectra of $\hat {R}^\prime (\lambda )$ and $\hat {a}^\prime \mathrm {_w}(\lambda )$ differed by more than $\pm$0.025 nm$^{-1}$, were omitted from further analysis here. This threshold value retained multiple adjacent wavebands per area of spectral convergence, while omitting the areas in between such features, i.e. where $a\mathrm {_w}(\lambda )$ does not determine the shape of $R(\lambda )$.

In relatively clear waters, a decrease in $R(\lambda )$ can often be observed around 605 nm. This is associated with the fifth O-H harmonic stretch vibration of the water molecule [18] and a doubling of $a\mathrm {_w}$ over a 20-nm interval. A similar change in $R(\lambda )$ can be caused in more productive waters by accessory photosynthetic pigment absorption, notably chlorophyll-c or phycocyanin, in the 615-625 nm region. In such cases, misinterpreting the shape of $R(\lambda )$ as being caused by $a\mathrm {_w}$ alone would lead to strong underestimation of $b\mathrm {_b(\lambda )}$ in this region compared to other (longer) wavelengths where $a\mathrm {_w}$ typically dominates absorption. A simple workaround (Step 4) is to detect the presence of phytoplankton absorption in the red spectral region. We apply a ’red edge’ detection to determine whether phytoplankton absorption is evident at wavelengths < 700 nm. If $R(700)$ exceeded $R(675)$ by 10 %, only the region > 700 nm was considered. In the case where $1.1\times R(675) > R(700)$, no red edge is detected and the wavelength interval is considered where $a\mathrm {_w}$ > 0.1 m$^{-1}$ (i.e., wavelengths > 583 nm for pure water at 15$^{\circ }$C).

Having identified a set of wavebands where $a\mathrm {_w}$ defines the shape of $R$, it is now possible to invert the original reflectance values (before normalization) at these wavebands into the bulk inherent optical properties (Step 5). Given that $a\mathrm {_w}(\lambda )$ must dominate absorption for the change in $R(\lambda )$ to correspond to the change in $a\mathrm {_w}(\lambda )$ over the small spectral intervals considered, and with $b\mathrm {_b(\lambda )}$ unknown, the solution for the latter only depends on how reflectance was recorded (above vs. below-surface). Rearranging the reflectance inversion model from [19] for subsurface radiance reflectance $R_L(0^-, \lambda )$, the backscattering coefficient is retrieved for a waveband where $a\mathrm {_w}$ dominates total absorption, as (see e.g. [20,21]):

Following these steps, a number of $b\mathrm {_b(\lambda )}$ wavebands will be available, located between the yellow-orange (for very clear waters) to red and near infra-red spectral region. Each waveband corresponds to a particular sunlight penetration depth, depending on variant volume absorption and scattering properties. In homogeneously mixed waters we should expect to observe only minor variations in $b\mathrm {_b}$ between wavebands where absorption is dominated by $a\mathrm {_w}$, while allowing for a gradual decrease with wavelength depending on the shape of the particle size distribution and local absorption peaks [26]. However, if the population of scatterers, or the absorption by dissolved matter, are not uniformly mixed over the sunlight penetration depth, this should show as larger inconsistencies in the retrieved values of $b\mathrm {_b(\lambda )}$. Step 6 in the proposed scheme therefore inspects the retrieved and filtered set of $b\mathrm {_b(\lambda )}$ for spectral consistency. The selection of suitable wavebands is expected to correspond to the molecular vibrations of the water molecule, because these are associated with (from shorter to longer wavelengths) increasingly efficient absorption properties. This property of the results distribution is used to further filter out results which are not consistent with the surrounding spectrum, by removing clusters with $\leq$ 3 results. In the examples given below, one-way analysis-of-variance (ANOVA, scipy v1.0.0) is used as a first-order test of consistency of $b\mathrm {_b(\lambda )}$ results between clusters of wavebands grouped by their nearest vibrational moment, while the Quartile Coefficient of Dispersion (QCD) is used to inspect how closely the full set of results is grouped. The QCD is calculated as half the interquartile range normalized to the midpoint of their range, i.e. ((Q3 – Q1)/2) / ((Q3 + Q1)/2) where Q1 and Q3 are the 25th and 75th percentiles, respectively.

2.2 In situ reflectance data

Hyperspectral reflectance and associated biogeochemical and optical sample analysis data were sourced from the community-owned Lake Bio-optical Measurements and Matchup Data for Remote Sensing (LIMNADES) database [27]. The datasets used here originate from campaigns in six lakes in the United Kingdom (Bassenthwaite Lake, Coniston Water, Derwent Water, Loch Leven, Windermere, Ullswater and Loch Lomond), three South African reservoirs (Loskop, Hartbeespoort, Theewaterskloof), Lake Geneva in Switzerland and the Lithuanian part of the Curonian Lagoon. All reflectance data were above-surface remote-sensing reflectance ($R\mathrm {_{rs}}$, units sr$^{-1}$) and ranged from the ultraviolet (350 nm) into the near infrared (800-1000 nm depending on instrument used) while satisfying the criteria for collection of above-water reflectance measurements set out in [28–30], obtaining radiance and irradiance signals of the water and sky either simultaneously or sequentially. Radiance sensors were pointed between 90 and 135$^{\circ }$ away from the solar azimuth and pointed 40-42$^{\circ }$ from the vertical at either the sky or the water surface. In the European lakes, two Satlantic HyperSAS spectroradiometers were positioned on the bow of a small vessel to simultaneously record water-leaving radiance and downwelling sky radiance, while downwelling irradiance was recorded with a cosine-corrected instrument. $R\mathrm {_{rs}}(\lambda )$ was calculated using the fingerprint method [31], which lets a scalar correction factor for skylight reflected on the water surface vary while optimizing the $R\mathrm {_{rs}}$ spectrum until atmospheric absorption features are minimized. The (ir)radiance spectra were interpolated to a common 3.3-nm interval prior to deriving reflectance. In all other cases, the skylight reflectance factor was assumed 0.0256, i.e. the Fresnel reflectance for a flat surface corresponding to a viewing angle 42$^{\circ }$ and refractive index of freshwater (1.333 relative to air), while instruments were hand-held, maintaining the optimal viewing geometry, and reported at a 1-nm spectral interval. For the South African reservoirs an Analytical Spectral Devices FieldSpec3 was used, sequentially recording water-leaving radiance, sky radiance and the radiance from a horizontally positioned 99 % Spectralon diffuse reflectance panel (calibrated to yield a signal equivalent to downwelling irradiance). In the Curonian Lagoon a WISP-3 handheld spectroradiometer was used to simultaneously obtain sky and water radiance and downwelling irradiance, spectrally binned to a common 1-nm grid. The WISP-3 internally averaged three sequential measurements on each channel before producing $R\mathrm {_{rs}}(\lambda )$. Further differences between the instruments, such as their field-of-view, spectral resolution (bandwidth), and post-processing of the measurements from each individual instrument are not considered relevant here.

An overview of the reflectance data from field measurements and those generated in optical simulations (described further below) is provided in Fig. 2. The data span a wide range in phytoplankton biomass (chlorophyll-a up to 350 mg m$^{-3}$) but do not include cases of very high turbidity due to mineral matter, such as may be found in and near rivers or shallow areas.

 

Fig. 2. Above-surface remote-sensing reflectance ($R\mathrm {_{rs}}(\lambda )$, sr $^{-1}$) spectra (A) obtained from in situ datasets and (B) generated in optical simulations, coloured by the log-transformed concentration of chlorophyll-a extracted from surface water samples.

Download Full Size | PDF

2.3 Validation and assessment

Optical and biogeochemical data were used here to demonstrate the waveband selection scheme over a range of optical conditions. In situ $b\mathrm {_b(\lambda )}$ are generally scarce and limited to a small set of wavebands, particularly from depth-profiles. Therefore, validation of $b\mathrm {_b(\lambda )}$ is done here in controlled optical simulations (described next). The impact of retrieving a band set of $b\mathrm {_b(\lambda )}$ is then assessed by comparing results from a semi-analytical red to near infra-red band ratio algorithm to retrieve the chlorophyll-a concentration [21]. In this algorithm, the term which retrieves $b\mathrm {_b}$ from a single near infra-red waveband (778 nm) was replaced with the median of retrieved $b\mathrm {_b(\lambda )}$ values from the proposed method. Retrieved chlorophyll-a results using both single-band and hyperspectral retrieval of $b\mathrm {_b}$ were then compared to extracted pigment concentrations from surface water samples. The variability in $b\mathrm {_b(\lambda )}$ was further propagated to chlorophyll-a estimates.

For the South African sites, chlorophyll-a was determined spectrophotometrically from extracts produced with boiling ethanol (95%) as described in [32]. Sampling and analytical procedures for pigments and absorption properties are described in further detail in [33]. For the European lakes sampled in the GloboLakes project, chlorophyll-a was similarly extracted in hot ethanol (90 %) according to ISO 10260:1992. For the Curonian Lagoon, pigments were extracted in 90 % acetone for 24 h at 4$^{\circ }$C and quantified spectrophotometrically according to [34], including sample acidification (1N HCl) to correct for phaeopigment.

EcoLight v5.2.0 (Sequoia scientific) was used to simulate $R\mathrm {_{rs}}(0^+,\lambda )$ over a wide concentration range (see below) of particles having phytoplankton absorption and scattering properties to evaluate whether uniformly mixed versus shallow and deep stratification conditions can be adequately handled by the inversion scheme. A water column of 25 m depth was simulated using only this scatterer and pure water. The scatterer was placed in discrete layers starting at either 0, 1, 3, 5, or 10 m depth, and extending downward over an interval ranging between 1 and 25 m depth, covering 20 variations in total (0-1 m, 0-2 m, 0-3 m, 0-5 m, 0-10 m, 0-15 m, 0-20 m, 0-25 m, 1-3 m, 1-15 m, 1-25 m, 3-5 m, 3-15 m, 3-25 m, 5-7 m, 5-15 m, 5-25 m, 10-12 m, 10-15 m, 10-25 m). Clear water was always included in the simulation below this layer of scatterers. The simulations were run at three chlorophyll-a concentrations (1, 5 and 10 mg m$^{-3}$) using phytoplankton absorption and (back)scattering properties recorded in the Baltic Sea during spring bloom [35]. The resulting reflectance (Fig. 2(B)) was inverted using the proposed procedure to yield the distribution of $b\mathrm {_b(\lambda )}$ in the selected wavebands. It should be noted that the optical properties of the scatterer used in the simulations is arbitrary. While chosen to reflect a realistic range of concentrations from relatively clear to more productive water bodies, the effect of a non-uniformly mixed scatterer on the distribution of retrieved $b\mathrm {_b(\lambda )}$ is expected to show for all types of scatterers.

3. Results

3.1 Inversion of in situ reflectance spectra

Results from the proposed inversion procedure are illustrated in further detail through a number of examples from optically varying water bodies (Figs. 3–6). The first three examples are from South African reservoirs, where a wide range of optical properties can be encountered and reflectance was measured in the 350-950 nm range. A further example from the Curonian Lagoon (Lithuania) uses reflectance in 400-800 nm range. These examples serve to demonstrate the retrieval capability across a range of optically complex water types which are not expected to be severely affected by non-uniformity of the water column.

 

Fig. 3. Stepwise retrieval of $b\mathrm {_b(\lambda )}$ from Loskop dam, 29 July 2011. (A) Normalized, smoothed reflectance and absorption by pure water, corresponding to processing Steps 1-2. The dashed line marks $a\mathrm {_w}$ = 0.1 m$^{-1}$. (B-C) Resulting $b\mathrm {_b}$ values obtained from Eq. (3), assuming dominant $a\mathrm {_w}$. Wavebands selected in processing Step 3 are overlaid in panel C as circular markers, grouped (marker colour) by their nearest molecular vibration wavelength of water (panel C legend). Statistically similar group means (Bonferroni-corrected t-test) are indicated by matching characters (A, B, C) at the top of panel C. Panel B provides a histogram and QCD of the retrieved $b\mathrm {_b(\lambda )}$ estimates. (D) Spectral slopes for reflectance and (inverted) $a\mathrm {_w}$ corresponding to processing Step 3A. Wavebands with slopes matching within $\pm$0.05 are marked. (E) Derivatives (Step 3B) of the slopes from panel B, with locations of matching derivatives within 0.025 nm$^{-1}$ marked.

Download Full Size | PDF

 

Fig. 4. Stepwise processing results from Hartbeespoort dam, 19 Oct 2010. Legends are as in Fig. 3 and red dashed line in panel A indicating that a red-edge was identified.

Download Full Size | PDF

 

Fig. 5. Stepwise processing results from Theewaterskloof dam, 26 Apr 2012. Legends are as in Figs. 3–4.

Download Full Size | PDF

 

Fig. 6. Stepwise processing results from the Curionian Lagoon, 3 Sep 2014. Legends are as in Figs. 3–5.

Download Full Size | PDF

The first example (Fig. 3) is from Loskop dam (25.453$^{\circ }$S, 29.288$^{\circ }$E) for an observation recorded on 29 July 2011 in relatively clear water (Secchi disk depth 2.35 m, chlorophyll-a 2.8 mg m$^{-3}$) with strong dissolved and non-phytoplankton particulate absorption determined spectrophotometrically from water samples (38 % and 43 % of absorption at 440 nm, respectively) and <20 % inorganic matter in a particulate dry weight of 2.63 g m$^{-3}$. The resulting reflectance spectrum has a distinct peak in the green region and no obviously discernible features corresponding to pigment absorption (no red edge was detected). When comparing normalized reflectance and $a\mathrm {_w}$ spectra (Fig. 3(A)), it is clear that significant drops in reflectance coincide with the stepwise increase of $a\mathrm {_w}$ around the H$_{2}$O vibrations at 606, 660, 739 and 836 nm. Locations along the spectrum where the reflectance and the inverse of $a\mathrm {_w}(\lambda )$ have similar local waveband ratios are seen to coincide with the same features in $a\mathrm {_w}(\lambda )$ (Fig. 3(D)) and notably include results at wavelengths as short as 600 nm. Further filtering by the derivative of the ratio spectra (Fig. 3(E)) removes wavebands associated with the weakest changes in reflectance (and higher noise), and leaves only three results in the region around 606 nm. The remaining wavebands (11 around 660 nm, 17 around 739 nm and 16 around 839 nm) are marked in Fig. 3(C), superimposed (for illustration purposes) on the spectrum of $b\mathrm {_b(\lambda )}$ that would be obtained when assuming all absorption is due to $a\mathrm {_w}(\lambda )$. A relatively high degree of consistency is observed between the values of $b\mathrm {_b}$ in the range 0.020 – 0.045 m$^{-1}$), in regions corresponding to the stepwise increase in absorption efficiency by water, although the 660 nm region includes a high outlier. ANOVA ($p$<0.0001) statistically separates the three clusters, despite an apparent close distribution of values in the 739 nm and 836 nm clusters. Bonferroni-corrected t-tests ($p$<0.017) for the three cluster pairs also reject close similarity. The QCD for this sample is 2.2 %.

The second example (Fig. 4) is from Hartbeespoort dam (25.7479$^{\circ }$S, 27.8635$^{\circ }$E), recorded on 19 October 2010 in eutrophic water where the optical properties were dominated by a high abundance of phytoplankton (Secchi disk depth 0.74 m, chlorophyll-a 3443.0 mg m$^{-3}$, suspended solids dry weight 500 g m$^{-3}$), in particular the cyanobacterium Microcystis aeruginosa. Under calm weather, the buoyancy of this species can overcome vertical mixing leading to accumulations at or near the surface. At the time of sampling, wind speed was 10 m s$^{-1}$ which would have introduced mixing. A red edge is clearly visible (Fig. 4(A)) and therefore only wavebands > 700 nm are considered. Wavebands suitable for inversion are identified only in the regions around 739 and 836 nm (n=34 and 20, respectively) and the two clusters belong to the same population if judged by t-test ($p$=0.20), with only few of the values in the 839 nm region clearly separated from the median at 0.46 m$^{-1}$. The QCD is at 0.8 % in this case.

The third example (Fig. 5) is from Theewaterskloof dam (34.029$^{\circ }$S, 19.209$^{\circ }$E), recorded on 26 April 2012 in moderately eutrophic but turbid conditions (Secchi disk depth 0.65 m, chlorophyll-a concentration 35.3 mg m$^{-3}$, suspended solids 17.13 g m$^{-3}$). Predominant cyanobacteria and diatom species included Anabaena ucrainica, Sphaerodinium fimbriatum, Asterionella formosa and Aulacoseria ambigua. Inorganic matter made up 45 % of particulate dry weight. The phytoplankton component is less pronounced and mineral particles contribute relatively more to light scattering compared to the previous example. Retrieved $b\mathrm {_b(\lambda )}$ was again from the 739 and 836 nm regions (n=32 and 19, respectively) and ranged 0.12 – 0.16 m$^{-1}$, with statistically significant differences (t-test, $p$<0.05) between these regions whereas the QCD was 1.1%. It is interesting to note that in all examples shown thus far, reflectance data spanned the wavelength range up to 950 nm but retrieved $b\mathrm {_b(\lambda )}$ was consistently lower in the region > 850 nm than at shorter wavelengths. In this and the previous example, results asssociated with the 836 nm region fall into two clusters, with those < 836 nm very similar to retrieval around 739 nm and those from > 850 nm distinctly lower.

The final example is from the Curionian Lagoon (Fig. 6). Sampling was carried out at Vente station (55.33568$^{\circ }$N, 21.18489$^{\circ }$E) using a WISP-3 (Water Insight, The Netherlands) to obtain reflectance in the 400-800 nm region which is narrower than previous examples. The sample shown was recorded on 3 September 2014, in similarly eutrophic, but less turbid, conditions compared to the previous example (Secchi disk depth 1.6 m, chlorophyll-a concentration 30.0 mg m$^{-3}$, suspended solids 9 g m$^{-3}$). A phycocyanin : chlorophyll-a ratio of 2.1 indicates that cyanobacteria were abundant. Inorganic matter was 36 % of particulate dry weight. $b\mathrm {_b(\lambda )}$ was, in this case, only retrieved from the 739 nm region (n=16), showing three sub-groups within a narrow distribution (Fig. 6(D), left panel) around 0.05 m$^{-1}$ and a low QCD of 1.2 %.

3.2 Validation with simulated reflectance

Optical simulations, with phytoplankton as the sole scatterers distributed over varying depth intervals, show that the range of retrieved $b\mathrm {_b(\lambda )}$ was narrowest for cases where the population of scatterers was evenly distributed from zero depth, or when a clear water layer of at least 3 m (for the lowest concentration) was overhead (Fig. 7). Shorter wavebands proved most sensitive to non-uniformity in the vertical distribution, departing further from the median. This tended towards underestimation when the scatterer layer was shallow (from 0 to several meters from the surface depending on concentration) and towards overestimation when the surface layer was uniform from 0 down to 5-25 meters in the case of the lowest phytoplankton concentration. The presence of these outliers around an otherwise relatively narrow distribution would suggests that the median of retrieved $b\mathrm {_b(\lambda )}$ provides a reasonable first order of estimate of $b\mathrm {_b(\lambda )}$ when mixing conditions are either uniform, or as long as the optical properties of the depth layer containing the scatterer are similar to the average of the observed layer.

 

Fig. 7. Retrieval of $b\mathrm {_b(\lambda )}$ for a range of depth distributions of a simulated phytoplankton population at three chlorophyll-a concentrations (top, middle and bottom panel). From left to right, scatterer distributions within the 25-m water column start at the surface over an increasing depth layer (0 – 1 m to 0 – 25 m) followed by sub-surface layers and deeper layers with clear water overhead. Boxplots (first and third quartiles) include the median (orange line). $b\mathrm {_b}$ coefficients retrieved at individual wavebands are overlaid as dots coloured according to the nearest H${_2}$O molecular vibration moment (legend in top-panel). The QCD is indicated on a scale of 0-1 at the top of each panel.

Download Full Size | PDF

When the phytoplankton layer was simulated below a layer of clear water of varying depth, but within the observed light penetration depth, results were more varied. This is marked by a broader interquartile range of the retrieved $b\mathrm {_b(\lambda )}$. In these cases, wavebands in the 600-700 nm range, where light penetration is presumed deepest, captured $b\mathrm {_b}$ by the scatterer, leading to $b\mathrm {_b(\lambda )}$ values comparable to shallower distributions, whereas longer wavebands captured weaker average backscattering by interacting less with the scatterer. The results for $b\mathrm {_b}$ around 970 nm are in some cases higher than expected given the efficient absorption by water. Whilst some of the retrieved values in this spectral region corresponded to clear water (low $b\mathrm {_b}$) as expected, others showed up in the mid-range of the distribution. This is most clearly seen in the lowest concentration range. At higher concentrations, fewer results were produced in this spectral region, which suggests relatively poor spectral matching of the reflectance and water absorption signatures. The QCD was in the order of 2-5 % when the scatterer extended from the surface. For subsurface (1-3 m clear water above the scatterer) and deeper distributions (10 m clear water above the scatterer), QCD exceeded 30 %. For deep layers (> 5m) at the lowest concentration, QCD ranged 15-19 %.

3.3 Impact on chlorophyll-a retrieval from the red to near-infrared spectrum

The influence of hyperspectral $b\mathrm {_b}$ retrieval on chlorophyll-a estimates was assessed using the semi-analytical algorithm of [21]. This algorithm was designed for use with narrow (in the order of 10 nm) satellite sensor wavebands, relating the ratio of infra-red (709 nm) over red (665 nm) bands to pigment and water absorption, and inverting the reflectance in the near infra-red (778 nm) to $b\mathrm {_b}$ under the assumption that non-water absorption at 778 nm is negligible. The algorithm is selected here because it has a wide validated concentration range of 10-200 mg m$^{-3}$ [20,36,37] and because it provides a straightforward opportunity to replace the $b\mathrm {_b}$ estimate with that retrieved from the hyperspectral approach. The in situ data set contained 215 stations for which both chlorophyll-a concentration and reflectance data were available. Given the range of applicability of the algorithm we only consider the 167 data points in the concentration range 5-1000 mg m$^{-3}$. Further removing stations where the single-band retrieval of $b\mathrm {_b}$(778) fails (yielding a negative result, n=73) and/or where scums were detected (n=24) resulted in 94 observation points with an associated chlorophyll-a concentration in the range 5.0 to 355 mg m$^{-3}$.

Substitution of the median $b\mathrm {_b(\lambda )}$ for $b\mathrm {_b}$(778) improved algorithm performance, as shown in Fig. 8(A). Chlorophyll-a retrieval using the original NIR-based single-band retrieval of $b\mathrm {_b}$ resulted in a coefficient of determination ($R^2$) of 0.77, root-mean square error (RMSE, calculated against regression fit) of 75.45 mg m$^{-3}$, bias 0.51 and mean absolute percentage error 74.2 % (MAPE, based on each data pair). The regression slope was 0.77 with an intercept of 11.67 mg m$^{-3}$. After substitution of $b\mathrm {_b}$(778) with the spectral median $b\mathrm {_b}$, results improved to $R^2$ = 0.82, RMSE 44.13 mg m$^{-3}$, bias 0.40 and MAPE = 62.74 %, regression slope 0.82 and intercept of 11.67 mg m$^{-3}$. Improvements are particularly relevant in eutrophic water bodies (Curonian Lagoon, Loskop dam), as shown in the per-lake regression statistics in Table 1, whereas performance in mesotrophic water bodies was comparable in most cases but worse performance was also observed (MAPE increased from 52.0 to 130.6 % in Theewaterskloof dam, whilst RMSE improved from 3.9 to 2.2 mg m$^{-3}$).

 

Fig. 8. (A) Chl-a retrieval adopted from [21] using the original method which resolves $b\mathrm {_b(778)}$ (red symbols and linear fit) and using the median $b\mathrm {_b}$ from hyperspectral inversion (blue symbols and linear fit). Regression fits (see main text) are based on data before log transformation and dashed line marks unity. (B) The mean absolute difference of Chl-a retrieval (vertical lines) around the in situ Chl-a sample concentration (markers as in panel A legend), plotted as a function of the mean (crosses) and standard deviation (horizontal lines) of all $b\mathrm {_b}$ values retrieved per sample, and their linear regression fit (before log transformation, drawn line). (C) Frequency of $b\mathrm {_b}$ results as a function of wavelength (all $b\mathrm {_b}$ wavebands, all samples). (D-F) Frequency plots of (D) $b\mathrm {_b}$ values from all bands, all samples; (E) spectral-median $b\mathrm {_b}$ per sample; (F) $b\mathrm {_b(778)}$ per sample.

Download Full Size | PDF

Table 1. Linear regression results for chlorophyll-a retrieval performance, per water body for which > 3 data points are returned, with single NIR-band $b\mathrm {_b}$ retrieval (’NIR’) or with the median of hyperspectrally derived $b\mathrm {_b}$ (’hyper’).

View Table

Uncertainty in the retrieval of chlorophyll-a can now be propagated from the variability in the retrieval of $b\mathrm {_b(\lambda )}$ from multiple wavebands. This is illustrated in Fig. 8(B) by showing both the variability in retrieved $b\mathrm {_b(\lambda )}$ as the standard deviation around the mean, and the associated spread in chlorophyll-a estimates, plotted around the in situ observed chlorophyll-a concentration. Chlorophyll-a estimates for a single reflectance sample are linearly proportional with $b\mathrm {_b(\lambda )}$ (not shown), and therefore a wider range of $b\mathrm {_b(\lambda )}$ is associated with a wider range of uncertainty in chlorophyll-a. For the dataset as a whole, a weak correlation is found between chlorophyll-a and mean $b\mathrm {_b(\lambda )}$ (linear regression slope = 0.16, intercept 40.10 mg m$^{-3}$, $R^2$ = 0.16, on data before log-transformation).

Figure 8(C-F) compare the retrieved distributions of $b\mathrm {_b(\lambda )}$ across the dataset. The first panel shows the frequency of wavebands contributing to the median $b\mathrm {_b(\lambda )}$ estimates used in this analysis, and panels D-F compare the resulting $b\mathrm {_b(\lambda )}$ including all bands, the spectral median, and the $b\mathrm {_b(778)}$ per sample, respectively. We note that observations from the Curionian Lagoon do not contain reflectance > 800 nm, which may explain a lower frequency in this range in Fig. 8(C). The number of results returned in the 600-700 nm range is lower because it is not used when a red edge shows in the reflectance spectrum. Distributions of the spectral median and the full $b\mathrm {_b(\lambda )}$ distribution are highly similar, with the majority of results < 0.1 $b\mathrm {_b(\lambda )}$ and remainder associated with the most eutrophic and turbid waterbodies (Loskop, Theewaterskloof dams and Curonian Lagoon). Comparing distributions of $b\mathrm {_b(778)}$ and the spectral median $b\mathrm {_b}$, however, suggests clearer separation of these two groups of waterbodies using the spectral median, than observed with $b\mathrm {_b(778)}$.

4. Discussion

Hyperspectral radiometry is now the established standard for above-water in situ water-leaving reflectance observations, ranging from handheld to shipborne and airborne platforms, and gaining ground in new spaceborne missions (e.g. PRISMA, EnMap, SHALOM, PACE, FLEX, HyspIRI and HYPXIM). These missions are complementary to multispectral ocean colour missions, adding enhanced spectral resolution at lower temporal frequency. Regardless of the scale of the application, there is a growing need to evolve reflectance spectroscopy from interpreting key reflectance wavebands to the exploitation of spectrally resolved measurements. A key difference between these multi- and hyperspectral approaches is that the latter allow inspection of consistency of information derived from a contiguous range of wavebands, such as demonstrated here in the retrieval of $b\mathrm {_b(\lambda )}$ allowing for initial spectral consistency checks in the orange to NIR part of the spectrum. At the same time, however, validation of hyperspectral $b\mathrm {_b}$ estimates are hampered by lack of suitable in situ instrumentation outside of key wavebands.

The only constant in the interpretation of water-leaving reflectance is the well-characterized absorption property of water, which is only altered by temperature and salinity (density). Salinity effects may be considered negligible in the spectral domain (> 600 nm) considered here. Temperature effects were not clearly observed here but may have influenced retrieval accuracy in the region around 750 nm where temperature effects on absorption by water are in the order of 0.5 % per degree, which is significant in this waveband but may be considered non-critical in the retrieval of $b\mathrm {_b(\lambda )}$ from a wide spectral range. A combination of simulated and field observation retrieval results, representing a range of relatively clear and optically highly complex waterbodies, confirm that $a\mathrm {_w}$ features are observable in reflectance spectra of natural waters from the yellow-orange domain where the absorption is in the order of 0.1 m$^{-1}$. For reference, a pigment absorption coefficient of 0.1 m$^{-1}$ at the red peak of chlorophyll-a would correspond to a concentration in the order of 5 mg m$^{-3}$ for a pigment-specific absorption coefficient in the order of 0.02 mg m$^{-2}$ [38].

Optical simulations show that the distribution of retrieved $b\mathrm {_b(\lambda )}$ widens when absorbing and scattering matter is not mixed homogeneously within the observed depth layer. This effect was most strongly observed when scatterers were concentrated in discrete layers with clear water overhead. In near-surface stratification conditions the divergence is expected to be narrower because the backscattered signal will predominantly originate from near the surface. The overall effect is wavelength-dependent, relating to the specific light penetration depth in each considered waveband. As the discrepancy widens, accurate inversion of water-leaving reflectance becomes unlikely, at least without employing further efforts to model vertical distribution effects.

It is highly desirable to define a test on the $b\mathrm {_b(\lambda )}$ value distribution which can consistently identify the presence of non-uniform mixing, but this proves challenging. Testing for modality (not shown) of the returned $b\mathrm {_b(\lambda )}$ distribution did not prove useful to distinguish uniformly mixed columns from heterogeneous mixing in either simulated or field data because the number of individual $b\mathrm {_b(\lambda )}$ results per sample is too small to generate continuous value distributions. Grouping the distribution of $b\mathrm {_b(\lambda )}$ retrievals by the vibration energy of water molecules is one way to discern differences between spectral regions where they are most likely to diverge, and we observe that the wavebands selected for retrieval indeed cluster around these regions. However, conventional tests to compare distributions (t-tests and ANOVA) did not yield a consistent approach. This may be caused by a significant range of $a\mathrm {_w}$ within these regions, and separating them further could be considered. Inspecting the results from individual wavebands, rather than in clusters, is of interest because a relationship between $a\mathrm {_w(\lambda )}$ and $b\mathrm {_b(\lambda )}$, across a spectral range where the absorption by water induces different light penetration depths, should not be expected when the water is homogeneously mixed. However, such a test is not straightforward because of the observed clustering, variable number of responses per region, and large data gaps in between (result not shown). Ultimately, the clustering of $b\mathrm {_b(\lambda )}$ per H$_2$O overtone does summarize which spectral domains contain significant water absorption signatures, as a proxy for the optical complexity of the sample, which could be used to determine further processing steps, including algorithm selection for biogeochemical properties. The sharp drop in $b\mathrm {_b(\lambda )}$ observed around 839 nm in two of the data processing examples from relatively turbid water conditions was common in the in situ data set and requires additional consideration. Overall, the in situ observations are not expected to include poor mixing conditions, prompting the question whether the wavelength range > 850 nm provides reliable estimates. Simulated conditions do not show this clear deviation, so signal saturation issues may need to be considered. Reflectance at longer wavebands, particularly in turbid conditions, will originate predominantly from shallow depths, and near-surface effects with lower concentrations of scatterers are not implausible (e.g. motile species moving away from the surface). On the other hand, bubbles or films of organic material are more likely to increase the backscattering signal in the same range. In clear but shallow water columns, the effect of bottom reflectance on the $b\mathrm {_b(\lambda )}$ would be worthwhile to investigate since it could potentially help to identify this condition.

In the search for a robust metric to determine (non-)uniform mixing, the quartile coefficient of dispersion included in the simulated and in situ data examples was one of the more robust indicators of spectral inconsistencies, suggesting that samples with QCD < 5 % were less likely to have inversion issues if the mode of the distribution is carried forward to the next processing step. 69 % of the in situ data fell within this range and no robust correlation between the QCD and chlorophyll-a retrieval accuracy was found (result not shown). The eutrophic (turbid) inland waterbodies included in the in situ dataset are generally more likely to show near-surface accumulations due to positive buoyant cyanobacteria, than deeper (e.g. density or temperature) stratification, and it is not yet clear whether such cases can be unambigously identified. In simulations, shallow surface layers were not associated with strongly elevated QCD, while deeper stratified layers had QCD > 15 %. There, the $b\mathrm {_b(\lambda )}$ is retrieved from multiple spectral regions (vibrations) with different light penetration depths (because water absorption dominates yet increases with each overtone). Thus, when the light penetration depth associated with one spectral region reaches through the stratification, the $b\mathrm {_b(\lambda )}$ distribution is seen to widen. Because the QCD is relatively insensitive to outliers it does not increase rapidly when only one of the spectral regions shows a different response, even if this response is consistent within that region. Some weighting on the number of results observed per spectral region should therefore be considered.

When inspecting the $b\mathrm {_b(\lambda )}$ distributions it is important to consider that, the larger the number of vibration regions considered, the wider the spectral region becomes from which $b\mathrm {_b(\lambda )}$ is drawn. Then, because $b\mathrm {_b(\lambda )}$ is not spectrally neutral but follows a slope determined by the particle size distribution, the $b\mathrm {_b(\lambda )}$ distribution should also be expected to widen when more regions are included. Any quality control thresholds, therefore, should take into account the nature of the spectral distribution of $b\mathrm {_b(\lambda )}$. Further simulations may prove useful, in order to control these experiments with known inherent optical properties, whilst these thresholds should also be tested against larger datasets including independently observed $b\mathrm {_b(\lambda )}$ and underwater depth profiles of inherent or apparent optical properties, to represent all measurement uncertainties. It is, however, noted that instrumentation to record $b\mathrm {_b(\lambda )}$ in situ has conventionally been multi-spectral (in the order of 3-9 wavebands), making direct validation of the concept presented here, challenging. It would be very useful to apply the proposed method to long data sets of high-frequency reflectance observations from moored stations where independent information on vertical distributions, e.g. from profiling measurements, is obtained.

Hyperspectral retrieval of $b\mathrm {_b(\lambda )}$ showed improvements in the retrieval of chlorophyll-a concentration with a semi-analytical algorithm, particularly in the more eutrophic and turbid range of the data set. Whilst direct validation of $b\mathrm {_b(\lambda )}$ is not feasible using (at best) multispectral in situ backscattering sensors, we observe that the median $b\mathrm {_b(\lambda )}$ better approximated $b\mathrm {_b}$ in the red and near-infrared bandset used by the algorithm compared to $b\mathrm {_b(778)}$ for the majority of samples considered. The hyperspectral retrieval procedure produced results in all cases where the reflectance spectrum was non-negative. The improved result may, of course, be an indication that retrieval of $b\mathrm {_b}$ from a single band is more likely prone to noise or measurement artefacts, which have a lower influence in the consensus solution drawn from the median $b\mathrm {_b(\lambda )}$. It follows that retrieval of (median) $b\mathrm {_b(\lambda )}$ from the whole NIR, where $a\mathrm {_w}$ is likely to dominate, will be similar without the need for a waveband selection scheme, and this is how full-spectral inversion methods already work. It should also be noted that the use of the 778 nm band in the original algorithm was guided by the availability of wavebands on medium resolution satellite sensors and from in situ datasets at the time of publication, rather than being considered optimal for all water conditions. The advantage of the selection scheme is, nevertheless, to draw information from wavebands outside this range whilst using consistent selection criteria. By spanning a wider spectral range (including regions nearer the band ratio used in the algorithm), extrapolation errors should be diminished. It is tempting to extrapolate the spectral slope of $b\mathrm {_b(\lambda )}$, although we note that it is unlikely that non-uniformity of mixing and the spectral slope can be unambiguously resolved together due to gaps between the selected wavelength ranges.

The distributions of per-sample $b\mathrm {_b(\lambda )}$ observed within the in situ dataset vary widely between waterbodies and across a wide phytoplankton biomass range. This suggests that there is a wealth of information contained in the selected $b\mathrm {_b(\lambda )}$ bandset, theoretically related to light penetration depth as well as potential observation artefacts such as suggested above (films, bubbles near the surface) as well as potential signal saturation or sensor inter-calibration issues resulting in noise. Stricter quality control of (above-surface) reflectance observations available in community data repositories, and harmonization of reflectance protocols used in the field, would help reduce some of the uncertainties related to potential signal saturation and sensor inter-calibration. Per definition, when relying on above-water spectroradiometry alone, these effects can not be fully isolated, whilst propagation of the uncertainty from $b\mathrm {_b(\lambda )}$ to substance concentration estimates should prove useful to determine application-specific quality control thresholds.

In conclusion, it is shown here that selection of a set of wavebands, where absorption by water changes rapidly with wavelength, can be effective to prioritise such a bandset for the initial inversion of reflectance, optimising retrieval procedures and providing additional analytical information on the spectral consistency of the backscattering coefficient. This helps address the observation uncertainty associated with the unknown extent of vertical mixing in remote sensing applications. Whether this spectral (in)consistency, shown here as proof-of-concept, can be interpreted clearly as uniform versus non-uniform mixing still requires further definition, through optical simulation with known depth profiles but ideally through field experiments with optical measurements resolved hyperspectrally over depth. Such experiments would help to further develop and validate this proof-of-concept. Limiting disturbance of vertical gradients during in situ observation from (large) research vessels is a particular challenge, and the use of experimental enclosures could therefore also be considered. The distribution of backscatttering retrieved from a single sample can be used to propagate uncertainties in subsequent algorithms, which in remote sensing imaging studies may theoretically be associated with the presence of fronts, near-surface blooms and potentially other types of stratification. Stepwise implementation of the inversion scheme allows testing for a range of conditions that prevent accurate interpretation of the spectrum over given spectral intervals, such as prevalence of absorption by other water constituents. These waveband intervals are then further ignored in the analysis. The presence of scums or floating vegetation invalidates the expected relation between reflectance and inherent optical properties entirely and should be determined (e.g. using a band ratio threshold as done here) to remove such samples. The method can be applied to reflectance from clear to turbid waters as shown with hyperspectral data from in situ observations from lakes, reservoirs and a lagoon, and to airborne or spaceborne hyperspectral instruments, following appropriate atmospheric correction.