## Abstract

We develop two algorithms for determining two inherent optical properties (IOPs) from radiometric measurements in vertically homogeneous waters. The first algorithm is for estimation of the ratio of the backscattering to absorption coefficients from measurements of only the vertically upward radiance and the downward planar irradiance at depths where the light field is in the asymptotic regime. The second algorithm enables estimation of the absorption coefficient from measurement of the diffuse attenuation coefficient in the asymptotic regime after use of the first algorithm. Multiplication of the two estimates leads to an estimate for the backscattering coefficient. The algorithms, based upon the use of a simplified phase function and the asymptotic eigenmode, are shown to potentially provide good starting conditions for iteratively determining the absorption and backscattering coefficients of a wide variety of waters. The uncertainty in the estimates defines a subspace for IOPs that may reduce ambiguity in such iterative solutions. Because of the ease of estimating the backscattering to absorption ratio from in-water measurements, this IOP deserves further investigation as a proxy for biogeochemical quantities in the open ocean.

© 2011 Optical Society of America

## 1. Introduction

The estimation of ocean inverse optical properties (IOPs) from radiometric measurements is a long-standing and important problem in ocean optics [1], forming the basis of remote sensing of ocean water constituents and biogeochemical properties. The solution for one or more IOPs from apparent optical properties (AOPs) is an example of such a problem. Using just the knowledge of the in-water light field, the inverse problem for IOPs can be solved by closed-form (explicit) methods or iterative (implicit) methods [2, 3].

Explicit methods do this by manipulation of the radiative transfer equation (RTE) or analytic solutions to derive formulas that estimate IOPs from AOPs. The explicit method is illustrated by the early work of Gordon *et al.* [4] who showed, assuming quasi-single scattering [5] and employing Monte Carlo simulations, that remote sensing reflectance was strongly correlated with the ratio of the backscattering to absorption coefficients, *b _{b}*/

*a*. Explicit algorithms also can be derived from solutions to the RTE expressed as a linear combination of eigenfunctions [6], or by carrying out a large number of radiative transfer simulations in order to parameterize relationships between IOPs and AOPs [7, 8].

Implicit methods estimate IOPs by repeatedly solving the RTE, starting with an initial estimate or guess of the IOPs. At each iteration, radiometric measurements (or AOPs derived therein) are compared to the RTE solution and then an objective function is computed that expresses the difference between the RTE solution and the measured values. The IOPs are then adjusted to reduce the objective function. Iteration is complete when the objective function is reduced to an acceptably small value. Implicit approaches to estimating absorption and backscattering in homogeneous waters have been developed that estimate IOPs from surface illumination and depth profiles of easily-measured radiometric quantities, while allowing for general multiple scattering as well as inelastic scattering [9–11].

Two inversion algorithms are developed here, one for estimating the *b _{b}*/

*a*ratio and the other for the absorption coefficient

*a*. For estimating

*b*/

_{b}*a*, we use a hybrid approach where we iteratively solve a system of equations derived from an explicit eigenfunction formulation of the RTE. Our algorithm for

*b*/

_{b}*a*requires one fitted parameter

*F*that corrects for an oversimplification of our assumed phase function. That phase function enables us to use a simplified solution of the radiative transfer equation which depends only on

*a*and

*b*. We also assume that the optical measurements are made in the asymptotic regime. (The depths at which the light reaches an asymptotic distribution depends on the IOPs and the surface illumination and occurs at shallower depths when the incident light is more diffuse [3, 12, 13].) With the values of

_{b}*F*and

*b*/

_{b}*a*, our second algorithm can be used to estimate

*a*from measurement of the in-water downward diffuse attenuation coefficient. Combination of the two estimates directly leads to an estimate of the backscattering coefficient

*b*. Alternatively,

_{b}*a*can be estimated from

*b*/

_{b}*a*, measurements of

*b*at one or more wavelengths, and a model for the spectral behavior of

_{b}*b*(

_{b}*λ*).

The resulting estimates of *a* and *b _{b}* can be used, for example, as the initial starting estimate, instead of a guess, for any implicit method, such as one using a nonlinear global optimization technique [11]. Furthermore, uncertainty in estimates from our algorithms can be used as inversion constraints [11, 14, 15].

## 2. Radiative transfer equation

For our investigation the radiance *L*(*z*,*μ*,*φ*) for *z* ≥ 0 is taken to be a function of the direction cosine *μ*, measured with respect to the depth *z* in the water, and *ϕ*, the azimuthal angle measured in the plane perpendicular to the *z*-axis. The water is taken to be spatially uniform and deep enough that bottom effects can be neglected. For scattering given by the phase function *β*̃(*μ*_{0}) = *β*̃(*μ*′,*φ*′,*μ*,*φ*) for *μ*_{0} the cosine of the angle between {*μ*′,*φ*′} and {*μ*,*φ*}, the radiative transfer equation for the semi-infinite medium is [3]

*c*=

*a*+

*b*is the sum of the absorption and scattering coefficients. Notice that this formulation and subsequent algorithm development ignores inelastic scattering. As a result, the algorithm developed here will not be applied to spectral regions in which inelastic scattering can make a significant contribution to the upward radiance [10, 16, 17].

Because the objective here is to develop algorithms for measurements done with a planar irradiance detector and a radiance detector aimed in the vertically upward direction for which *μ* = −1, there is no azimuthal dependence of the data so the analysis will be done for the azimuthally averaged radiance

We will derive an approximate equation for determining *b _{b}*/

*a*from data for the vertically upward radiance

*L*(

_{u}*z*) =

_{m}*L*(

*z*, −1) at measurement depth

_{m}*z*and the downward planar irradiance

_{m}*E*(

_{d}*z*),

_{m}*b*/

_{b}*a*that is expressed in terms of the in-water remote sensing ratio, Following that we will obtain the algorithm for determining

*a*from the estimated value of

*b*/

