Inferring inherent optical properties and water constituent profiles from apparent optical properties

The BP09 experiment conducted by the Centre for Maritime Research and Experimentation in the Ligurian Sea in March 2009 provided paired vertical profiles of nadir-viewing radiances Lu(z) and downward irradiances Ed(z) and inherent optical properties (IOPs, absorption, scattering and backscattering coefficients). An inversion algorithm was implemented to retrieve IOPs from apparent optical properties (AOPs, radiance reflectance RL, irradiance reflectance RE and diffuse attenuation coefficient Kd) derived from the radiometric measurements. Then another inversion algorithm was developed to infer vertical profiles of water constituent concentrations, including chlorophyll-a concentration, non-algal particle concentration, and colored dissolved organic matter from the retrieved IOPs based on a bio-optical model. The algorithm was tested on a synthetic dataset and found to give reliable results with an accuracy better than 1%. When the algorithm was applied to the BP09 dataset it was found that good retrievals of IOPs could be obtained for sufficiently deep waters, i.e. for Lu(z) and Ed(z) measurements conducted to depths of 50 m or more. This requirement needs to be satisfied in order to obtain a good estimation of the backscattering coefficient. For such radiometric measurements a correlation of 0.88, 0.96 and 0.93 was found between retrieved and measured absorption, scattering and backscattering coefficients, respectively. A comparison between water constituent values derived from the measured IOPs and in-situ measured values, yielded a correlation of 0.80, 0.78, and 0.73 for chlorophyll-a concentration, non-algal particle concentration, and absorption coefficient of colored dissolved organic matter at 443 nm, respectively. This comparison indicates that adjustments to the bio-optical model are needed in order to obtain a better match between inferred and measured water constituent values in the Ligurian Sea using the methodology developed in this paper. © 2015 Optical Society of America OCIS codes: (010.4450) Oceanic optics; (010.5620) Radiative transfer; (290.0290) Scattering; (280.4991) Passive remote sensing. #237170 Received 30 Mar 2015; revised 2 Jul 2015; accepted 6 Jul 2015; published 22 Jul 2015 © 2015 OSA 27 Jul 2015 | Vol. 23, No. 15 | DOI:10.1364/OE.23.00A987 | OPTICS EXPRESS A987 References and links 1. H.R. Gordon and G.C. Boynton, “Radiance–irradiance inversion algorithm for estimating the absorption and backscattering coefficients of natural waters: homogeneous waters,” Appl. Opt. 36, 2636–2641 (1997). 2. H.R. Gordon and G.C. Boynton, “Radiance–irradiance inversion algorithm for estimating the absorption and backscattering coefficients of natural waters: vertically stratified water bodies,” Appl. Opt. 37, 3886–3896 (1998). 3. V. Sanjuan Calzado, D. McKee and C. Trees, “Multi and single cast radiometric processing and merging in the Ligurian Sea,” Optics of Natural Waters 2011 Saint Petersburg, September 2011. 4. J.R.V. Zaneveld, J.C. Kitchen, and C.M. Moore, “The scattering error correction of reflecting-tube absorption meters,” Proc. SPIE 2258 Ocean Optics XII, 44–55 (1994). 5. D. McKee, J. Piskozub, R. Röttgers, and R. Reynolds, “Evaluation and improvement of an iterative scattering correction scheme for in situ absorption and attenuation measurements,” J. Atmos. Ocean. Tech. 30, 1527–1541 (2013). 6. R. Röttgers, D. McKee, and S.B. Woźniak, “Evaluation of scatter corrections for ac-9 absorption measurements in coastal waters,” Methods in Oceanography 7, 21–39 (2013). 7. Z. Jin and K. Stamnes, “Radiative transfer in nonuniformly refracting media such as the atmosphere/ocean system,” Appl. Opt. 33, 431–442 (1994). 8. G.E. Thomas and K. Stamnes, Radiative Transfer in the Atmosphere and Ocean, 2nd Edition (Cambridge University, 2002) 9. B. Hamre, S. Stamnes, K. Stamnes and J.J. Stamnes, “C-DISORT: A versatile tool for radiative transfer in coupled media like the atmosphere-ocean system,” AIP Conference Proceedings 1531, 923 (2013). 10. B. Hamre, S. Stamnes, K. Stamnes, and J.J. Stamnes, “A versatile tool for radiative transfer simulations in the coupled atmosphere-ocean system: Introducing AccuRT,” Ocean Optics XXII, Portland, ME, 26-31 Oct. 2014. 11. C.D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic, 1994) 12. D. McKee, A. Cunningham, J. Slater, K.J. Jones, and C.R. Griffiths, “Inherent and apparent optical properties in coastal waters: a study of the Clyde Sea in early summer,” Estuarine, Coastal and Shelf Science 56(2), 369–376 (2003). 13. A. Morel and B. Gentili, “Diffuse reflectance of oceanic waters:II bi-directional aspects,” Appl. Opt. 32, 6864– 6879 (1993). 14. H.R. Gordon, “Can the Lambert-Beer law be applied to the diffuse attenuation coefficient of ocean water?,” Limnol. Oceanogr. 34, 1389–1409 (1989). 15. G. Fournier and J.L. Forand, “Analytic phase function for ocean water,” Proc. SPIE 2258 Ocean Optics XII, 194–201 (1994). 16. C.D. Mobley, L.K. Sundman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. 41, 1035– 1050 (2002). 17. T.L. Petzold, “Volume scattering functions for selected ocean waters,” Technical Report SIO 7278 Scripps Institute of Oceanography, San Diego, Calif. (1972). 18. W. Li, K. Stamnes, R. Spurr, and J.J. Stamnes, “Simultaneous retrieval of aerosols and ocean properties: A classic inverse modeling approach II. SeaWiFS case study for the Santa Barbara channel” , Int. J. Remote Sens.29, 5689– 5698 (2008). 19. Y.X. Hu, B. Wielicki, B. Lin, G. Gibson, S.C. Tsay, K. Stamnes, and T. Wong, “Delta-fit: A fast and accurate treatment of particle scattering phase functions with weighted singular-value decomposition least squares fitting,” J. Quant. Spectrosc. Radiat. Transf. 65, 681–690 (2000). 20. K. Ruddick, “DUE Coastcolour Round Robin Protocol, version 1.2,” Coastalcolour Round Robin, Oct. 2010. 21. C. Rodgers, Inverse Methods for Atmospheric Sounding (World Scientific, 2000). 22. D.S. Broomhead and D. Lowe, “Radial basis functions, multi-variable functional interpolation and adaptive networks,” Complex Systems 2, 321–355 (1988). 23. A. Bricaud, A. Morel, M. Babin, K. Allali, and H. Claustre, “Variations in light absorption by suspended particles with chlorophyll concentration in oceanic (case 1) waters: Analysis and implications for bio-optical models,” J. Geophys. Res.103, 31033–31044 (1998). 24. H. Loisel and A. Morel, “Light scattering and chlorophyll concentration in case 1 waters: a re-examination,” Limnol. Oceanogr. 43, 847–857 (1998). 25. A. Morel, D. Antoine, and B. Gentili, “Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function,” Appl. Opt. 41, 6289–6306 (2002). 26. M. Babin, A.D. Stramski, G.M. Ferrari, H. Claustre, A. Bricaud, G. Obelesky, and N. Hoepffner, “Variations in the light absorption coefficients of phytoplankton, nonalgal particles and dissolved organic matter in coastal waters around Europe,” J. Geophys. Res. 108(C7):3211 (2003). 27. M. Babin, A. Morel, V. Fournier-Sicre, F. Fell, and D. Stramski, “Light scattering properties of marine particles in coastal and open ocean waters as related to the particle mass concentration,”Limnol. Oceanogr. 28, 843–859 (2003). #237170 Received 30 Mar 2015; revised 2 Jul 2015; accepted 6 Jul 2015; published 22 Jul 2015 © 2015 OSA 27 Jul 2015 | Vol. 23, No. 15 | DOI:10.1364/OE.23.00A987 | OPTICS EXPRESS A988 28. A.M. Baldridge, S.J. Hook, C.I. Grove and G. Rivera, “The ASTER Spectral Library Version 2.0,” Remote Sens. Environ. 113, 711–715 (2009). 29. F. X. Kneizys, L. W. Abreu, G.P. Anderson, J.H. Chetwynd, E.P. Shettle, A. Berk, L.S. Bernstein, D.C. Roberson, P. Acharya, L.S. Rothman, J.E.A. Selby, W.O. Gallery, and S.A. Clough, “The MODTRAN2/3 report and LOWTRAN 7 model,” Phillips Laboratory, Hanscom AFB (1996). 30. L.G. Henyey and J.L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). 31. D. McKee, R. Röttgers, G. Neukermans, V. Sanjuan Calzado, C. Trees, M. Ampolo-Rella, C. Neil, and A. Cunningham, “Impact of measurement uncertainties on determination of chlorophyll-specific absorption coefficient for marine phytoplankton,” J. Geophys. Res.-Oceans 119, 9013–9025 (2014).


