Water-leaving radiance is subject to depth variability of the water constituents. The optical penetration depth is strongly dependent on the wavelength λ, which allows to retrieve a non-uniform vertical profile of an optically-active constituent CTSM(z) from remote-sensing reflectance Rrs(λ,Cz). We define the apparent particle concentration CTSM,app(λ) of a vertically homogeneous water column whose Rrs(λ,Cconst) matches Rrs(λ,Cz). Subsequently, we define a vertically-weighted averaged particle concentration CTSM,ave(λ), only dependent on CTSM(z), and retrieve CTSM(z) by minimizing the error between CTSM,app(λ) and CTSM,ave(λ) with genetic algorithms. We conclude that the retrieval is excellent if the sub-surface maximum lays close to the surface or the background concentration of CTSM(z) is low. Conversely, results worsen for opposite conditions, due to insufficient signal strength from superimposed sub-surface maxima.
© 2014 Optical Society of America
Inherent optical properties (IOPs) vary with depth in almost all natural waters. In open ocean water, phytoplankton finds a trade-off between stress by UV at the surface, depth-decreasing light availability and depth-increasing nutrient concentrations at the so-called deep chlorophyll maximum [1–4]. In coastal and inland waters, this feature exists as well for the same reasons , but light availability is in addition affected by other constituents like inorganic particles, which also have inhomogeneous vertical distributions [6, 7].
Remote-sensing reflectance (Rrs(λ)) is used to retrieve the IOPs and subsequently concentrations of the water’s optically active constituents, usually chlorophyll-a (Cchl-a), dissolved organic matter (quantified by its light absorption at a reference wavelength ag(λ0)) and total suspended matter (CTSM)  and references therein. Rrs(λ) represents a vertically-weighted average of water-leaving radiances across wavelength-dependent penetration depths. But in order to prevent under-determined numerical inversion, vertical homogeneity is usually assumed for the constituent retrieval algorithms. Dissolving this assumption poses a critical challenge for remote sensing, given the ecological significance of the stratification of different water constituents.
The problem of vertical ambiguity means that different vertical distributions of TSM (CTSM(z)) can result in similar surface signals. Contrariwise, each surface signal represents a rough average for different vertical stratifications . But light penetration depth varies strongly with wavelength, therefore the depth considered in this average varies across the spectral Rrs(λ,Cz). If this variation is properly parameterized, there is a chance to retrieve the vertical structure of the IOPs. Several model studies have showed that Rrs(λ) is sensitive to variations in the vertical structure of the IOPs, like CTSM(z) [6, 7] and the pigments chlorophyll-a  and phycocyanin .
In this article, we propose a new approach for the reconstruction of a vertical profile CTSM(z) from its corresponding Rrs(λ,Cz) by taking into account the different penetration of light at different wavelengths. Our procedure consists of three steps: (1) Knowing Rrs(λ,Cz), we define for each wavelength the apparent TSM concentration CTSM,app(λ), calculated as the vertically-constant CTSM,const that provides the corresponding reflectance Rrs(λ,Cconst) equal to Rrs(λ,Cz). (2) We define another function CTSM,ave(λ) as a vertically-weighted average of CTSM(z) so that CTSM,ave(λ) can be used as a proxy for CTSM,app(λ). Therefore, a mathematical link between Rrs(λ,Cz) and CTSM(z) is established. (3) We reconstruct the profile CTSM(z) by minimizing the difference between CTSM,ave(λ), calculated by applying the calibrated g(λ,z) from step (2) and CTSM,app(λ), calculated in step (1) by using genetic algorithms (GA) . The model presented in this paper retrieves one independent vertical profile and needs the knowledge of the specific inherent optical properties (SIOPs) of TSM. Table 1 shows all symbols used in this paper, in order of appearance.
2.1 Construction of the simulated data set
The starting point for our model study is the construction of a simulated data set (SDS) consisting of many pairs of CTSM(z), Rrs(λ,Cz). We arbitrarily consider all TSM profiles CTSM(z) as Gaussian curves:Table 2.This scheme leads to 3024 different profiles. The values were chosen for a case study from clear to moderately turbid waters and should be adapted to the range of values to be expected in each real application.
For each CTSM(z), Rrs(λ,Cz) is calculated using the software package Ecolight 5.1.2 , which solves the azimuthally-averaged radiative transfer equation in the water given the IOPs and boundary conditions.
Total absorption and scattering are decomposed as:17] between 400 nm and 700 nm and those by Smith and Baker  from 700 nm to 800 nm.
Absorption by phytoplankton is modeled as proportional to chlorophyll-a concentration and to a mass-specific spectrum:13].
Absorption by dissolved organic substances is modeled by a spectral exponential shape:16], pp. 26-28.
Wavelength band limits are set from 400 nm to 800 nm in steps of δλ = 4 nm (optical bandwidth also 4 nm). The output is averaged for each band and referenced to the central wavelength. All sources of inelastic scattering are considered in the simulation. The quantum yield of fluorescence is set to 0.02. The Raman scattering coefficient is set to 2.6 × 10−4 m−1 at the reference wavelength of 488 nm. The time is set to 25th May 2012, 12 h UTC and the coordinates to 47.52 °N, 9.59 °E (Lake Constance), which results in a sun zenith angle of 28 °. Normalized sky radiances are computed using the sky model of Harrison and Coombes , for an average cloud cover of 0.3. Diffuse and direct sky irradiances are computed using the RADTRANX model, included in the Ecolight software. The clear-sky irradiances are adjusted for cloudiness by the model of Kasten and Czeplak . Other parameters are: atmospheric pressure (101,253 Pa), air mass (type 3), relative humidity (80%), precipitable water (2.5 cm), wind speed (4 m s−1), visibility (15 km), ozone concentration (330.1 Dobson units) and aerosol optical thickness (0.261 at 550 nm). The index of refraction of water is calculated for an average water temperature of 14 °C and zero salinity. The water column is considered infinitely deep.
2.2 Construction of the apparent TSM concentrations by spectral matching to a look-up table
The apparent TSM concentration CTSM,app(λ) is defined as the concentration CTSM,const retrieved from Rrs(λ,Cz) by spectral matching with Rrs(λ,Cconst). Prior to this spectral matching, we construct a look-up table (LUT) for a range of CTSM,const from 0 to 20 mg l−1. We vary CTSM,const in steps of 0.05 mg l−1 (Fig. 1(a)), so that the LUT has 401 elements. The first 25 elements of the LUT for CTSM,const from 0 to 1.2 mg l−1 are plotted in Fig. 1.
Absorption and scattering of TSM are modelled as:Fig. 1(b)) are simulated using Ecolight.
CTSM,const (Fig. 1(a)) and Rrs(λ,Cconst) (Fig. 1(b)) have a biunivocal relationship for every wavelength. That is, increasing CTSM,const correspond to increasing Rrs(λ,Cconst) for every λ. Therefore, the correspondence Rrs(λ,Cz) ↔ CTSM,app(λ) is also biunivocal. Overlapped in Fig. 1(a) (bold line) is an example of an inhomogeneous CTSM(z) profile and its corresponding Rrs(λ,Cz), in Fig. 1(b). As shown, the range of CTSM,const from 0 to 1.2 mg l−1 corresponds to a set of Rrs(λ,Cconst) that cover all the dynamic range of Rrs(λ,Cz). An important feature is that it is not possible to fit Rrs(λ,Cz) to any of the Rrs(λ,Cz) in the LUT. Instead, Rrs(λ,Cz) crosses different Rrs(λ) spectra in the LUT. For instance, Rrs(406,Cz) = Rrs(406,Cconst) corresponds to line #11, with the corresponding concentration CTSM,const(line #11) = 0.5 mg l−1. Rrs(413,Cz) = Rrs(413,Cconst) (line #12) with CTSM,const(line #12) = 0.55 mg l−1, and so forth. Then, we define for Rrs(λ,Cz) the apparent CTSM,app(406) = 0.5 mg l−1 and CTSM,app(413) = 0.55 mg l−1. For the wavelengths between 406 nm and 413 nm, we interpolate linearly between 0.5 mg l−1 and 0.55 mg l−1. After following the procedure for all wavelengths, we obtain CTSM,app(λ) (Fig. 1(c)) corresponding to CTSM(z) from Fig. 1(a) (bold line) and Rrs(λ,Cz) (Fig. 1(b), bold line). CTSM,app(λ) tells how much TSM we see when we look from above at a given wavelength. As the light penetration depth changes strongly with λ, CTSM,app(λ) also changes in magnitude. For instance, at central wavelengths CTSM,app(λ) reaches the maximum because the clear water layer above is highly transparent and most of the water-leaving signal is from the turbid layer. At long wavelengths, absorption by water is so high that the water-leaving signal represents only the top centimeters of the water column, which leads to CTSM,app(λ) = 0, because the surface concentration is also 0 in this example. For other examples of CTSM(z), CTSM,app(λ) shows different spectral features, as next section will show.
2.3 Vertical averaging of a TSM profile
In this section, we present and calibrate a weighting function to vertically average a given profile CTSM(z) and obtain CTSM,ave(λ), with the condition that it matches CTSM,app(λ) as good as possible for every wavelength. Under the assumption of a biunivocal correspondence CTSM(z) ↔ CTSM,ave(λ), the aim is to retrieve CTSM(z) if CTSM,ave(λ) can be determined.
The only quantity that is vertically averaged is light. The averaging of any remote sensing quantity is an artificial procedure, so that the result of the averaging has to be defined a-priori (for instance, that it matches the result of a given remote sensing algorithm). In the context of oceanic chlorophyll, the weighting function g(λ,z) has been a matter of discussion in several papers. Gordon and Clark  proposed the round-trip attenuation:21] showed that Eq. (7) worked well in the vertical averaging of chlorophyll-a if absorption and scattering co-varied, and worked poorly if they did not. In our model study, there is strong covariance between absorption and scattering, so Eq. (7) is expected to be a good starting point for the search of a weighting function of TSM. Nevertheless, Eq. (7) gives preference to upper layers in the water column. In the presence of clear surface water and a turbid sub-surface layer, the latter has higher relative influence, as modeled by Eq. (7). For this reason, Zaneveld et al.  proposed the following formula instead, based on the derivative of Eq. (7):Eqs. (7) and (8):Eq. (10) adjust the contribution of the former to KTSM:Equation (13) should hold for every CTSM(z). However, as Eq. (11) is a heuristic model, we can only expect a reasonable match. Therefore, we proceed to minimize the error between both terms in Eq. (13).
A good calibration of Eq. (13) is expected to hold for all 3024 elements of the SDS. However, the computation time required for the calibration is proportional to the number of elements used, and would be too high if we tried to fit Eq. (13) for all the SDS. For this reason, we select 32 representative CTSM(z)i, with i = 1,…,32 from Table 2 for calibration. For each one of CTSM(z)i, we simulate Rrs(λ)i and the corresponding CTSM,app(λ)i. Every CTSM(z)i is determined by a given xi = (Cbg, Cmax, σ, zmax)i, thereby we make the dependence of CTSM,ave(λ,p,xi) on xi explicit as well. We minimize the RMS error between CTSM,app(λ)i and CTSM,ave(λ,p,xi), normalized to the λ-average of CTSM,app(λ)i and sum for all 32 profiles, as expressed in Eq. (13). A wavelength average is also performed, so that fp(p) only depends on p. The lower and upper bounds of the wavelength are, respectively, λmin = 400 nm and λmax = 800 nm and the wavelength range is Δλ = λmax - λmin = 400 nm.Eq. (12):Equation (15) shows that, for the case of absorption dominating over backscattering (our case, using the brown earth), KTSM ≈a roughly in Eq. (10). Note also that, the solution κsol = 4.51 makes our weighting function closer to Eq. (7) than to Eq. (8).
Using the calibrated values of Eq. (15), an excellent match is found between CTSM,app(λ) and CTSM,ave(λ) in this subset (results not shown). However the calibration in Eq. (14) is dependent on the used subset, therefore we must validate the calibration using an independent data set.
The calibrated parameters of Eq. (15) are now applied to a validation data set of another 32 TSM profiles chosen randomly from the SDS. The results (Fig. 2) are very satisfactory and no loss in fitness is observed. Hereafter, we use CTSM,ave(λ) as a proxy for CTSM,app(λ). The reader should note the different scale for each graph. For the cases CTSM,app(λ) is almost constant at a non-zero value, a fluctuating pattern can be appreciated. This error is around five times less than the precision of the LUT, e.g., <0.01 mg l−1 in the worst case we found, which is a negligible difference in real applications. A study in detail of Ecolight’s numerical output showed that the fluctuations come from the numerical precision of Ecolight, whose output reflectances are given to the smallest order of 10−6 sr−1. Another option would be to construct CTSM,app(λ) by taking the nearest element of the LUT instead of linearly interpolating. That would completely eliminate this effect, but would give CTSM,app(λ) a step-wise shape with a precision of 0.05 mg l−1.
2.4 Retrieval of an unknown TSM profile
In this section, we reconstruct CTSM(z) from a given Rrs(λ,Cz). Without loss of applicability for further studies, we suppose here CTSM(z) having an unknown Gaussian shape according to Eq. (1), thereby depending on the four parameters x = (Cbg, Cmax, σ, zmax). Ideally, the aim is to find x so that:Eq. (16), we make the dependence of CTSM,ave(λ,x) explicit on x, and omit the dependence on p, since p = psol has been fixed in Eq. (15).
Equation (16) does not have an analytical solution on x, thereby needing a numerical approach. For this sake, we use again GA. Cbg is searched between the bounds 0 and 5 mg l−1, Cmax between 0 and 10 mg l−1, σ between 0 and 3 m and zmax between 0 and 10 m. The number of generations is set to 25 and the population is set to 30. A goal function fx(x) to be minimized is also defined:Figure 3 shows how fx(x) (Eq. (17)) is constructed. Based on the established LUT, we estimate the apparent TSM concentration CTSM,app(λ) from a measured Rrs(λ,Cz). Independently, a guess of the Gaussian set x is used to construct a TSM profile, which is vertically averaged by Eq. (11) to obtain CTSM,ave(λ,x). Finally, the goal function is constructed as the Euclidean distance between CTSM,ave(λ,x) and CTSM,app(λ) along the entire λ range. No normalizations are applied in Eq. (17) because zero or near-zero values of CTSM,ave(λ,x) are to be expected for some cases of the SDS.
3. Results and discussion
We perform individual retrievals for each case of SDS. To quantify the performance of our model, we define an error function between the results of our model and the original CTSM(z) profiles (which, in a realistic situation, are not observed) for each case. In our study, all profiles are assumed to be Gaussian, so we define the following non-dimensional distance between the real values and the retrieved (subscript ret):Eq. (18) increases from zero as the differences between the Gaussian parameters increase. Note that, in a different situation, where CTSM(z) were not Gaussian, the distance d could be defined differently. An option would be a norm between real and retrieved profiles. Figure 4 shows selected original (blue) and retrieved (red) profiles, with the corresponding distance d for each retrieval. By visual inspection of several retrieved profiles, we set the threshold of a good retrieval at d = 0.4. From the plots, the retrievals look excellent if the maximum lays right at the surface (upper panels). However, if the maximum is deep in the water, other factors seem to influence the quality of the retrieval (lower panels).
We study in Fig. 5 which of the Gaussian parameters influence the goodness of the profile retrieval, and what are the interdependencies among them in the subset of the well-retrieved profiles (d < 0.4). Each graph shows, for each value of one of the Gaussian parameters (horizontal axis), the arithmetic mean of the other three parameters in the subset of the SDS that holds the condition d < 0.4. Figure 5(a) illustrates that, although zmax has a mean value for all the SDS of 5 m (cyan, dashed), the retrieval is good (cyan, continuous) only down to a depth of zmax ≈3 m in the best case, which is Cbg = 0 mg l−1. However, as Cbg increases up to 5 mg l−1, zmax decreases down to ~1 m. The explanation is the following: as Cbg increases, light penetration decreases, so that the TSM maximum has to be placed at upper layers so that d can remain lower than 0.4. In contrast, the sensitivity of Cmax to Cbg is weak. The same pattern is observed for σ. In Fig. 5(b), Cmax increases in the horizontal axis. The blue line (Cbg) reveals the same effect as in Fig. 5(a): as Cmax increases, Cbg must decrease, so that d can remain lower than 0.4. The parameter σ remains almost constant and zmax shows not a clear trend. The thickness σ of all profiles increases in the horizontal axis of Fig. 5(c) and variations in the average values of Cbg, Cmax and zmax, are plotted in the vertical axis. As all lines in this graph are fairly horizontal, σ is not a parameter that conditions the success of the retrieval. Last, the depth of the sub-surface maximum zmax is the independent parameter in Fig. 5(d). The strong decrease of Cbg (blue) indicates the explained effect: the deeper the profile is, the clearer the background water has to be so that light reaches the depth where TSM is concentrated.
In summary, Fig. 5 conveys the following conclusions: (1) in average, zmax < 3 m is needed in order to achieve good retrievals in the most favorable cases of clear background water. (2) If condition (1) applies, zmax and Cbg are negatively correlated: as Cbg increases from 0 mg l−1 to 5 mg l−1, zmax decreases on average from 3 m to 1 m. Conversely, as zmax deepens from [0,1] m to [6,10] m depth, Cbg decreases on average from ~2 to ~0 mg l−1. (3) The quality of the retrieval shows no sensitivity to Cmax and σ.
The reader must note that these numerical results are dependent on the particular IOPs of each case study. We focused on a case study of moderately turbid waters and a certain background chlorophyll-a and CDOM. For the case of cleaner waters, good performance down to greater depths is to be expected. On the contrary, for the case of very turbid surface waters, no vertical information at all can be extracted.
The algorithm’s performance using different numbers of bands is assessed in Fig. 6.We perform the whole procedure explained in this article for different numbers of bands (keeping the optical bandwidth at 4 nm). In black, we show results using all bands available (δλ = 4 nm), whereas the other colors correspond to the three wavelength steps δλ = 12 nm (brown), δλ = 36 nm (red) and δλ = 132 nm (green). Results deteriorate progressively as we use fewer bands, because the information in the shape features of the spectra is lost. However, the effect is not extreme. The spectra evaluated in our algorithm are CTSM,ave(λ) and CTSM,app(λ) (see Fig. 2), whose spectral features are basically the same as of Rrs(λ,Cz), and ultimately, of the IOPs. Our results suggest that for the retrieval of variations in CTSM, a 4 band multispectral sensor has almost the same potential as medium resolution ocean color sensors, provided an equally good radiometric resolution and spectral band characterization. This finding is not transferrable to Cchl-a retrieval, where narrow bands at appropriate wavelengths are crucial.
Although light is integrated on its vertical path through the water column, we showed that the vertical structure of a TSM profile leaves signatures in the remote sensing reflectance. Application of this model to our simulated data set showed that the most favorable case for the success of the retrieval is when the TSM profiles have high differences between maximum and background, and when the peak concentrations are close to the surface. The numerical results have to be interpreted for our specific simulated data set, consisting of low to moderate turbid waters with a background chlorophyll-a and CDOM. We do not discard room for further improvements in the method which could increase the sensitivity to deep layers. Besides this, there is always the physical limitation: That is, if the effect of a deeper layer on the water leaving radiance is tiny, that layer is not likely to be detected.
The retrieval of profiles by optical remote sensing is an ambitious and challenging task. Hence limited results are expected if the method is used as standalone. However, we suggest the integration of this model into a scheme of data assimilation of hydrodynamic simulations. Unlike other schemes that assimilate remote sensing products like TSM, this scheme would assimilate the reflectance, which contains the integrated information of the water inherent optical properties. Together with some ancillary data (e.g. the SIOPs) and other source of vertical water quality parameters distribution (e.g. field sampling or ecological modeling) this may lead to reasonably accurate retrieval of the vertical profiles.
A further step would be to address the more difficult problem of the retrieval of several water constituents. However, the success in this approach depends on an accurate knowledge of the IOPs after intensive field and laboratory work, which is also the requirement of physically-based remote sensing algorithms.
The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement n° 263287 (FRESHMON project). We thank Karin Schenk and Thomas Heege of EOMAP for the management of the project and two anonymous reviewers for their constructive review reports. The first author wants also to thank his new colleagues at ISAC Rome for their interest in this investigation. Their comments and suggestions on both the theoretical development and the practical applicability of this model were incorporated to a revised version of this article.
References and links
1. K. Fennel and E. Boss, “Subsurface maxima of phytoplankton and chlorophyll: Steady-state solutions from a simple model,” Limnol. Oceanogr. 48(4), 1521–1534 (2003). [CrossRef]
2. A. B. Ryabov, L. Rudolf, and B. Blasius, “Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer,” J. Theor. Biol. 263(1), 120–133 (2010). [CrossRef] [PubMed]
4. M. R. Clegg, U. Gaedke, B. Boehrer, and E. Spijkerman, “Complementary ecophysiological strategies combine to facilitate survival in the hostile conditions of a deep chlorophyll maximum,” Oecologia 169(3), 609–622 (2012). [CrossRef] [PubMed]
5. D. Odermatt, F. Pomati, J. Pitarch, J. Carpenter, M. Kawka, M. Schaepman, and A. Wüest, “MERIS observations of phytoplankton blooms in a stratified eutrophic lake,” Remote Sens. Environ. 126, 232–239 (2012). [CrossRef]
6. P. Forget, P. Broche, and J.-J. Naudin, “Reflectance sensitivity to solid suspended sediment stratification in coastal water and inversion: A case study,” Remote Sens. Environ. 77(1), 92–103 (2001). [CrossRef]
7. Q. Yang, D. Stramski, and M.-X. He, “Modeling the effects of near-surface plumes of suspended particulate matter on remote-sensing reflectance of coastal waters,” Appl. Opt. 52(3), 359–374 (2013). [CrossRef] [PubMed]
8. D. Odermatt, A. Gitelson, V. E. Brando, and M. Schaepman, “Review of constituent retrieval in optically deep and complex waters from satellite imagery,” Remote Sens. Environ. 118, 116–126 (2012). [CrossRef]
10. M. Stramska and D. Stramski, “Effects of a nonuniform vertical profile of chlorophyll concentration on remote-sensing reflectance of the ocean,” Appl. Opt. 44(9), 1735–1747 (2005). [CrossRef] [PubMed]
11. T. Kutser, L. Metsamaa, and A. G. Dekker, “Influence of the vertical distribution of cyanobacteria in the water column on the remote sensing signal,” Estuar. Coast. Shelf Sci. 78(4), 649–654 (2008). [CrossRef]
12. L. Davis, Handbook of genetic algorithms (Van Nostrand Reinhold, New York, 1991).
13. P. Gege, “Characterization of the phytoplankton in Lake Constance for classification by remote sensing (with 6 figures and 2 tables),” in Lake Constance, Characterization of an Ecosystem in Transition, E. Baeuerle and U. Gaedke, eds. (E. Schweizerbart'sche Verlagsbuchhandlung (Nägele und Obermiller), 1998), pp. 179–194.
16. C. D. Mobley and L. K. Sundman, Hydrolight 5 Technical Documentation (Sequoia Scientific, Inc., 2008), http://www.hydrolight.info.
19. A. W. Harrison and C. A. Coombes, “An opaque cloud cover model of sky short wavelength radiance,” Sol. Energy 41(4), 387–392 (1988). [CrossRef]
20. F. Kasten and G. Czeplak, “Solar and terrestrial radiation dependent on the amount and type of cloud,” Sol. Energy 24(2), 177–189 (1980). [CrossRef]