_{b}*a*and multiple measurements of

*E*(

_{d}*z*). Both algorithms will be developed with a value of

_{m}*F*that is obtained by fitting data from previous calculations or experiments to situations where

*b*and

_{b}*a*are known.

There are two assumptions needed to develop the algorithms. The first is that water beneath measurement depth *z _{m}* is sufficiently spatially uniform and deep enough that

*z*is approaching the so-called “asymptotic regime” where the angular distribution of the radiance will not perceptibly change with increasing depth. Another interpretation of this first assumption is that

_{m}*z*must be deep enough that the directional nature of the incident surface illumination does not affect the directional nature of the underwater light field. The second assumption is that the angular distribution of the volume scattering function is so strongly forward that an “isotropic plus delta forward” synthetic scattering phase function model can be used,

_{m}*b̃*=

_{b}*b*/

_{b}*b*is the backscattering ratio and

*F*is the adjustment factor to be determined from numerical simulations or experimental data. Our phase function satisfies the proper normalization condition A primary advantage of this phase function is that it leads directly to a radiative transfer equation that depends only on

*a*and

*b*, but not

_{b}*b*, which has been shown to be a significant advantage for remote sensing applications [4, 18].

The phase function of Eq. (7) is a slight extension of the simplest form of the classic “delta-Eddington” scattering model [19] to which it degenerates if *F* = 1. Although this elementary model for the scattering phase function does not account for the directional dependence of backscattering, such a dependence has been shown to give errors in underwater radiances and planar irradiances on the order of 10% or less [20].

Substitution of the phase function of Eq. (7) into Eq. (3), followed by changes in the variables to obtain a dimensionless optical distance,

results in the radiative transfer equationIt is the parameter *F* that first must be obtained from either Hydrolight [21] calculations with known values of *b _{b}*/

*a*or from experimental data in which both

*r*(

_{rs}*z*) and

_{m}*b*/

_{b}*a*have been measured. Once obtained, the resulting

*F*subsequently can be used in the algorithm for determining

*b*/

_{b}*a*for waters where no or limited IOP measurements are available. Notice that with this model the objective is to avoid the dependence on

*b*by having the scattering depend instead on

*b*, which is consistent with the quasi-single-scattering approximation [3–5].

_{b}## 3. Algorithm for backscattering-to-absorption ratio

The initial step in the derivation of the algorithm for *b _{b}*/

*a*follows by combining Eq. (10) written for

*μ*and for −

*μ*in a way that the derivative with respect to

*τ*can be eliminated [22]. This is done by first multiplying the equation for

*μ*by

*μL*(

*τ*,−

*μ*) and the one written for −

*μ*by

*μL*(

*τ*,

*μ*), integrating both equations over −1 ≤

*μ*≤ 1, and subtracting the results to obtain

*μ*, integrate over −1 ≤

*μ*≤ 1, and rearrange the result to obtain the so-called diffusion equation,

*G*(

*τ*)/d

*τ*over

*τ*from an arbitrary measurement depth

*τ*or

_{m}*z*to infinity, where the light field vanishes, allows us to write

_{m}*E*(

_{d}*z*)/2

_{m}*π*]

^{2}and with Eqs. (5) and (12), Eq. (16) can be written as

The integrals in Eq. (17) involve radiance data that are not available because *L*(*z _{m}*,

*μ*) is not measured, so we introduce two new radiometric functions

*r*(

_{rs}*z*) value as

_{m}Equations (18) and (19) can be evaluated by invoking the second assumption for the algorithm that the measurements are made at any *asymptotic* depth *z _{as}*. With the approximation of

*L*(

*z*,

_{as}*μ*) to evaluate Ω(

*z*) and Λ(

_{m}*z*), from Eq. (21) the final form for the

_{m}*b*/

_{b}*a*algorithm is

To evaluate *L*(*z _{as}*,

*μ*) in Eqs. (23) and (24), an analytic representation can be used for the radiance at depths where it is asymptotic. The eigenmode expansion method [23,24] shows that the asymptotic radiance far from the surface or bottom consists of only the dominant eigenmode [25] corresponding to the dominant eigenvalue

*ν*

_{0},

*A*(

*ν*

_{0}) is a constant that depends on the surface illumination.) Substitution of Eq. (25) into Eq. (10) gives

*ν*

_{0}≤ ∞, the eigenfunction

*ϕ*(

*ν*

_{0},

*μ*) can be normalized with so from Eqs. (27) and (28),

*ν*

_{0}as the positive root of the transcendental equation From Eq. (11) it follows that

*ϖ*depends on

*F*, so

*ν*

_{0}determined from Eq. (30) also does. Another very useful equation from the eigenfunction method is obtained by integrating Eq. (27) over −1 ≤

*μ*≤ 1 to obtain, after use of Eqs. (11) and (28),

*L*(

*z*,

_{as}*μ*) from Eqs. (26) and (29) it follows from Eq. (18) that

To summarize, for a specified value of *b _{b}*/

*a*,

*F*must be determined consistent with Eqs. (11), (22), (30), (32), and (33).

## 4. Algorithm for absorption and backscattering coefficients

The algorithm for determining the absorption coefficient *a* ultimately follows from the Gershun equation

*K*(

_{d}*z*) is measured instead of

_{m}*K*(

_{E}*z*), it is again necessary to invoke the assumption that measurements are made at depths

_{m}*z*→

_{m}*z*where

_{as}*K*(

_{E}*z*) =

_{as}*K*(

_{d}*z*). Thus an algorithm for determining

_{as}*a*can be obtained from where

*K̂*(

_{d}*z*), an estimate of the diffuse attenuation coefficient, follows from Eq. (35) as

_{as}*K*(

_{d}*z*),

_{n}*n*= 1 to 3, where the depths

*z*are such that the radiance is approaching the asymptotic regime, an estimate of

_{n}*K̂*(

_{d}*z*) can be obtained from [6, 26]

