Retrieval of vertical particle concentration profiles by optical remote sensing : a model study

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 OCIS codes: (010.4450) Oceanic optics; (010.7340) Water; (280.0280) Remote sensing and sensors. 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). 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). 3. J. P. Mellard, K. Yoshiyama, E. Litchman, and C. A. Klausmeier, “The vertical distribution of phytoplankton in stratified water columns,” J. Theor. Biol. 269(1), 16–30 (2011). 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). 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). 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). 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). 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). 9. H. R. Gordon, “Remote sensing of optical properties in continuously stratified waters,” Appl. Opt. 17(12), 1893– 1897 (1978). 10. M. Stramska and D. Stramski, “Effects of a nonuniform vertical profile of chlorophyll concentration on remotesensing reflectance of the ocean,” Appl. Opt. 44(9), 1735–1747 (2005). 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). #207765 $15.00 USD Received 6 Mar 2014; revised 9 Apr 2014; accepted 11 Apr 2014; published 18 Apr 2014 (C) 2014 OSA 5 May 2014 | Vol. 22, No. S3 | DOI:10.1364/OE.22.00A947 | OPTICS EXPRESS A947 12. L. Davis, Handbook of Genetic Algorithms (Van Nostrand Reinhold, 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. 14. H. R. Gordon and D. K. Clark, “Remote sensing optical properties of a stratified ocean: an improved interpretation,” Appl. Opt. 19(20), 3428–3430 (1980). 15. J. R. Zaneveld, A. Barnard, and E. Boss, “Theoretical derivation of the depth average of remotely sensed optical parameters,” Opt. Express 13(22), 9052–9061 (2005). 16. C. D. Mobley and L. K. Sundman, Hydrolight 5 Technical Documentation (Sequoia Scientific, Inc., 2008), http://www.hydrolight.info. 17. R. M. Pope and E. S. Fry, “Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,” Appl. Opt. 36(33), 8710–8723 (1997). 18. R. C. Smith and K. S. Baker, “Optical properties of the clearest natural waters (200-800 nm),” Appl. Opt. 20(2), 177–184 (1981). 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). 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). 21. H. R. Gordon, “Diffuse reflectance of the ocean: influence of nonuniform phytoplankton pigment profile,” Appl. Opt. 31(12), 2116–2129 (1992).