Introduction
Simultaneous measurements of the upward and downward irradiances, E u (z) and E d (z), or the nadir-viewing radiance L u (z) combined with E d (z), where z is the water depth, are among the most common data acquisition products in hydrologic optics.The ratios E u (z)/E d (z) and L u (z)/E d (z) are apparent optical properties (AOPs) that depend not only on the water and its constituents including embedded particles, but also on the illumination conditions, i.e. the angular distribution of the radiance throughout the water column.In contrast the inherent optical properties (IOPs), i.e. the absorption coefficient, a(z), the scattering coefficient, b(z), and the scattering phase function, p(z, Θ), where Θ is the scattering angle, depend exclusively on the properties of the water and its constituents including embedded particles, and are independent of the illumination.The radiative transfer equation (RTE) provides a connection between IOPs and radiances and irradiances, which are commonly measured, and from which AOPs including E u (z)/E d (z) and L u (z)/E d (z) are obtained.Hence, for a specified set of IOPs, the solution of the RTE provides the information required to determine the AOPs.Since the radiance reflectance defined as is obtained from L u (z) and E d (z), which are relatively easy to measure compared to the IOPs, a desirable goal is to infer the IOPs from measured profiles of the AOPs, in general, and from the irradiance reflectance R E (z) = E u (z)/E d (z) or the radiance reflectance R L (z), in particular.This nonlinear inversion problem is much more challenging than solving the forward problem involving the linear RTE. Gordon and Boynton [1] showed that it is possible to invert the radiances and irradiances to retrieve the IOPs for a homogeneous water column, and that an algorithm based on E u (z) and E d (z) profiles was more robust than one based on profiles of L u (z) and E d (z).An extension to make the method work for stratified waters was described in another paper by Gordon and Boynton [2].The IOPs of the water column consist of contributions from the pure sea water as well as from organic and inorganic suspended particles, i.e. pigmented (chlorophyll) and mineral particles, and colored dissolved organic matter or CDOM for short.For simplicity we will refer to these optically active constituents due to all embedded particles and dissolved matter as water constituents in the remainder of this paper.Once the IOPs are inverted from the irradiance reflectance R E (z) = E u (z)/E d (z) or the radiance reflectance R L (z) = L u (z)/E d (z), a further analysis can be performed to determine the contribution of each of the water constituents and the results will give a vertical distribution of water constituents throughout the water column.
In this paper we will describe our attempt to numerically implement the algorithm described in the papers by Gordon and Boynton [1,2], and the application to field measured L u (z) and E d (z) data for deep ocean waters.We then describe an algorithm developed to retrieve vertical profiles of water constituents from the retrieved IOP data.