_{as}*z*

_{2}= 2

^{−1}(

*z*

_{1}+

*z*

_{3}).

To estimate *μ*̄(*z _{as}*) in order to obtain the absorption coefficient from Eq. (37), it follows from Eqs. (26), (28), and (31) that

*b*/

_{b}*a*has been obtained from Eq. (22) and used, along with the previously determined

*F*, to compute the

*ν*

_{0}.

A second approach for determining the absorption coefficient from knowledge of *b _{b}*/

*a*and

*F*is to take the ratio of measurements in the asymptotic regime at two locations

*z*

_{1}and

*z*

_{2}and to use

*z*

_{1}and

*z*

_{2}in the asymptotic radiance regime this leads to the approximation

*L*(

_{u}*z*), however, this algorithm is expected to be prone to numerical difficulties. Another possibility, of course, is if

_{as}*b*is measured or can be estimated at the wavelengths of interest, then

_{b}*a*can be directly estimated from

*b*/

_{b}*a*as

Finally, the estimated values of *b _{b}*/

*a*from Eq. (22) and

*a*from Eq. (41) or (43) can be combined to give an estimated value of

*b*from

_{b}## 5. Data and methods

#### 5.1. Synthetic data sets

We assessed the ability of our extended delta-Eddington model to estimate the scattering phase function using synthetic data. Radiometric quantities *L _{u}* and

*E*were simulated at 490 nm using a simple Case 1 chlorophyll-based bio-optical model embedded in Hydrolight [21] to generate inherent optical properties. Simulations were performed using chlorophyll concentrations

_{d}*Chl*(0.01, 0.03, 0.1, 0.3, 1, 3, 10 mg m

^{−3}) and three zenith angles (15°,30°,60°), yielding 7x3 = 21 total simulations. While simulation output was requested every 5 m from 0 m to 50 m, radiometric data for three output depths (5, 10, 30 m) from each simulation was included in this synthetic data set, resulting in 63 simulated data points.

One synthetic data set was created allowing Hydrolight to select a Fournier-Forand scattering phase function *β*̃_{FF} [27] based on a specified particle backscattering ratio *b̃ _{bp}* =

*b*/

_{bp}*b*, which resulted in a distinct scattering phase function selection for each chlorophyll concentration. To assess the ability of our algorithm to retrieve

_{p}*b*/

_{b}*a*in the absence of measurement noise and natural variability, the entire first synthetic data set was used both as a

*training set*to determine the fitted parameter

*F*and also as a

*validation set*to compare values of

*b*/

_{b}*a*generated by the Case 1 bio-optical model to those estimated from corresponding simulated values of

*r*=

_{rs}*L*/

_{u}*E*.

_{d}A second, somewhat independent, synthetic validation set was created using the same cross product of chlorophyll concentrations, depths, and zenith angles, but using a single Petzold scattering phase function *β*̃_{P} [28] for all simulations of radiometric quantities. This second synthetic data set was used solely as a validation set to test the sensitivity of our algorithm’s estimates of *b _{b}*/

*a*to any assumed phase function effects in the parameter

*F*fitted using the first training set.

#### 5.2. NAB08 data set

We also tested our algorithm’s performance using depth profiles of *L _{u}*(

*z*),

*E*(

_{d}*z*),

*a*(

*z*), and

*b*(

_{b}*z*) from a calibration campaign carried out as part of the 2008 North Atlantic Bloom Experiment (NAB08) in support of long term (51 d) autonomous physical, optical, and radiometric measurements from a Lagrangian float [29]. The data were obtained at eight process cruise stations during the spring bloom of May 2008, representing two distinct phytoplankton bloom communities: a large diatom-dominated community (

*Chl*∼ 2 – 5 mg m

^{−3}) during the bloom and mixed ciliate/picoeukaryote community (

*Chl*< 1 mg m

^{−3}) following the bloom peak.

### 5.2.1. NAB08 data collection and processing

Temperature, salinity, pressure (CTD), and bio-optical profiles to approximately 80 m depth were performed during a cruise on the R/V Knorr from 1–22 May 2008. A Satlantic free-falling optical profiler was used to measure downward spectral irradiance, *E _{d}* (

*z*,

*λ*), and upward spectral radiance,

*L*(

_{u}*z*,

*λ*) at 3.3 nm increments from 350 to 800 nm with a spectral accuracy of 0.3 nm and a spectral resolution of 10 nm. Profiles were taken within 1.5 h of local noon. Radiometric data was processed using 1 m bins with ProSoft 8-RC5 software (Satlantic, Inc.) to determine of

*E*(

_{d}*z*,

*λ*),

*L*(

_{u}*z*,

*λ*), and

*K*(

_{d}*z*,

*λ*).

*In situ* IOPs were obtained within 15 minutes of radiometric measurements. An AC-9 measured *a*(*z*); a BB2F measured *b _{b}*(

*z*) at 470 and 700 nm (both WET Labs, Inc.); a Sea-Bird Electronics SBE25 CTD measured temperature, salinity, and pressure. The instruments were factory calibrated prior to field deployment. Manufacturer recommended protocols were used to track instrument calibration during the process cruise. Daily clean water calibrations were conducted whenever possible; however, sampling schedules did not consistently allow for daily calibration. Under these cases, the most recent water calibration was used. The maximum period without a calibration was three days. The drift during this study was negligible (< 4%) and the precision of the AC-9 data was ±0.01 m

^{−1}, which corresponds to an average uncertainty of about 17% in the waters under consideration. Absorption data were subsequently corrected for temperature and salinity [30] and the absorption coefficient was corrected for scattering using the wavelength-dependent method [31]. Absorption and backscattering data were binned into 1 m intervals and then averaged within each bin. Backscattering at wavelengths other than those measured was estimated by fitting the particulate fraction

*b*using a power-law

_{bp}*b*(

_{bp}*λ*) =

