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_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 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_b. Alternatively, a can be estimated from b_b/a, measurements of b_b at one or more wavelengths, and a model for the spectral behavior of 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 β̃(μ₀) = β̃(μ′,φ′,μ,φ) for μ₀ the cosine of the angle between {μ′,φ′} and {μ,φ}, the radiative transfer equation for the semi-infinite medium is [3]

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_m, −1) at measurement depth z_m and the downward planar irradiance E_d(z_m),

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_m 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 z_m 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,

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,

It 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_m) and 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_b, which is consistent with the quasi-single-scattering approximation [3–5].

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

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

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_m) and Λ(z_m), from Eq. (21) the final form for the 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 ν₀,

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

To estimate μ̄(z_as) in order to obtain the absorption coefficient from Eq. (37), it follows from Eqs. (26), (28), and (31) that

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₁ and z₂ and to use

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_b from

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_d 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 Chl (0.01, 0.03, 0.1, 0.3, 1, 3, 10 mg m⁻³) 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_p, which resulted in a distinct scattering phase function selection for each chlorophyll concentration. To assess the ability of our algorithm to retrieve 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⁻³) during the bloom and mixed ciliate/picoeukaryote community (Chl < 1 mg m⁻³) 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⁻¹, 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_bp using a power-law b_bp(λ) = b_bp(λ₀)(λ₀/λ)^γ 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.

 

Fig. 1 An IOP and a radiometric profile from the NAB08 dataset. (A) Typical IOPs were nearly homogeneous in the lit portion of the water column. (B) Radiometric profiles of L_u and E_d at 488 nm show little curvature until 40 m, another indication of vertically homogeneous water.

Download Full Size | PDF

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_b(470) and b_b(700) as described above and b_b/a is computed from a(488) and the estimated b_b(488).

When processing the measured data, no attempt was made to incorporate possible instrumentation inaccuracies. Radiometric measurements were discarded for L_u(z) < 10⁻⁴ μW cm⁻² sr⁻¹ or E_d(z) < 10⁻² μW cm⁻².

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_as). 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 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.

Table 1. Example training and validation sets of NAB08 calibration profiles for k = 1,2 used for cross validation.

View Table | View all tables in this article

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_b, and 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_F, where

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_d computed from simulated or field radiometric data were used to estimate b_b/a. The contemporaneous measured or simulated values of b_b and 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

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

The error analysis was carried out using the NAB08 radiometric and IOP profile data for the M_k validation sets complementing the M_k training sets used to determine fitting factor F. We summarize estimation of the errors in three ways for an implicitly denoted k:

For the two synthetic data sets, ψ̄ and $\overline{\left|\psi \right|}$ were computed as averages of individual errors ψ_i(z_j) and |ψ_i(z_j)| in IOPs estimated from the 63 simulated data points.

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̃_bp was insignificant with respect to the phase function we assumed.

Table 2. Average relative error ψ̄ and absolute error $\overline{\left|\psi \right|}$ (in percent) for estimates of b_b/a from radiative transfer simulations at 490 nm. Chl in mg m⁻³. Parenthetical values indicate the range of $\overline{\left|\psi \right|}$ over all depths.

View Table | View all tables in this article

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⁻³ where b̃_bp was closest in both data sets. This is because 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).

 

Fig. 2 Estimates of F using r_rs = L_u/E_d and IOP measurements from the NAB08 dataset. A) Values of fitted parameter F at 488 nm for all training set combinations (black x’s) and the linear fit for each value of k. Color indicates number of profiles k used in the training set. (B) fitted parameter F at selected wavelengths (color), created by linear fit for all k = 6 training set combinations.

Download Full Size | PDF

Table 3. Slope, intercept and R² for fitted parameter F = p₁(b_b/a) + p₂ for AC-9 wavelengths ≤ 532 nm with a linear least squares fit using all k = 6 training set combinations.

View Table | View all tables in this article

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.

 

Fig. 3 (A) Matchups of estimated and measured b_b/a values show convergence towards 1:1 line (black dashed line) as k increases; each colored dot represents a single matchup. (B) Depth-averaged absolute error for each validation set ( $\overline{\left|{\psi }_i\right|}$, colored dots) remains variable irrespective of the number of profiles k used to determine F. The overall average absolute error $\overline{\left|\psi \right|}$ for the algorithm (black dashed line) decreases most from k = 1 to 4. (C) Similarly, relative error depth-averaged for each validation set ( $\overline{{\psi }_i}$, colored dots) shows variability but little average bias ψ̄ (black dashed line). (D) Absolute error averaged across all validation set combinations at each depth $\overline{\left|\psi \left(z_j\right)\right|}$ decreases markedly in deep water.

Download Full Size | PDF

Table 4. Relative and absolute error (in percent) for estimates of b_b/a using different numbers of profiles k in the training set for fitted parameter F

View Table | View all tables in this article

Table 5. Relative and absolute error (in percent) for estimates of a using Eqs. (41) and (44) using different numbers of profiles k in the training set for fitted parameter F

View Table | View all tables in this article

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 ν₀ (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, ν₀, and K_d(z_as) using Eq. (41) (Fig. 4A–C) and (2) absorption a_est (488) estimated from b_b/a and measurements of b_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/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_as) from E_d data; estimates of K_d (z_as) 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 b_b/a and similar 20% uncertainty in K_d (z_as) to result in ∼30% uncertainty in 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.

 

Fig. 4 Values of absorption a (A) and average error (B–C) at 488 nm, estimated from F, b_b/a, ν₀, and K_d (z_as), using Eq. (41), and (D–F) estimated from b_b/a and measurements of b_b(470) and b_b(700) using Eq. (44). Colors indicate the number of profiles k used to determine the value of F; remaining profiles are used to estimate b_b/a, K_d (z_as) and a.

Download Full Size | PDF

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.

 

Fig. 5 Values of backscattering b_b (A) and average error (B) and (C) at 488 nm, estimated from b_b/a, F, and K_d (z_as). Colors indicate the number of profiles k used to determine the value of F; remaining profiles are used to estimate b_b/a, K_d (z_as) and a.

Download Full Size | PDF

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_b are available, even when a bio-optical model is used to estimate b_b(λ) at wavelengths other than the measurement wavelengths (470 and 700 nm) (Table 5).

Table 6. Relative and absolute error (in percent) for estimates of b_b using different numbers of profiles k in the training set for fitted parameter F

View Table | View all tables in this article

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_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].

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_m) by the “importance” of the photons as given by the adjoint radiance 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_bp 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 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_u and E_d, this IOP deserves further investigation as a bio-optical proxy in the open ocean.