AOP and IOP dataset
The AOP data derived from L u (z) and E d (z) measurements and paired IOP data that will be used in this paper are from the BP09 experiment which was conducted in the Ligurian Sea in March 13-26, 2009 by the NATO Centre for Maritime Research and Experimentation (CMRE) and collaborating institutions.Figure 1 shows 14 locations where measurements were conducted.

Ligurian Sea
Fig. 1.Location of 14 L u (z), E d (z), and IOP measurement conducted in the BP09 experiment used in this paper, see Table 1 for details.
The L u (z) and E d (z) measurements were performed with a Satlantic Hyper ProII radiometer package working in profiling mode.The free-falling profiler was configured with two hyperspectral radiometric sensors measuring downwelling irradiance E d (z) and upwelling radiance L u (z).Another set of radiometric sensors was mounted on the deck measuring the solar irradiance.The measurements were only performed in cloud free daylight conditions down to a maximum depth of 90 m.The measured L u (z) and E d (z) data have a spectral resolution of 3.3 nm from 400 to 750 nm and a depth resolution of 0.1 m.The L u (z) and E d (z) measurements were performed in the multicast mode in the top 10-15 m and then in a single-cast mode sampling further down in the water column until the irradiance was reduced to 1% of the surface irradiance, otherwise down to either the ocean bottom or 80 m due to limited cable length.The L u (z) and E d (z) data were normalized and rescaled in a post-processing step [3] to account for fluctuations in the solar irradiance at the ocean surface during the period of the measurements, and the multicast and single cast data were then merged into one profile.An interpolation between wavelengths was also performed during the data merging process and the merged L u (z) and E d (z) measurements have a spectral resolution of 5 nm from 400 nm to 750 nm.
The IOP measurements were performed with a WET Labs AC9+ and a BB9 spectrophotometer mounted on the University of Strathclyde IOP frame, in addition to a CTD.The AC9+ measures absorption, a(z), and attenuation, c(z), coefficients at 9 wavelengths centered at 412, 440, 488, 510, 532, 555, 650, 676 and 715 nm.The BB9 instrument measures the backscattering coefficient, b b (z), at the same wavelengths except that 555 nm is replaced by 595 nm.The backscattering coefficient at 555 nm was interpolated from the measurements at 532 and 595 nm.Calibration was performed using ultrapure water from a number of sources, and the measured IOPs agree well among these sources.The raw IOP dataset was then corrected for temperature, salinity and averaged over the multicast measurements in post-processing.A modified version of the commonly used "proportional" scattering correction method [4] was applied to the AC9+ data to correct for errors due to scattering in the absorption coefficients.However, recent studies by McKee et al. [5] and Röttgers et al. [6] showed that this scattering correction method leads to an underestimation of absorption coefficients which is more significant at longer wavelengths and under turbid water conditions.
The IOP [a(z), c(z), b b (z)] and E d (z), L u (z) measurements do match exactly in wavelength.IOP measurements are available at 412, 440, 488, 510, 532, 555, 650, 676 and 715 nm, while L u (z) and E d (z) measurements are available between 400 nm and 750 nm with 5 nm spectral resolution.In our study, we chose the L u (z) and E d (z) data whose wavelengths are closest to the wavelengths of the IOP data and we only used data from the six shortest wavelengths at 412, 440, 488, 510, 532 and 555 nm.We did not use data for longer wavelengths because (i) the radiative transfer model used in our IOP inversion algorithm does not take into account inelastic (Raman) scattering, and (ii) the uncertainties in the IOP measurements at the longer wavelengths are large.
The calibration data for the AC9 instrument show that the random measurement uncertainty for all wavelengths was ±0.003 m −1 for absorption coefficients and ±0.002 m −1 for attenuation coefficients.In the BP09 dataset, measured absorption coefficients a(z) and attenuation coefficients c(z) have overall average values of 0.06 m −1 and 0.7 m −1 , which means that the uncertainty is ±5% for absorption coefficients and less than ±1% for attenuation coefficients.The uncertainty in the L u (z) and E d (z) measurements are mostly due to surface waves, cloud cover and variation in solar zenith angle, but the post-processing of the L u (z) and E d (z) data reduced the uncertainty by merging multicast data and rescaling the data to account for the fluctuation of the solar irradiance.Therefore, the overall uncertainty of the L u (z) and E d (z) measurements is ±5% of the signal.The processed L u (z) and E d (z) data after multicast and merging are within the 95% confidence interval of the signal.