*b*(

_{bp}*λ*

_{0})(

*λ*

_{0}/

*λ*)

*wavelength dependence [14]. Measured IOPs were typically homogeneous from 0 to 30 m (e.g., Fig. 1A). Some IOP profiles showed sharp changes deeper than 30 m, corresponding to the bottom of the mixed layer; outliers from these regions with different IOPs appeared as outliers in IOP estimates.*

^{γ}For each profile, radiometric and IOP data were sampled at eight depths (5, 10, 15, 20, 25, 30, 35, 40 m) at AC-9 wavelengths (412, 440, 488, 510, 532 and 555 nm). For wavelengths ≥ 555 nm we expect our algorithms to work poorly because our formulation of the radiative transfer equation (Eq. 1) does not include inelastic scattering [10, 17]. The 40 m depth was omitted for two profiles where radiometric measurements were below the noise level of the radiometer. Example contemporaneous radiometric and IOP profiles at 488 nm for a single station are shown in Fig. 1; *b _{b}*(488) is estimated from

*b*(470) and

_{b}*b*(700) as described above and

_{b}*b*/

_{b}*a*is computed from

*a*(488) and the estimated

*b*(488).

_{b}When processing the measured data, no attempt was made to incorporate possible instrumentation inaccuracies. Radiometric measurements were discarded for *L _{u}*(

*z*) < 10

^{−4}

*μ*W cm

^{−2}sr

^{−1}or

*E*(

_{d}*z*) < 10

^{−2}

*μ*W cm

^{−2}.

### 5.2.2. Cross-validation approach

We used the NAB08 data set to both determine the fitted parameter *F* as well as test the algorithm. We were also interested in answering the question: “How many calibration profiles are needed to determine *F* for a specified accuracy?” To answer this question as fairly as possible, we employed a cross-validation approach [32] to estimate error in our predictions of IOPs from *F*, *r _{rs}* and

*K*(

_{d}*z*). Cross-validation is a resampling technique used to estimate how accurately a predictive model will perform in practice. Our cross-validation scheme iteratively partitions the

_{as}*n*= 8 profiles into a

*training set*used to determine the fitted parameter

*F*and a

*validation set*used to estimate the IOPs. A training set consists of a

*k*-element subset of profiles,

*k*= 1 ... 7, chosen from the set of eight profiles, resulting in ${M}_{k}=\left(\begin{array}{c}n\\ k\end{array}\right)$ possible training set combinations for a given

*k*. A validation set of profiles is simply the complementary set of profiles for any training set combination, assuring that training data is never used as validation data. Examples for profile numbers included for

*k*= 1,2 are shown in Table 1.

Traditional cross-validation yields a nearly unbiased estimate of prediction error by removing one data point at a time from the data set used for prediction. A weakness of our approach is that the number of members of each set and the number of sets vary with each iteration *k*.

#### 5.3. Determination of fitted parameter F

For each training set, *F* was found using the contemporaneous values of *r _{rs}* =

*L*/

_{u}*E*,

_{d}*b*, and

_{b}*a*available from Hydrolight simulation or field measurement. The problem was solved as a system of three equations

**G**=

_{F}**0**with three unknowns

**X**, where

_{F}**X**that minimizes the cost function

_{F}*g*=

_{F}**G**

_{F}^{T}

**G**, where the first two elements of

_{F}**G**correspond to Eqs. (22) and (30), and the final element provides the necessary third constraint by matching the estimated and measured values of

_{F}*b*/

_{b}*a*:

**X**

_{F}^{0}, was chosen as Finally, using the results for all the depths of given profile and wavelength, a simple linear least squares fit was made to estimate

*F*as a linear function of

*b*/

_{b}*a*,

#### 5.4. Determination of b_{b}/a

The validation sets of profiles were used to test the algorithm. The linear equation Eq. (49) for the fitted parameter *F*(*b _{b}*/

*a*) determined from the training set and

*r*=

_{rs}*L*/

_{u}*E*computed from simulated or field radiometric data were used to estimate

_{d}*b*/

_{b}*a*. The contemporaneous measured or simulated values of

*b*and

_{b}*a*were reserved for subsequent error analysis (see section 5.5). Similar to the approach used to determine fitted parameter

*F*in section 5.3, this problem was solved as a system of two equations

**G**=

**0**and two unknowns

**X**, where

*g*=

**G**

^{T}

**G**, where the first two elements of

**G**again correspond to Eqs. (22) and (30):

**X**

^{0}, was chosen as

#### 5.5. Error analysis

Assessment of systematic biases and typical uncertainty between estimated IOPs and measured IOPs is made using relative percent difference *ψ* and its absolute value |*ψ*|, respectively. In the formulation of the relative percent differences, measured IOP data are considered the reference values, so for a paired estimate and measurement the relative error is

*b*/

_{b}*a*,

*a*, or

*b*,

_{b}*IOP*

_{i,est}(

*z*) is an estimated value, and

_{j}*IOP*

_{i,meas}(

*z*) is the corresponding measured value for a given validation set

_{j}*i*= 1 ...

*M*at depth

_{k}*z*. Similarly, the average absolute percent difference is computed with |

_{j}*ψ*(

_{i}*z*)|. Outliers were excluded to prevent biased estimates by removing single

_{j}*ψ*(

_{i}*z*) values exceeding the average plus or minus three times the standard deviation

_{j}*σ*of the total number of values within each average.

The error analysis was carried out using the NAB08 radiometric and IOP profile data for the *M _{k}* validation sets complementing the

*M*training sets used to determine fitting factor

_{k}*F*. We summarize estimation of the errors in three ways for an implicitly denoted

*k*:

- Average absolute error $\overline{\left|\psi \left({z}_{j}\right)\right|}$ of all
*M*validation set combinations at each depth_{k}*z*in order to explore our assumption of asymptotic radiance._{j} - Depth-averaged absolute error $\overline{\left|{\psi}_{i}\right|}$ and relative error $\overline{{\psi}_{i}}$ for each validation set
*i*$$\overline{\left|{\psi}_{i}\right|}=\frac{1}{8}\sum _{j=1}^{8}\left|{\psi}_{i}\left({z}_{j}\right)\right|,$$with which to explore the variability of the*M*different validation set combinations, assuming a vertically homogeneous regime._{k}

For the two synthetic data sets, *ψ*̄ and
$\overline{\left|\psi \right|}$ were computed as averages of individual errors *ψ _{i}*(

*z*) and |

_{j}*ψ*(

_{i}*z*)| in IOPs estimated from the 63 simulated data points.

_{j}## 6. Results

#### 6.1. Synthetic data sets

The results for the first synthetic data set constructed using *β*̃_{FF} are shown in Table 2. For each value of chlorophyll concentration *Chl*, IOPs at 490 nm were generated. The parameter *F* was fitted for values of *b _{b}*/

*a*using all the simulated radiometric and IOP data and was essentially constant, varying from 0.074 to 0.076 while

*b*/

_{b}*a*varied from 0.037 to 0.11. Using the simulated radiometric data to test our

*b*/

_{b}*a*algorithm, the average relative and absolute errors were found to be small, demonstrating that once the fitted parameter

*F*was determined, recovery of the original data is possible to within ∼4%. This implies that the variation of

*β*̃

_{FF}due to changes in the backscattering ratio

*b*̃

*was insignificant with respect to the phase function we assumed.*

_{bp}When we replaced the validation set by repeating the simulation except for a change from *β*̃_{FF} to *β*̃_{P} for all values of *Chl*, there was a greater error in the estimates of *b _{b}*/

*a*(right half of Table 2). As one might expect, the smallest errors corresponded to

*Chl*= 0.10 to 1.0 mg m

^{−3}where

*b*̃

*was closest in both data sets. This is because*

_{bp}*F*has been “tuned” to the set of Fournier-Forand scattering phase functions used in the first data set. The Petzold data set represents a different bio-optical regime determined by a different scattering phase function; we should not expect

*F*tuned to one bio-optical regime to do as well in another.

#### 6.2. NAB08 data set

For the field data at 488 nm, Fig. 2A shows the evolution of the fitted parameter *F* as the number of training profiles *k* increases. The slope and intercept of *F*(488) show convergence towards stable mean values, as one would expect as more of the data set is used to determine *F*; more of the natural variability of *b _{b}*/

*a*is captured as additional profiles cover the experimental area and duration of the bloom. A linear model of

*F*(

*b*/

_{b}*a*) as in Eq. (49), using data for all

*k*= 6 training set combinations at each wavelength, demonstrates the wavelength dependence of

*F*beyond 488 nm (Fig. 2B, Table 3).

Figure 3 shows the evolution of the estimates of *b _{b}*/

*a*and average retrieval errors at 488 nm for

*k*= 1 ...7. As the number of training profiles increases, estimated values of

*b*/

_{b}*a*converge toward their measured values (Fig. 3A, Table 4) with a corresponding decrease in average absolute error $\overline{\left|\psi \right|}$ (Fig. 3B, Table 4). For

*b*/

_{b}*a*, the average retrieval error $\overline{\left|\psi \right|}$ for

*k*≥ 4 is approximately ±20%, with some retrieval errors $\overline{\left|{\psi}_{i}\right|}$ as high as 100%. The range of relative error $\overline{{\psi}_{i}}$ for each validation set combination

*i*remains large across all values of

*k*, but the relative errors show little bias when averaged as

*ψ*̄ (Fig. 3C, Table 5). As expected, we can conclude that taking multiple measurements within a vertically homogeneous water column will reduce the average absolute error in the IOP estimates.

The average absolute error at each depth
$\overline{\left|\psi \left({z}_{j}\right)\right|}$ is shown in Fig. 3D and decreases markedly after 25 m; this depth represents 3–5 optical depths for most profiles. This indicates that the radiance from 0–25 m is not necessarily in an asymptotic regime and a solution to the radiative transfer equation governed by the dominant eigenvalue *ν*_{0} (Eq. 25) is therefore approximate.

As the asymptotic regime is reached and its concomitant assumptions are satisfied, the average absolute error is less than 10% for *k* ≥ 4. While these field results indicate that measurable radiance and irradiance data are available within the asymptotic regime, previous studies indicate that asymptotic optical depths vary with IOPs and the angular distribution of the surface illumination [3, 12, 13]. Therefore, a significant contributer to the variability in average absolute error
$\overline{\left|\psi \right|}$ can be the failure to achieve an asymptotic light field in some surface illumination conditions; this is a limitation of our approach.

The estimates of absorption at 488 nm and average retrieval errors are shown in Fig. 4 for two different algorithms: (1) absorption *a*_{Kd}(488) estimated from *F*, *b _{b}*/

*a*,

*ν*

_{0}, and

*K*(

_{d}*z*) using Eq. (41) (Fig. 4A–C) and (2) absorption

_{as}*a*(488) estimated from

_{est}*b*/

_{b}*a*and measurements of

*b*(470) and

_{b}*b*(700) using Eq. (44) and the bio-optical model to estimate

_{b}*b*(488) (Fig. 4D–E). As we saw with estimates of

_{b}*b*/

_{b}*a*, both absorption estimates improve as

*k*increases, based on improved

*F*and

*b*/

_{b}*a*estimates. The absorption estimate

*a*

_{Kd}(488) (Fig. 4A–C) requires the estimation of the asymptotic quantity

*K*(

_{d}*z*) from

_{as}*E*data; estimates of

_{d}*K*(

_{d}*z*) from radiometric profile samples below 25 m assured that we were close to the asymptotic regime, but in some cases IOPs changed below 30 m that challenge this assumption. Applying standard error propagation techniques, we can expect a 20% uncertainty in