Introduction
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][2][3][4].In coastal and inland waters, this feature exists as well for the same reasons [5], but light availability is in addition affected by other constituents like inorganic particles, which also have inhomogeneous vertical distributions [6,7].
Remote-sensing reflectance (R rs (λ)) is used to retrieve the IOPs and subsequently concentrations of the water's optically active constituents, usually chlorophyll-a (C chl-a ), dissolved organic matter (quantified by its light absorption at a reference wavelength a g (λ 0 )) and total suspended matter (C TSM ) [8] and references therein.R rs (λ) represents a verticallyweighted 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 (C TSM (z)) can result in similar surface signals.Contrariwise, each surface signal represents a rough average for different vertical stratifications [9].But light penetration depth varies strongly with wavelength, therefore the depth considered in this average varies across the spectral R rs (λ,C z ).If this variation is properly parameterized, there is a chance to retrieve the vertical structure of the IOPs.Several model studies have showed that R rs (λ) is sensitive to variations in the vertical structure of the IOPs, like C TSM (z) [6,7] and the pigments chlorophyll-a [10] and phycocyanin [11].
In this article, we propose a new approach for the reconstruction of a vertical profile C TSM (z) from its corresponding R rs (λ,C z ) by taking into account the different penetration of light at different wavelengths.Our procedure consists of three steps: (1) Knowing R rs (λ,C z ), we define for each wavelength the apparent TSM concentration C TSM,app (λ), calculated as the vertically-constant C TSM,const that provides the corresponding reflectance R rs (λ,C const ) equal to R rs (λ,C z ).(2) We define another function C TSM,ave (λ) as a vertically-weighted average of C TSM (z) so that C TSM,ave (λ) can be used as a proxy for C TSM,app (λ).Therefore, a mathematical link between R rs (λ,C z ) and C TSM (z) is established.(3) We reconstruct the profile C TSM (z) by minimizing the difference between C TSM,ave (λ), calculated by applying the calibrated g(λ,z) from step (2) and C TSM,app (λ), calculated in step (1) by using genetic algorithms (GA) [12].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.K TSM (λ,z) Attenuation coefficient defined in Eq. ( 10) to calculate the vertical weighting function in Eq. ( 9) Inverse of the exponent of K TSM (λ,z) in Eq. ( 9), which adjusts the vertical average from the exponential model of Gordon and Clark [14] its derivative, as proposed by Zaneveld et al. [15] non-dim α Coefficient, to adjust the relative weight of a(λ) in Eq. (10)

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 C TSM (z), R rs (λ,C z ).We arbitrarily consider all TSM profiles C TSM (z) as Gaussian curves: ( ) The Gaussian profile depends on four independent parameters, which we vary in discrete steps, according to 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 C TSM (z), R rs (λ,C z ) is calculated using the software package Ecolight 5.1.2[16], which solves the azimuthally-averaged radiative transfer equation in the water given the IOPs and boundary conditions.
Total absorption and scattering are decomposed as: Absorption a w (λ) and scattering b w (λ) of pure water are parameterized using the measurements of Pope and Fry [17] between 400 nm and 700 nm and those by Smith and Baker [18] from 700 nm to 800 nm.Absorption by phytoplankton is modeled as proportional to chlorophyll-a concentration and to a mass-specific spectrum: As this model study focuses on the retrieval of TSM, we set arbitrarily C chl-a = 1 µg l −1 and a* ph (λ) being the average of the specific spectra of the five dominant algal species in Lake Constance, as proposed by Gege [13].
Absorption by dissolved organic substances is modeled by a spectral exponential shape: We arbitrarily fix the reference absorption a g (440) = 0.2 m −1 and the decay constant S = 0.014 nm −1 .Spectral absorption by non-algal particles and scattering by TSM are modeled as proportional to C TSM (z) and the corresponding mass-specific spectra: The mass-specific spectra a* NAP (λ) and b* TSM (λ) are from the brown earth example provided with Ecolight [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 [19], 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 [20].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.

Construction of the apparent TSM concentrations by spectral matching to a look-up table
The apparent TSM concentration C TSM,app (λ) is defined as the concentration C TSM,const retrieved from R rs (λ,C z ) by spectral matching with R rs (λ,C const ).Prior to this spectral matching, we construct a look-up table (LUT) for a range of C TSM,const from 0 to 20 mg l −1 .We vary C TSM,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 C TSM,const from 0 to 1.2 mg l −1 are plotted in Fig. 1.
Absorption and scattering of TSM are modelled as: The mass-specific spectra a* NAP (λ) and b* TSM (λ) and the rest of IOPs and boundary conditions are set equal as in the SDS.The corresponding R rs (λ,C const ) (Fig. 1(b)) are simulated using Ecolight.C TSM,const (Fig. 1(a)) and R rs (λ,C const ) (Fig. 1(b)) have a biunivocal relationship for every wavelength.That is, increasing C TSM,const correspond to increasing R rs (λ,C const ) for every λ.Therefore, the correspondence R rs (λ,C z ) ↔ C TSM,app (λ) is also biunivocal.Overlapped in Fig. 1(a) (bold line) is an example of an inhomogeneous C TSM (z) profile and its corresponding R rs (λ,C z ), in Fig. 1(b).As shown, the range of C TSM,const from 0 to 1.2 mg l −1 corresponds to a set of R rs (λ,C const ) that cover all the dynamic range of R rs (λ,C z ).An important feature is that it is not possible to fit R rs (λ,C z ) to any of the R rs (λ,C z ) in the LUT.Instead, R rs (λ,C z ) crosses different R rs (λ) spectra in the LUT.For instance, R rs (406,C z ) = R rs (406,C const ) corresponds to line #11, with the corresponding concentration C TSM,const (line #11) = 0.5 mg l −1 .R rs (413,C z ) = R rs (413,C const ) (line #12) with C TSM,const (line #12) = 0.55 mg l −1 , and so forth.Then, we define for R rs (λ,C z ) the apparent C TSM,app (406) = 0.5 mg l −1 and C TSM,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 C TSM,app (λ) (Fig. 1(c)) corresponding to C TSM (z) from Fig. 1(a) (bold line) and R rs (λ,C z ) (Fig. 1(b), bold line).C TSM,app (λ) tells how much TSM we see when we look from above at a given wavelength.As the light penetration depth changes strongly with λ, C TSM,app (λ) also changes in magnitude.
For instance, at central wavelengths C TSM,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 C TSM,app (λ) = 0, because the surface concentration is also 0 in this example.For other examples of C TSM (z), C TSM,app (λ) shows different spectral features, as next section will show.

Vertical averaging of a TSM profile
In this section, we present and calibrate a weighting function to vertically average a given profile C TSM (z) and obtain C TSM,ave (λ), with the condition that it matches C TSM,app (λ) as good as possible for every wavelength.Under the assumption of a biunivocal correspondence C TSM (z) ↔ C TSM,ave (λ), the aim is to retrieve C TSM (z) if C TSM,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 [14] proposed the round-trip attenuation: Gordon [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. [15] proposed the following formula instead, based on the derivative of Eq. ( 7): We propose the following formula for the vertical averaging of TSM.The constant κ adjusts the weight of sub-surface maxima relative to the upper layers and adds a degree of freedom between the approaches of Eqs. ( 7) and ( 8): ( ) ( ) We parameterize K TSM as a function of absorption and backscattering.The constants α and β in Eq. ( 10) adjust the contribution of the former to K TSM : We define the vertical-weighted average of TSM as: We define the vector of unknowns: ( ) And impose the condition: We have made the dependence C TSM,ave (λ,p) on the vector of unknowns p explicit here.Equation ( 13) should hold for every C TSM (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 C TSM (z) i , with i = 1,…,32 from Table 2 for calibration.For each one of C TSM (z) i , we simulate R rs (λ) i and the corresponding C TSM,app (λ) i .Every C TSM (z) i is determined by a given x i = (C bg , C max , σ, z max ) i , thereby we make the dependence of C TSM,ave (λ,p,x i ) on x i explicit as well.We minimize the RMS error between C TSM,app (λ) i and C TSM,ave (λ,p,x i ), normalized to the λ-average of C TSM,app (λ) i and sum for all 32 profiles, as expressed in Eq. (13).A wavelength average is also performed, so that f p (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.
To solve for p, we apply a minimization routine using GA, included in the Global Optimization Toolbox of MATLAB (The Mathworks, Inc.).The lower bounds of κ, α, and β for parameter search are set to 0.2, 0.2 and 1, respectively, and the upper bounds, to 10, 10 and 10.The number of generations was set to 100 with a population of 50 individuals.The calibration procedure leads to the optimal choices of the parameters of Eq. ( 12): ( ) Equation shows that, for the case of absorption dominating over backscattering (our case, using the brown earth), K TSM ≈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 C TSM,app (λ) and C TSM,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.Fig. 2. C TSM,app (λ) (blue) and C TSM,ave (λ,p sol ) (red) calculated for 32 independent TSM profiles randomly chosen from the 3024 simulations of the SDS of Table 2.
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 C TSM,ave (λ) as a proxy for C TSM,app (λ).The reader should note the different scale for each graph.For the cases C TSM,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 C TSM,app (λ) by taking the nearest element of the LUT instead of linearly interpolating.That would completely eliminate this effect, but would give C TSM,app (λ) a step-wise shape with a precision of 0.05 mg l −1 .

Retrieval of an unknown TSM profile
In this section, we reconstruct C TSM (z) from a given R rs (λ,C z ).Without loss of applicability for further studies, we suppose here C TSM (z) having an unknown Gaussian shape according to Eq. ( 1), thereby depending on the four parameters x = (C bg , C max , σ, z max ).Ideally, the aim is to find x so that: Note that, in Eq. ( 16), we make the dependence of C TSM,ave (λ,x) explicit on x, and omit the dependence on p, since p = p sol 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.C bg is searched between the bounds 0 and 5 mg l −1 , C max between 0 and 10 mg l −1 , σ between 0 and 3 m and z max between 0 and 10 m.The number of generations is set to 25 and the population is set to 30.A goal function f x (x) to be minimized is also defined: Figure 3 shows how f x (x) (Eq.( 17)) is constructed.Based on the established LUT, we estimate the apparent TSM concentration C TSM,app (λ) from a measured R rs (λ,C z ).Independently, a guess of the Gaussian set x is used to construct a TSM profile, which is vertically averaged by Eq. ( 11) to obtain C TSM,ave (λ,x).Finally, the goal function is constructed as the Euclidean distance between C TSM,ave (λ,x) and C TSM,app (λ) along the entire λ range.No normalizations are applied in Eq. ( 17) because zero or near-zero values of C TSM,ave (λ,x) are to be expected for some cases of the SDS.
Vert. aver.Eq. ( 11) C TSM,ave (λ,x) x RMS error, Eq. ( 17) Fig. 3. Construction of the goal function f x (x) to retrieve the shape x of the unknown profile.
From R rs (λ,C z ), C TSM,app (λ) is obtained by LUT search.In parallel, a guess of solution x is used to build a vertical TSM profile, which is vertically averaged by Eq. ( 11) to obtain C TSM,ave (λ,x).Finally, the goal function is constructed as the wavelength average of the difference between C TSM,ave (λ,x) and C TSM,app (λ).

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 C TSM (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): The non-dimensional distance d of Eq. ( 18) increases from zero as the differences between the Gaussian parameters increase.Note that, in a different situation, where C TSM (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).Fig. 4. Selected original (blue) and retrieved (red) profiles, with the corresponding distance d (Eq.( 18)) for each retrieval.Axes are equal for all plots (bottom left) and are omitted for readability.Red: σ.Cyan: z max .In the same colors, the mean values of each parameter for all the SDS are plotted (dashed lines).Note that, in graph (d), all z max values are grouped in intervals.
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 z max 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 z max ≈3 m in the best case, which is C bg = 0 mg l −1 .However, as C bg increases up to 5 mg l −1 , z max decreases down to ~1 m.The explanation is the following: as C bg 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 C max to C bg is weak.The same pattern is observed for σ.In Fig. 5(b), C max increases in the horizontal axis.The blue line (C bg ) reveals the same effect as in Fig. 5(a): as C max increases, C bg must decrease, so that d can remain lower than 0.4.The parameter σ remains almost constant and z max 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 C bg , C max and z max , 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 z max is the independent parameter in Fig. 5(d).The strong decrease of C bg (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, z max < 3 m is needed in order to achieve good retrievals in the most favorable cases of clear background water.(2) If condition (1) applies, z max and C bg are negatively correlated: as C bg increases from 0 mg l −1 to 5 mg l −1 , z max decreases on average from 3 m to 1 m.Conversely, as z max deepens from [0,1] m to [6,10] m depth, C bg decreases on average from ~2 to ~0 mg l −1 .(3) The quality of the retrieval shows no sensitivity to C max 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 C TSM,ave (λ) and C TSM,app (λ) (see Fig. 2), whose spectral features are basically the same as of R rs (λ,C z ), and ultimately, of the IOPs.Our results suggest that for the retrieval of variations in C TSM , 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 C chl-a retrieval, where narrow bands at appropriate wavelengths are crucial.

Conclusion
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.

C− 1 C− 1 a 1 a 1 S 1 δλ 1 b
TSM,app (λ) Apparent TSM associated to R rs (λ,C z ), defined for each wavelength as the corresponding C TSM,const that causes R rs (λ,C z ) = R rs (λ,C const ) mg l TSM,ave (λ) Vertically-averaged TSM profile, calculated with Eq.(11).Dependences of C TSM,ave (λ) on x and p are made explicit only when needed for clarity reasons mg l * ph (λ) Mass-specific phytoplankton absorption spectrum, calculated by Gege[13] as the average phytoplankton absorption spectrum of the five most abundant algal species in Lake Constance m 2 mg −1 a g (440) CDOM absorption spectrum at the reference λ = 440 nm m −Spectral decay of the exponential model of a g (λ) nm −1 a* NAP (λ) Mass-specific non-algal absorption spectrum, taken from the brown earth example included in Hydrolight / Ecolight m 2 g −1 b* TSM (λ) Mass-specific TSM scattering spectrum, taken from the brown earth example included in Hydrolight / Ecolight m 2 g −Wavelength step (discretization) nm g(λ,z) Weighting function for depth averaging of a remotely sensed optical parameter non-dim K d (λ,z) Diffuse attenuation coefficient of the downwelling irradiance m −

Fig. 5 .
Fig.5.For every value of the Gaussian parameters C bg , C max , σ and z max (horizontal axis of each graph) the arithmetic mean of the other three parameters, belonging to the fraction of wellretrieved profiles of the SDS (d < 0.4), are shown (continuous lines).Blue: C bg .Green: C max .Red: σ.Cyan: z max .In the same colors, the mean values of each parameter for all the SDS are plotted (dashed lines).Note that, in graph (d), all z max values are grouped in intervals.

Fig. 6 .
Fig. 6.Fraction of well-retrieved profiles (d < 0.4) discriminated for the values of the Gaussian parameters C bg , C max , σ and z max .The different colors correspond to different values of the wavelength step δλ.