Radiative transfer model
Radiative transfer describes the interaction of the radiation field (i.e. the incoming solar radiation) with particles and molecules in the atmosphere and ocean, and provides a connection between IOPs and AOPs, such as the irradiance reflectance R E (z) = E u (z)/E d (z) or the radiance reflectance R L (z) = L u (z)/E d (z).In a slab geometry, the radiative transfer equation (RTE) can be written as [8]: where L(τ, µ, φ ) is the radiance, µ is the cosine of the polar angle, dτ(z is the single scattering albedo, and p(τ; µ , φ ; µ, φ ) is the phase function describing scattering from direction (µ , φ ) to (µ, φ ).The second term on the right hand side of the RTE is called the source term: S * (τ, µ, φ ) = ϖ(τ)F 0 4π p(τ; µ 0 , φ 0 ; µ, φ )e −τ/µ 0 , where F 0 is the solar irradiance, µ 0 and φ 0 are the cosine of the solar zenith and the corresponding azimuth angle, respectively.The third term on the right hand side is due to multiple scattering.
Solving the RTE for the azimuthally-averaged radiance L(τ, µ), and evaluating it at µ = 1 (nadir direction), one obtains the upward radiance L u (τ) in the nadir-viewing direction, and integrating L(τ, µ) multiplied by µ over all upward or downward directions, one obtains the upward irradiance E u (τ) and downward irradiance E d (τ): (3) Similarly, integrating the radiance L(τ, µ) over all directions, one obtains the scalar irradiance E 0 (τ): There are many methods available to solve the RTE.The forward radiative transfer model (RTM) used by Gordon and Boynton [1,2] was based on Monte Carlo simulations.In their original papers, Gordon and Boynton tested the algorithm at the top 5 meters of the water body.The Monte Carlo forward model is very computer-intensive and becomes impractical for operational use.To circumvent this problem we used AccuRT, which is a RTM for the coupled atmosphere-ocean system, based on the discrete-ordinate method [7][8][9][10].The AccuRT RTM is accurate, well-tested, reliable, and orders of magnitude faster than Monte Carlo simulations.Hence, we can easily extend the application to arbitrarily deep waters and test the algorithm for such cases.

IOP inversion algorithm
The IOP inversion algorithm employs the entire vertical profile of L u (z) and E d (z) to retrieve the profile of IOPs, a(z), b(z) and b b (z) iteratively, wavelength by wavelength.The theoretical basis of the IOP inversion algorithm described in the papers by Gordon and Boynton [1,2] can be briefly summarized as follow: 1.An approximate value of the absorption coefficient, a(z), is obtained from Gershun's Law [11]: where the average cosine µ(z) 2. Gordon and Boynton [2] proposed an approximation of the backscattering coefficient given by: b where and R E (z) ≡ E u (z)/E d (z). #237170 3. The scattering coefficient, b(z), can be computed from b b (z), if we know the scattering phase function p(cos Θ), as follows: bb (z) = 2π π π/2 p(cos Θ) sin(Θ)dΘ is called the backscatter ratio, or the backscatter fraction.
In the papers of Gordon and Boynton [1,2], both the E u -E d and the L u -E d algorithms were tested.Since the AOP data we have available are radiance reflectances R L (z) = L u (z)/E d (z), we will focus on the L u -E d algorithm and some modifications are required to make it more suitable for application to real (as opposed to synthetic) data.A flow chart of the L u -E d algorithm is shown in Fig. 2. To start the algorithm, a first guess of IOPs (a, b, and the scattering phase function) must be provided.We have found that a good first guess really helps reduce the number of iterations.In the absence of measured E u (z) and E 0 (z), Gordon and Boynton suggested that µ(z) can be replaced by µ w = cos θ w , the cosine of the solar zenith angle θ w measured in water, and that an approximate value of the upward irradiance, E m u (z), can be estimated from the measured radiance, L u (z), using the Q factor, Q(z) ≡ E u (z)/L u (z), and they suggested using Q(z) = π as a starting value.We have found that the Q factor may vary from π for isotropic scattering to about 10 for a strongly forwardscattering medium.Since the water column is inhomogeneous and the scattering in the water is anisotropic, using π as the Q factor at all depths may not be a good choice.In our algorithm we used an alternative approach to select a first guess, by using the radiance reflectance R L (z) given by Eq. ( 1) and the diffuse attenuation coefficient K d (z), given by McKee et al. [12] suggested a method of retrieving absorption and backscattering coefficients by combining the proposal made by Morel and Genitili [13] and Gordon [14].In 1989, Gordon [14] suggested a relation between the diffuse attenuation coefficient K d (z) and the absorption and backscattering coefficients: where θ w is the solar zenith angle measured in water.In 1993, Morel and Genitili [13] suggested a relation between the backscattering to absorption ratio (b b /a) and the radiance reflectance: By combining Eqs. ( 11) and ( 12), an estimation of a and b b can be obtained: McKee et al. [12] used Eqs.( 13) and ( 14) to retrieve IOPs from L u (z) and E d (z) measurements.
In our study, we found the estimation from these two equations to be insufficient for accurate inversion, however, they provide good starting values for our algorithm.From measured L u (z) and E d (z) profiles, we compute R L (z) from Eq. ( 1), K d (z) from Eq. ( 10), and θ w from the time and location of each measurement.Then we use Eqs.( 13) and ( 14) to provide a first guess of the IOPs, i.e. a (0) (z) and b (0) b (z).An appropriate scattering phase function must be provided to the radiative transfer model.The scattering properties of the water column are a combination of scattering by water molecules and by suspended particles.The fraction of scattering particles versus molecules varies with depth and the composition of the particles changes throughout the water column.Therefore, the scattering properties will depend also on wavelength.The angular variation of scattering by water molecules can be described by a Rayleigh phase function, but the scattering by suspended particles depends on the particle size distribution, refractive index, and wavelength.Hence it is difficult to find one phase function that fits all water types.In our study, we found that a good estimation of the backscattering ratio is critical to achieve a good estimation of scattering coefficients.The choice of the phase function is less important if the backscattering ratio is correct.In our algorithm we implemented a combined phase function of Rayleigh, Petzold and Fournier-Forand (FF) [15].We assumed that there are two types of particles in the water, algal particles (algae or phytoplankton) and non-algal particles (detritus, heterotrophic organisms and minerals).To describe the angular variation of the algal particle scattering, we used the FF phase function [15,16]: where ν = 0.5(3 − γ), and γ is the slope of the particle size distribution function (assumed to be a Junge or power law distribution), which typically varies between 3.0 and 5.0, u = 2 sin(Θ/2), 3(m−1) 2 , Θ is the scattering angle, and m is the refractive index.In a previous study, Li et al. [18] used m = 1.0686, and γ = 3.38, which correspond to a backscattering ratio of 0.0056.As noted by Mobley et al. [16], this choice of [m, γ] values is consistent with a certain mixture of living microbes.For non-algal particles, we chose the turbid water Petzold phase function, which was derived from measurements [17] made in San Diego Harbor and tabulated by Mobley [11].The water in San Diego Harbor is a mixture of microbes and sediments, with a backscattering ratio of 0.021.This phase function is a reasonable choice for the non-algal particles.
We used the Rayleigh phase function to represent molecular scattering by pure water: where f = 1−ρ 1+ρ , and ρ is the depolarization ratio, attributed to the anisotropy of the scatterer, which was set to be ρ = 0.0899.
The moment-fitting method of Hu et al. [19] was used to compute Legendre expansion coefficients χ ,PET and χ ,FF for the Petzold and FF scattering phase functions.The total IOP scattering phase function Legendre expansion coefficients χ are given by where f FF and f PET are the fractions of particles that scatter according to the FF and Petzold phase functions.These fractions were tuned so that the total backscattering ratio matches the measured one.
The backscattering ratio can be derived from IOP measurements, since bb (z) = b b (z)/b(z), and we have considered it to be known information.To find the fraction of the FF and Petzold phase function, we first subtract the scattering of the pure water and compute the backscattering ratio of suspended particles: from which a fraction of the FF phase function, f FF , can be derived to satisfy The fractions of the FF and Petzold phase function can be written as: The algorithm consists of the following steps: 1.The first guess IOPs (a, b and phase function) are used in the RTM to generate: 2. From the RTM generated radiance and irradiance values, new estimates of µ (1) (z) and Q (1) (z) are obtained.
3. Then Q (1) (z) will yield a new estimate of E m u (z) and K m v (z), and eventually a new estimate of a (1) (z).Once we have a new estimation of E m u (z), we can compute a new R m (z) and therefore a new X m (z).Meanwhile, by using the RTM-generated E u (z) and E d (z) we can compute R(z) and X(z), and the difference between X m (z) and X(z

The new estimate of the backscattering coefficient b
, where ε is a number between 0 and 1 that has been introduced to stabilize the algorithm.Once again b (1) (z) can be computed using Eq. ( 9).The new estimates of a(z) and b(z), along with the assumed phase function p(Θ) are then put back into the RTM to generate new radiance and irradiance values.
In the original paper of Gordon and Boynton, the algorithm stops when the residual error between measured and simulated values of E d and L u , δ (n) reaches a minimum, where They computed the difference between "measured" and simulated values of E d and L u directly and used the depth-averaged value as a convergence criterion in their algorithm because all the data were synthetic and they assumed a transparent sky condition, i.e. no atmosphere.However, the data used in this paper were based on measurements obtained in the presence of an atmosphere.Therefore, knowledge of the atmospheric condition when the measurements were performed is required if we were to use the same convergence criterion, but such information is not available in the dataset.We changed the convergence criterion by not computing the error of the "measured" L u (z) and E d (z) values directly but instead computing the error of the radiance reflectance R L (z), because our simulations showed that the radiance reflectance (an AOP) is insensitive to atmospheric conditions.We define a residual error between measured and simulated radiance reflectances, δ R L (n) , as The algorithm will stop when this error reaches a prescribed tolerance or the number of iterations reaches the maximum allowed value.

Water constituent inversion algorithm
The water constituent inversion algorithm employs the the multi-spectral IOPs, [a ], represents the wavelengths, to retrieve the water constituents iteratively, depth by depth.The IOPs retrieved from the radiance reflectance R L (z) data are the total water IOPs, which can be further divided into two parts: the IOPs due to pure water and the IOPs due to water constituents.We can subtract the pure water IOPs from the total, so that the remaining IOPs are solely contributed by water constituents, and these IOPs were used in our water constituent inversion algorithm to retrieve the water constituents.The vertical profile of water constituents is obtained after the inversion is completed.
The water constituents and the IOPs are connected by a bio-optical model.In our study, we assumed that there are two types of water constituent particles, algal particles and non-algal particles, in addition to the colored dissolved organic matter (CDOM), present in the water column.
Here we introduce the bio-optical model used in our study, the CCRR bio-optical model, where CCRR stands for CoastColor: Round Robin [20].In the CCRR bio-optical model, the spectral variation of the absorption and scattering coefficients for the algal particles is parameterized in terms of the chlorophyll-a concentration (CHL), and a Fournier-Forand phase function with backscattering ratio of 0.0056 is used as the scattering phase function.The spectral variation of the absorption and scattering coefficients for the non-algal particles is parameterized in terms of a mineral particle concentration (MIN), and we replaced the average Petzold phase function that was originally proposed in the CCRR bio-optical model with a turbid harbor Petzold phase function which has a slightly higher backscattering ratio of 0.021.The spectral variation of the CDOM absorption coefficients is given by an exponentially decreasing function.A more detailed description of the CCRR model can be found in Appendix A.
Based on the CCRR bio-optical model, the 3 water constituent parameters that will be retrieved from the IOPs are: chlorophyll-a concentration (CHL), the mineral particle concentration (MIN) and the CDOM absorption coefficient at 443 nm (CDOM).We developed an inversion algorithm based on a nonlinear least-squares optimal estimation with Levenberg-Marquardt regularization [18,21] that infers the 3 water constituents iteratively from the IOPs.
We first created a training synthetic dataset using the CCRR IOP model, where the 3 water constituents of the CCRR bio-optical model were randomly sampled in the following ranges: The output from the CCRR bio-optical model is the absorption, scattering and backscattering coefficients at 6 wavelengths: 412, 440, 488, 510, 532 and 555 nm.The dataset has 20,000 randomly sampled retrieval parameter [CHL, MIN, CDOM] and the derived IOPs [a i , b i , b bi ], where i = 1, . . ., 6, represents the 6 wavelengths.This dataset was used to train a Radial Basis Function Neural Network (RBF-NN) [22] that computes the 18 IOPs (a i , b i , b bi , i = 1, . . ., 6) from the 3 water constituents.The trained RBF-NN is used to replace the CCRR bio-optical model with the following single equation: where i takes the value of 1 to 18, N is the number of neurons, N in is the number of input water constituents (CHL, CDOM, MIN), which in this case is N in = 3, p k denote the input data, a i j , b, c jk and d i are the coefficients resulting from the training.The Jacobians K needed in the optimal estimation [see Eq. ( 26) below] are obtained by calculating the partial derivatives with respect to the input parameter p k : Our algorithm starts with a randomly picked initial guess of the 3 water constituents (CHL, CDOM, and MIN) at the top depth.Then these 3 constituents are input into the RBF-NN model to generate the IOPs, and an error between the true IOPs and the model generated IOPs at all wavelengths is computed.If the error is larger than a prescribed tolerance, the algorithm will update the 3 water constituents and start the next iteration.Otherwise, the algorithm stops, saves the 3 water constituents and proceeds to the next depth.At each iteration, the 3 water constituents are updated as follow: where x is the state vector, and in our algorithm x = [CHL, MIN, CDOM].γ is the Levenberg-Marquardt parameter.When γ = 0, Eq. ( 26) becomes standard Gauss-Newton optimal estimation, and when γ is a large number, Eq. ( 26) tends to the steepest gradient descent method.x a and S a are the a priori state vector and covariance matrix, respectively.K are the Jacobians and S m is the measurement error covariance matrix.y m and y i are the measured and simulated IOPs, respectively.
To speed up the algorithm, the first guess of the 3 water constituents is randomly picked only at the first depth, i.e. top of the ocean.For all depths below, we use the inverted 3 water constituents at the previous depth as a first guess.This approach appears to work well because the difference between each depth is only 1 meter and we assume that the 3 water constituents change gradually with depth.It also decreases the number of iterations needed to reach the preset error threshold, since the first guess is very close to the true value.

Synthetic data test of IOP inversion algorithm
To test the performance of the IOP inversion algorithm, we first generated a synthetic dataset at 488 nm using AccuRT, our RTM.The IOPs of the atmospheric gases are computed from a band model [29].The ocean can be stratified into a maximum number of 80 layers with arbitrary layer thickness, and a bottom reflection spectrum can also be added to simulate reflection from the ocean floor.
To generate synthetic L u (z) and E d (z) data we used an aerosol free standard U.S. atmosphere with solar zenith angle of 30 • coupled with a 500 m deep ocean.The top 80 m of the ocean was assumed to be stratified with a 2 m depth resolution.The remaining 420 m of water below 80 m was assumed to constitute a homogeneous layer, and the bottom of the ocean was assumed to be black (totally absorbing).The IOPs of the water column were assumed to vary with depth in the follow manner: We used a Henyey-Greenstein (HG) phase function [30] with asymmetry factor g = 0.9 which yields a backscatter ratio of 0.0229.Fig. 4. IOPs inverted from synthetic radiance reflectance R L (z) data at 488 nm.The upper panels show comparison between retrieved and synthetic IOPs, where filled circles (black) represent synthetic data, solid lines (blue) represent the inverted IOPs using the correct phase function (g = 0.9), and dashed lines(red) represent the inverted IOPs using incorrect phase function (g = 0.8).Note that when using correct phase function, the retrieved IOPs are so close to the synthetic data that the blue and black curves overlapped.The lower panels show the corresponding percentage errors.The IOPs inverted from the synthetic L u (z) and E d (z) data and the percentage error are shown in Fig. 4. When we use the same phase function as the one used to create the synthetic data, the inverted IOPs are very close to the synthetic data, with depth-averaged absolute percentage errors of 0.69%, 0.68% and 0.68% for absorption, scattering and backscattering coefficients, respectively.If we use a phase function that is different from the one used to create the synthetic data, we can still get reasonable inversion for absorption and backscattering coefficients, but the inverted scattering coefficient is incorrect.For example, when we used a HG phase function with g = 0.9 to create synthetic data and g = 0.8 in the inversion algorithm to invert the IOPs, the inverted absorption and backscattering coefficients have a depth-averaged absolute percentage error of 2.82% and 7.88%, respectively, but the inverted scattering coefficient had an error of 51.3%.These results agree with the original findings of Gordon and Boynton [1,2].The large error in the scattering coefficient is due to the use of an "incorrect" phase function, which yields a wrong backscattering ratio.The HG phase function used to create synthetic data has a backscattering ratio of 0.0229 (g = 0.9), while the "incorrect" phase function used to retrieve the IOPs has a backscattering ratio of 0.0507 (g = 0.8), which leads to an error in the backscattering ratio of 221%.This example shows that the algorithm produces better backscattering coefficient than scattering coefficient, i.e. it is more sensitive to b b than b, when an incorrect phase function is used.

Synthetic data test of water constituent inversion algorithm
The performance of the water constituent inversion algorithm was also tested with a synthetic dataset before we applied it to measured data.We created a validation synthetic dataset of 5,000 randomly sampled CHL, MIN, and CDOM combinations in the same way as we created the training dataset, as discussed in Section 4. The validation dataset is first used to test the accuracy of the neural network training.We take the 5,000 water constituents in the validation dataset and use the trained neural network to recompute the IOPs.Compared with the model generated IOPs, the correlation for the 18 IOPs [a(λ i ), b(λ i ) and b b (λ i ), i = 1, . . ., 6] had a minimum value of 0.99998, which means that the RBF-NN works very well, and is accurate enough to replace the CCRR bio-optical model.We then take the 5,000 validation data and put them into our water constituent inversion algorithm to test its performance.The 18 model generated IOPs (a(λ i ), b(λ i ) and b b (λ i ), i = 1, . . ., 6), were used to retrieve water constituents, and the results were then compared with the #237170 modeled data as shown in Fig. 5.We computed the percentage error for each retrieved combination of CHL, MIN and CDOM compared with the modeled data, and the errors are indicated by the color of the dots in Fig. 5. Except for a few cases in which the algorithm failed, most of the retrieved data match the model data very well.The correlation between retrieved and model data had R 2 -values of 0.9998, 0.9946 and 0.9995 for CHL, MIN and CDOM, respectively.The L u (z) and E d (z) data processing (merged multi-cast data corrected for solar irradiance fluctuations) and pairing with IOP data are still underway, and we have a total number of 14 measurements available.The time, location and mean ocean depth of all 14 measurements are shown in Table 1.There are six cases for which the measurements were obtained in shallow waters where the ocean depth is less than 30 m; 2 cases were obtained in very deep waters where the ocean is more than 2,000 m deep, and the remaining six cases were obtained in waters where the ocean depth varies from 50 m to 300 m.The paired IOP measurements, in general, has temporal difference of less than 1.5 hours with the L u (z) and E d (z) measurements.When applying the algorithm to the measurements, we assumed a 14 layer aerosol free U.S. standard atmosphere since the atmospheric condition has minimal impact on the retrieved IOPs.The ocean was stratified based on the mean ocean depth and the depth at which continuous nonzero L u (z) and E d (z) measurements were obtained.The layers were selected in the following way: (i) the layers were set with 1 m layer depth from the top of the ocean to the maximum depth where continuous non-zero L u (z) and E d (z) measurements were obtained; (ii) if the mean ocean depth is less than 200 m, we expanded the lowermost layer that was set in (i) down to the ocean bottom by changing the lower boundary of that layer; (iii) if the ocean depth is greater than 200 m, we first expand the lowermost layer that was set in (i) down to 200 m, then added another layer of pure water from 200 m to the bottom of the ocean.For example, ST28 has a mean ocean depth of 2425 m, and non-zero L u (z) and E d (z) measurements were obtained down to 56 m, where the radiance reaches the 1% light level of the radiance at the surface.So in the top 56 layers, the depth of the each of the first 55 layers was set with 1 m and the lower boundary of the 56th layer was set to 200 m which yields a layer depth of 145 m.Finally, we added a 2225 m thick layer of pure water, so the total ocean depth agrees with geographic data, yielding a total of 57 layers for this case.We expand the last layer in order to avoid abrupt changes in L u (z) and E d (z) profiles.The bottom of the ocean was assumed to reflect light according to the reflectance spectrum of loamy sand [28], although for these deep waters reflection from the ocean floor is negligible.Fig. 6.Retrieved IOPs (red) of ST15 compared with in-situ measurements (blue).The top 6 panels, from left to right, show comparisons of the absorption coefficients at 412, 440, 488, 510, 532 and 555 nm, respectively, whereas the middle and bottom 6 panels show the same for the scattering and backscattering coefficients.
The IOP inversion algorithm was then applied to the measured L u (z) and E d (z) datasets individually to retrieve IOP profiles wavelength by wavelength.However, as noted by Gordon and Boynton, the integral in Eq. ( 8) needs to be carried out to a depth where dR E (z)/dz → 0 in order to provide a good estimate of b b (z).This requirement is not satisfied unless the water is deep enough.So we applied the algorithm focusing on the deep water cases, i.e.ST10, ST13, ST14, ST15, ST19, ST27 and ST28.We excluded ST18 because L u (z) and E d (z) measurements were obtained only in the top 15 m of the water column.
For all the deep water cases, the algorithm required 15 iterations on average before the prescribed tolerance (δ R L <0.001) was satisfied, and the radiance reflectances R L (z) derived from L u (z) and E d (z) values that were reconstructed by our RTM using retrieved IOPs matched the ones derived from the L u (z) and E d (z) measurements with a depth averaged absolute percentage error (PE) given by of maximum 0.23%.Table 2 shows the number of iterations (i) and the absolute percentage error of the radiance reflectances R L (z) at each wavelength for the deep water cases.The retrieved IOPs generally agree well with the measured IOPs for the deep water cases.Table 3 summarizes the correlation and bias given by bias  We used different colors to distinguish between wavelengths in Fig. 8 in order to show how the algorithm performs at different wavelengths.The algorithm introduced a small bias to the retrieved IOPs with a tendency to overestimate absorption coefficients.This overestimation of the absorption coefficient occurs at all wavelengths, but has a larger impact at shorter wavelengths, i.e. 412 nm and 440 nm.In extreme cases, the algorithm overestimates the absorption coefficient at 412 nm by 50%.Since the absorption coefficient was determined by Gershun's Law, we suspect that even though the radiance reflectances R L (z) fit the measurements very well, with depth averaged absolute percentage error less than 0.23%, the irradiance reflectance R E (z) ≡ E u (z)/E d (z) may still have large errors, which could be traced to an inaccurate phase function.As pointed out by Mobley et al. [16], even for the same backscattering ratio, a different shape of the phase function can yield quite different L u (z) and E d (z) profiles.This issue could be addressed by including E u (z) in the radiometric measurements.Another possible reason, as discussed in Section 2, is the underestimation in the measured absorption coefficients due to imperfect scattering correction.This uncertainty may also explain some of the apparent overestimation in absorption coefficients.The errors in the retrieved scattering and backscattering coefficients seem to be wavelength independent, and a possible reason is the error in the evaluation of X(z).Even for the deep water cases, the requirement that dR E (z)/dz → 0 at z max is not completely satisfied, so there are still errors in X(z).

Retrieval of water constituents from IOPs
The water constituent inversion algorithm (IOPs → {CHL, CDOM, MIN}, see Section 4) was then applied to the BP09 measurement data.There are two sets of IOP data available for retrieval of water constituent profiles.One set is the 14 in-situ measured IOP data, and the other set consists of the IOP data retrieved from L u (z) and E d (z) measurements by our IOP inversion algorithm for the 7 deep water stations analyzed in this paper.Since the retrieved IOP profiles contain some errors, we first applied the water constituent inversion algorithm to the in-situ measured IOPs directly to examine the performance of the algorithm.We applied the algorithm to all the measured IOP profiles to retrieve profiles of the 3 water constituents, CHL, MIN and CDOM.The retrieved water constituents were compared with the in-situ measured values.There are in-situ measurements of CHL, MIN and CDOM available in the BP09 dataset.The MIN values were measured only in the top levels, i.e. water samples were taken 1 m below the ocean surface.The CHL and CDOM values were measured at both this top level and at a few additional depths (7 m, 14 m, 30 m, 50 m and 125 m) depending on location.Figure 9 shows a comparison of retrieved CHL, MIN and CDOM from the in-situ measured IOPs against the measured values for all the stations.The correlation between the two datasets has R 2 values of 0.80, 0.78 and 0.73 for CHL, MIN and CDOM respectively.This low correlation is clearly due quite visible in the MODIS standard CHL a product (Fig. 11).The measurements at ST27 and ST28 were performed right in the area of the algae bloom where a higher chlorophyll concentration can be found from the MODIS image, and our retrieved water constituents show clearly that the water column was dominated by chlorophyll.The vertical distribution also matches the measured CHL data well with a peak in the chlorophyll concentration around 10 m.The measurements at ST10, ST13, ST14, ST15, and ST19 were performed in less productive areas, where the chlorophyll concentration was found to be lower, and the water column is more dominated by non-algal particles.The water constituent profiles retrieved from inverted IOPs (dashed lines in Fig. 10) were slightly different from those retrieved from measured IOPs due to errors in the inverted IOPs.

( 1 )
b (z) is obtained by changing the old values of b (0) b (z) by a small amount ∆b b (z), where ∆b b (z) = a(z)∆X(z).

Fig. 5 .
Fig. 5. Retrieved CHL, MIN and CDOM compared with model data for the synthetic dataset.The color of the dots shows the percentage error in the retrieved data.

Fig. 7 .
Fig.7.Retrieved IOPs (red) of ST28 compared with in-situ measurements (blue).The top 6 panels, from left to right, show comparisons of the absorption coefficients at 412, 440, 488, 510, 532 and 555 nm, respectively, whereas the middle and bottom 6 panels show the same for the scattering and backscattering coefficients.

Fig. 8 .
Fig. 8.Comparison between retrieved and measured IOPs that combined all depth, wavelengths and stations.The 3 panels, from left to right, show comparisons of absorption, scattering and backscattering coefficients, respectively.The color of the dots indicates wavelength as shown on the color bar.

Fig. 9 .
Fig. 9. Retrieved CHL, MIN and CDOM from in-situ measured IOPs compare with measured values for all stations.

Fig. 11 .
Fig. 11.Standard MODIS CHL a product of the Ligurian Sea area on 17 March 2009.

Table 1 .
Time, location and mean ocean depth of the 14 measurements.The * symbol indicates the deep ocean cases tested in our IOP inversion algorithm.

Table 2 .
Number of iterations (i) needed to invert the IOPs from L u (z) and E d (z) measurements and absolute percentage error (PE) of the radiance reflectances R L (z) for deep water cases.

Table 3 .
Correlation and bias of retrieved IOPs for deep water cases.