_{as}*b*/

_{b}*a*and similar 20% uncertainty in

*K*(

_{d}*z*) to result in ∼30% uncertainty in

_{as}*a*

_{Kd}; this is close to the average absolute error $\overline{\left|\psi \right|}$ shown in Fig. 4B for

*k*≥ 4. The use of measured values of backscattering improve the estimation of

*a*(488) to a lower average absolute error $\overline{\left|\psi \right|}$ of 15% for

*k*≥ 4. The improvements are most evident at lower values of absorption.

Finally, estimates of backscattering at 488 nm are shown in Fig. 5. While the average absolute error is ∼20% (Fig. 5B), large individual errors can be seen in groups of matchups (Fig. 5A). Large overestimates in *b _{b}*(488) are the result of overestimates in both

*b*/

_{b}*a*and

*a*

_{Kd}; these cases are found in profiles taken near the bloom peak where both higher values of absorption and backscattering were common.

The estimates of average error for any given validation set
$\overline{\left|{\psi}_{i}\right|}$ are roughly twice that of the average for all validation sets
$\overline{\left|\psi \right|}$. For example, the estimate of absorption *a*(488) has an average uncertainty of
$\overline{\left|\psi \right|}=15-16\%$ for *k* ≥ 4 (Table 5), with uncertainty for individual validation sets
$\overline{\left|{\psi}_{i}\right|}$ exhibiting as much as 36% error (Fig. 4B). This information can be used to create a subspace of absorption values ranging from 0.64*ā* to 1.36*ā*, where *ā* represents the water column average absorption estimate. This water column average could be used as a starting point and the subspace values as bound constraints for a more sophisticated multi-parameter nonlinear inversion algorithm [11].

The validation results for other wavelengths below 555 nm are consistent with those at 488 nm (Tables 4, 5, 6): (1) little or no average bias *ψ*̄ is evident in relative errors; (2) the average absolute error
$\overline{\left|\psi \right|}$ in IOP estimates decreases as *k* increases, i.e., as more training profiles are used to estimate the fitted parameter *F*; (3) while parameterization of *F* as a function of *b _{b}*/

*a*(Eq. 49) varies with wavelength (Fig. 2B), there is no discernible wavelength trend in IOP errors for wavelengths less than 555 nm; and (4) retrievals of absorption

*a*nearly always improve when measurements of

*b*are available, even when a bio-optical model is used to estimate

_{b}*b*(

_{b}*λ*) at wavelengths other than the measurement wavelengths (470 and 700 nm) (Table 5).

The 555 nm retrieval of *b _{b}*/

*a*is poor even when most of the profiles are used as training data (Table 4). This is due to Raman scattering which our algorithm does not incorporate. The slope of fitted parameter

*F*(Eq. 49) does vary with wavelength (Fig. 2, Table 3). At wavelengths below 555 nm, we suspect

*F*can help compensate for some Raman scattering at these wavelengths, while at 555 nm our assumption of no inelastic scattering in the formulation of the radiative transfer equation (Eq. 1) breaks down.

## 7. Discussion

We have developed a simple analytical algorithm that, when tuned using contemporaneous IOP and radiometric measurements, produces an IOP subspace defined by a mean IOP estimate and its associated uncertainty. Following the suggestion of Defoin-Platel and Chami [34], the resulting IOP subspace can be used as a divide and conquer strategy to reduce or eliminate ambiguities of IOPs estimated from a wider set of radiometric measurements in the same waters where IOP measurements are limited or not available. For example, the predicted IOP subspace can be used to constrain a more sophisticated inversion algorithm that incorporates, for example, models of Raman scattering and chlorophyll fluorescence and/or more accurate scattering phase functions as well as computationally efficient radiative transfer models (e.g., Ecolight-S [35]) that do not require the assumptions of an asymptotic light regime employed in the algorithms presented here.

Our single-parameter *b _{b}*/

*a*algorithm converges quickly and reliably; no Monte Carlo simulations [36] or iterative solutions with forward problem solvers [3,37] are required. The algorithm works best when light field measurements are made in the asymptotic regime and the parameter

*F*is tuned to the local environment (Fig. 3). With simulations, we have shown that the

*b*/

_{b}*a*algorithm’s performance degrades as the scattering phase function of the radiometric measurements deviates from that used to generate the data (Table 2). Consistent with these results, we found that applying the parameterization of

*F*developed for our synthetic data produced very poor estimates when applied to the NAB08 field data set. For the NAB08 calibration campaign, numerous radiometric and IOP profiles were available that adequately spanned the variability of waters. Estimating

*F*directly from this calibration campaign data (Fig. 2) effectively tunes our simple scattering phase function

*β*̃ (Eq. 7) to the North Atlantic waters under study. While we did not develop an algorithm to directly estimate the number of profiles

*k*required for such a calibration campaign, the results presented here using the bio-optical data set (eight calibration profiles) from the NAB08 process cruise can be used as a guide.

If simple IOP measurements (e.g., *b _{b}*(470) and

*b*(700)) can be made at the time of radiometric measurements, we can improve the mean absorption estimate and decrease its overall uncertainty, further reducing potential ambiguity and increasing the performance of more computationally intensive IOP estimation approaches that benefit from reduced bound constraints. This use of ancillary IOP data is an example of an “enrichment strategy” suggested by Defoin-Platel and Chami [34].

_{b}The results presented here do not account for instrument error. If there is an expected uncertainty of ±17% in backscattering [38] and ±17% in absorption, as mentioned previously, this results in a relative instrumental error of 34% in *b _{b}*/

*a*, which is within the range of the average absolute errors shown in Table 4. Similarly, the resulting estimates for absorption and backscattering (Tables 5 and 6, respectively) are also comparable to the instrument error.

One might wonder, based on the assumptions introduced and the simplicity of the algorithm for *b _{b}*/

*a*, why it works as well as it does. One interpretation is that the algorithm is a “bilinear inverse method” that weights the measured data for

*r*(

_{rs}*z*) by the “importance” of the photons as given by the adjoint radiance

_{m}*L*

^{†}(

*z*,

_{m}*μ*) ≡

*L*(

*z*, −

_{m}*μ*) (e.g., Eq. 24) of the radiative transfer equation [39, 40]. Radiance, considered from the point of view of forward radiative transfer methods, conceptually traces every photon from the source to detector. Adjoint radiance can be visualized as tracing pseudo photons received by the detector backwards to the source [41]; these pseudo photons provide insight to the importance of “how much a given position-direction matters for a given radiometric observation” [42].

As the ratio of two inherent optical properties (IOPs), *b _{b}*/

*a*is itself an inherent optical property. Particulate backscattering

*b*has been shown to correlate to particulate organic carbon [38, 43, 44]; several bio-optical models relate chlorophyll concentration to particulate absorption [45, 46]. The ratio

_{bp}*b*/

_{b}*a*might, therefore, be close to a proxy for a carbon-to-chlorophyll ratio [47], and thus be a diagnostic IOP for plankton community structure. Given the ease of estimating

*b*/

_{b}*a*from in-water measurements of

*L*and

_{u}*E*, this IOP deserves further investigation as a bio-optical proxy in the open ocean.

_{d}## Acknowledgments

This work was supported by the National Science Foundation grants OCE0628379 and OCE0628107, and NASA grants NNX08AL92G and NNX-10AP29H. This work would not have been possible without the help of Emily Kallin, who assisted with the collection of all of the NAB08 bio-optical profile data used here and Mary Jane Perry and Emmanuel Boss who loaned their bio-optical and Satlantic radiometric profilers. We wish to thank Eric D’Asaro for statistical insight, Curtis Mobley for guidance with radiative transfer simulations, and the crew and technicians of the R.V. Knorr.

## References and links

**1. **H. R. Gordon, “Inverse methods in hydrologic optics,” Oceanologia **44**, 9–58 (2002).

**2. **N. J. McCormick, “Inverse radiative transfer problems: a review,” Nucl. Sci. Eng. **112**, 185–198 (1992).

**3. **C. D. Mobley, *Light and Water: Radiative Transfer in Natural Waters* (Academic, 1994).

**4. **H. R. Gordon, O. B. Brown, and M. M. Jacobs, “Computed relationships between the inherent and apparent optical properties of a flat homogeneous ocean,” Appl. Opt. **14**, 417–427 (1975). [CrossRef] [PubMed]

**5. **H. Gordon, “Simple calculation of the diffuse reflectance of the ocean,” Appl. Opt. **12**, 2803–2804 (1973). [CrossRef] [PubMed]

**6. **N. J. McCormick, “Analytical transport theory applications in optical oceanography,” Ann. Nucl. Energy **23**, 381–395 (1996). [CrossRef]

**7. **H. Loisel and D. Stramski, “Estimation of inherent optical properties of natural waters from the irradiance attenuation coefficient and reflectance in the presence of Raman scattering,” Appl. Opt. **39**, 3001–3011 (2000). [CrossRef]

**8. **M. Stramska, D. Stramski, B. G. Mitchell, and C. D. Mobley, “Estimation of the absorption and backscattering coefficients from in-water radiometric measurements,” Limnol. Oceanogr. **45**, 628–641 (2000). [CrossRef]

**9. **H. R. Gordon and G. C. Boynton, “A radiance-irradiance inversion algorithm for estimating the absorption and backscattering coefficients of natural waters: homogeneous waters,” Appl. Opt. **36**, 2636–2641 (1997). [CrossRef] [PubMed]

**10. **G. C. Boynton and H. R. Gordon, “Irradiance inversion algorithm for estimating the absorption and backscattering coefficients of natural waters: Raman-scattering effects,” Appl. Opt. **39**, 3012–3022 (2000). [CrossRef]

**11. **E. Rehm, “Inverting light with constraints,” Poster session presented at Ocean Optics XIX, Barga, Italy (2008).

**12. **R. A. Leathers, C. S. Roesler, and N. J. McCormick, “Ocean inherent optical property determination from in-water light field measurements,” Appl. Opt. **38**, 5096–5103 (1999). [CrossRef]

**13. **B. D. Piening and N. J. McCormick, “Asymptotic optical depths in source-free ocean waters,” Appl. Opt. **42**, 5382–5387 (2003). [CrossRef] [PubMed]

**14. **S. Maritorena, D. A. Siegel, and A. R. Peterson, “Optimization of a semianalytical ocean color model for global scale applications,” Appl. Opt. **41**, 2705–2714 (2002). [CrossRef] [PubMed]

**15. **IOCCG, “Remote sensing of inherent optical properties: fundamentals, tests of algorithms, and applications,” *Reports of the International Ocean-Colour Coordinating Group*, 5, Z.-P. Lee, ed. (IOCCG, Dartmouth, Canada, 2006).

**16. **T. G. Peacock, K. L. Carder, C. O. Davis, and R. G. Steward, “Effects of fluorescence and water Raman scattering on models of remote sensing reflectance,” in *Ocean Optics X*, R. W. Spinrad, ed., Proc. SPIE 1302, 303–319 (1990).

**17. **H. R. Gordon, “Contribution of Raman scattering to water-leaving radiance: a reexamination,” Appl. Opt. **38**, 3166–3174 (1999). [CrossRef]

**18. **A. Morel and L. Prieur, “Analysis of variations in ocean color,” Limnol. Oceanogr. **22**, 709–722 (1977). [CrossRef]

**19. **J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci. **33**, 2452–2459 (1976). [CrossRef]

**20. **C. D. Mobley, L. K. Sundman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. **41**, 1035–1050 (2002). [CrossRef] [PubMed]

**21. **C. D. Mobley and L. K. Sundman, *Hydrolight 5.0, Ecolight 5.0 Technical Documentation* (Sequoia Scientific, Inc., 2008).

**22. **N. J. McCormick, “Analytic inverse radiative transfer equations for atmospheric and hydrologic optics,” J. Opt. Soc. Am. A **21**, 1009–1017 (2004). [CrossRef]

**23. **K. M. Case and P. F. Zweifel, *Linear Transport Theory* (Addison-Wesley, 1967).

**24. **N. J. McCormick and I. Kuščer, “Singular eigenfunction expansions in neutron transport theory,” in *Advances in Nuclear Science and Technology*7, E. J. Henley and J. Lewins, eds. (Academic, 1973), pp. 181–282.

**25. **N. J. McCormick, “Asymptotic optical attenuation,” Limnol. Oceanogr. **37**, 1570–1578 (1992). [CrossRef]

**26. **J. R. V. Zaneveld, “An asymptotic closure theory of irradiance in the sea and its inversion to obtain the vertical structure of inherent optical properties,” Limnol. Oceanogr. **34**, 1442–1452 (1989). [CrossRef]

**27. **G. R. Fournier and J. L. Forand, “Analytic phase function for ocean water,” in *Ocean Optics XII*, J. S. Jaffe (ed), Proc. SPIE 2258, 194–201 (1994).

**28. **T. J. Petzold, “Volume scattering functions for selected ocean waters,” *Tech. Rep. SIO 72–78* (Scripps Institution of Oceanography, La Jolla, Calif., 1972), pp. 1–78.

**29. **E. A. D’Asaro, C. Lee, M. Perry, K. Fennel, E. Rehm, A. Gray, N. Briggs, and K. Gudmundsson, “The 2008 North Atlantic bloom experiment I: overview and strategy,” Eos, Transactions, American Geophysical Union89(53), Fall Meeting Supplement, Abstract OS24A-08 (2008), http://adsabs.harvard.edu/abs/2008AGUFMOS24A..08D.

**30. **W. S. Pegau, D. Gray, and J. R. V. Zaneveld, “Absorption and attenuation of visible and near-infrared light in water: dependence on temperature and salinity,” Appl. Opt. **36**, 6035–6046 (1997). [CrossRef] [PubMed]

**31. **J. R. V. Zaneveld, J. C. Kitchen, and C. C. Moore, “Scattering error correction of reflecting tube absorption meter,” in *Ocean Optics XII*, J. S. Jaffe (ed), Proc. SPIE 2258, 44–55 (1994).

**32. **B. Efron and G. Gong, “A leisurely look at the bootstrap, the jackknife, and cross-validation,” Am. Stat. **37**, 36–48 (1983). [CrossRef]

**33. ***MATLAB® Optimization Toolbox™ 5 User’s Guide*, (The MathWorks Inc., 2010). [PubMed]

**34. **M. Defoin-Platel and M. Chami, “How ambiguous is the inverse problem of ocean color in coastal waters?,” J. Geophys. Res. **112**, C03004 (2007). [CrossRef]

**35. **C. D. Mobley, “Fast light calculations for ocean ecosystem and inverse models,” Opt. Express, submitted (2011). [CrossRef] [PubMed]

**36. **H. R. Gordon and O. B. Brown, “irradiance reflectivity of a flat ocean as a function of its optical properties,” Appl. Opt. **12**, 1549–1551 (1973). [CrossRef] [PubMed]

**37. **G. E. Thomas and K. Stamnes, *Radiative Transfer in the Atmosphere and Ocean* (Cambridge University Press, 1999). [CrossRef]

**38. **G. DallOlmo, T. K. Westberry, M. J. Behrenfeld, E. Boss, and W. H. Slade, “Significant contribution of large particles to optical backscattering in the open ocean,” Biogeosciences **6**, 947–967 (2009). [CrossRef]

**39. **J. J. Duderstadt and W. R. Martin, *Transport Theory* (Wiley, 1979).

**40. **G. I. Bell and S. Glasstone, *Nuclear Reactor Theory* (Van Nostrand-Reinhold, 1970).

**41. **Q. Min and L. C. Harrison , “An adjoint formulation of the radiative transfer method,” J. Geophys. Res. **101**, 1635–1640 (1996). [CrossRef]

**42. **A. B. Davis and Y. Knyazikhin, “A Primer in Three-Dimensional Radiative Transfer,” in *Three-Dimensional Radiative Transfer in the Cloudy Atmosphere*, A. Davis and A. Marshak, eds. (Springer-Verlag, 2005), pp. 153–242. [CrossRef]

**43. **H. Loisel, J. Nicolas, P. Deschamps, and R. Frouin, “Seasonal and inter-annual variability of particulate organic matter in the global ocean,” Geophys. Res. Lett. **29**(24), 2196 (2002). [CrossRef]

**44. **D. Stramski, R. A. Reynolds, M. Babin, S. Kaczmarek, M. R. Lewis, R. Röttgers, A. Sciandra, M. Stramska, M. S. Twardowski, B. A. Franz, and H. Claustre, “Relationships between the surface concentration of particulate organic carbon and optical properties in the eastern South Pacific and eastern Atlantic Oceans,” Biogeosciences **5**, 171–201 (2008). [CrossRef]

**45. **L. Prieur and S. Sathyendranath, “An optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, dissolved organic matter, and other particulate materials,” Limnol. Oceanogr. **26**(4), 671–689 (1981). [CrossRef]

**46. **A. Morel, “Light and marine photosynthesis: a spectral model with geochemical and climatological implications,” Prog. Oceanogr. **26**, 263–306 (1991). [CrossRef]

**47. **M. J. Behrenfeld and E. Boss, “The beam attenuation to chlorophyll ratio: an optical index of phytoplankton physiology in the surface ocean?,” Deep-Sea Res I **50**, 1537–1549 (2003).