Optical inversions of the water column based on glider measurements

Catherine Mitchell; Howard R. Gordon; Bruce Bowler; Dave Drapeau; William M. Balch

doi:10.1364/OE.26.032824

1. Introduction

One of the major challenges in ocean optics is the so-called “inverse problem”: the estimation of the inherent optical properties (IOPs) from the underwater visible light field [1–3]. IOPs depend on the composition and relative concentration of any optically significant materials (OSMs) present in the water column, and consequently can be used as a starting point for many studies including the estimation of (1) individual concentrations of any OSMs such as phytoplankton, suspended mineral particles or colored dissolved organic matter (CDOM) [4–6], (2) primary productivity [7,8], (3) particulate carbon concentrations [9], (4) depth of penetration of solar radiation [10] and (5) particle dynamics in shelf seas [11,12].

There have been several in-water studies that use radiative transfer modelling to develop inversion methods for estimating IOPs from radiometric quantities and/or apparent optical properties (AOPs), and evaluate these methods using profiling instruments [13–19]. The studies vary in the inversion method, from a relatively simple set of equations based on radiative transfer theory [19], to numerically solving the full radiative transfer equation [18,20]. All these studies have shown reasonable recoveries of IOPs in different regions across the globe; however they all rely on multiple instruments, often with separate platforms for the radiometric and IOP instruments as opposed to a depth profile of AOPs and IOPs on a single platform.

Recent advances in sensor technology have enabled autonomous underwater vehicles (AUVs) to become ocean observing platforms equipped with a suite of different oceanographic instruments [21]. Consequently, AUVs have been used to study the full water column in detail [22–25], but only a few of these studies exploit the potential of optical sensors [26–28], and even fewer use AUVS for inversion methods.

The first study to mount a radiometer on an AUV demonstrated three different approaches to retrieve the absorption coefficients for phytoplankton and nonalgal colored matter [29]. While Brown et al. [29] had success quantifying and mapping bio-optical properties along the coastline of New Jersey, Moline et al. [30] applied the same method to the Santa Barbara Channel with marginal success, and in Hawaiian waters unsuccessfully. Moline et al. [30] attributed these results to the assumption of a homogenous uniform light field in the inversion model developed by Brown et al. [29], whereas their sites had variable spectral downwelling irradiance due to broken cloud cover and wave focusing effects. In both these cases, there were no other sensors on the AUVs to make a direct comparison with the bio-optical properties which were estimated.

Since these pioneering studies, AUVs have been developed which are equipped with sensors to measure the spectral downwelling planar irradiance, $E_{d} (λ)$ ( $μ$ W cm $^{- 2}$ nm $^{- 1}$ ), and the spectral upwelling radiance, $L_{u} (λ)$ ( $μ$ W cm $^{- 2}$ sr $^{- 1}$ nm $^{- 1}$ ), of the light field, in addition to the backscattering coefficient, $b_{b} (λ)$ (m $^{- 1}$ ), at one or more wavelengths. Therefore, these so-called bio-optical underwater gliders provide a unique platform for developing and evaluating inversion methods of the euphotic zone.

In this study, we demonstrate how gliders can be used to estimate IOPs of the water column from glider-measured $E_{d} (λ)$ and $L_{u} (λ)$ profiles. The inversion method used here is a slightly modified version of Gordon and Boynton [31], which deals with vertically stratified waters and was originally verified on simulated data and experimentally on hyperspectral field data [32]. We will give an overview of the modified inversion procedure and present a protocol to quality control (QC) glider-measured radiometry profiles. Finally, we will assess the performance of (1) the inversion method on modeled data and (2) the combined QC-inversion procedure on glider data collected within the Gulf of Maine.

2. Data and methods

2.1 Inversion method

The inversion method iteratively estimates IOPs from full profile measurements of $E_{d} (λ_{0})$ and $L_{u} (λ_{0})$ for a given wavelength, $λ_{0}$ , and requires no extra information from $E_{d} (λ)$ or $L_{u} (λ)$ at other wavelengths. The steps of the inversion are given below, with modifications to the method of Gordon and Boynton [31] highlighted. Note in the following description the superscript $m$ indicates a measured quantity and $(i)$ (where $i$ is an integer) indicates the $i$ -th iteration.

1. Estimate the measured upwelling irradiance profile, $E_{u}^{m} (z)$ : assuming first that $L_{u}$ is totally diffuse such that $Q = E_{u} / L_{u} = π$ , hence $E_{u}^{m} (z) = π L_{u}^{m} (z)$ , where $L_{u}^{m} (z)$ is the measured upwelling radiance.
2. Estimate the absorption coefficient (for the 0-th iteration), $a^{(0)} (z)$ (m $^{- 1}$ ): $a^{(0)} (z) = μ_{0} K_{E}^{m} (z) = \cos (θ_{0 w}) [- \frac{d}{d z} \ln (E_{d}^{m} - E_{u}^{m})]$
where $K_{E}^{m}$ is the diffuse attenuation coefficient for the vector irradiance, $E_{d}^{m}$ is the measured downwelling irradiance, and $μ_{0}$ is the cosine of $θ_{0 w}$ , the solar zenith angle measured beneath the surface. Note the solar zenith angle above the surface is calculated based on position and time for each measurement.
3. Estimate the backscattering coefficient, $b_{b}^{(0)} (z)$ : use the following, $\begin{array}{l} b_{b}^{(0)} (z) & = & a^{(0)} (z) X^{m} (z) \\ = & 3 a^{(0)} (z) {R^{m} (z) - \frac{d R^{m} (z)}{d z} \int_{z}^{z_{m a x}} {[\frac{E_{d}^{m} (z^{'})}{E_{d}^{m} (z)}]}^{2} d z^{'}} \end{array}$
where $R^{m} (z) = E_{u}^{m} (z) / E_{d}^{m} (z)$ , which can be approximated in vertically stratified waters by $R = 〈 X 〉 / 3$ for $X = b_{b} / a$ . The approximation, $R = 〈 X 〉 / 3$ , can be solved for X, leading to the above approximation for $b_{b}^{(0)} (z)$ (see Gordon and Boynton [31] for details).
4. Assume a scattering phase function. At this point, Gordon and Boynton [31] calculates the scattering coefficient, $b (λ)$ (m $^{- 1}$ ), based on $b_{b}^{(0)}$ from step 3 and the assumed phase function. But the radiative transfer code we use (see step 5 for details) only requires the backscattering coefficient and the scattering phase function (which we assume to be a Fournier-Forand phase function based on a backscatter fraction [33]). Hence, we do not explicitly calculate the scattering coefficient as it is done automatically by the selection of the phase function.
5. Use the IOP estimates in the radiative transfer equation to calculate the corresponding $L_{u} (z)$ , $E_{d} (z)$ , $E_{u} (z)$ and $E_{0} (z)$ values (where $E_{0} (z)$ is the scalar irradiance). We use EcoLight-S (Sequoia Scientific [34],) for the radiative transfer calculations. It uses Fournier-Forand scattering phase functions based on the backscattering fraction, $\tilde{b} = b_{b} / b$ , which is provided by the user based on regional knowledge. Gordon and Boynton [31] used Monte Carlo radiative transfer code as presented in [35] and [36].
6. Re-calculate $R (z)$ and $Q (z)$ based on $E_{d}^{(1)} (z)$ , $E_{u}^{(1)} (z)$ and $L_{u}^{(1)} (z)$ (the results from the radiative transfer calculations).
7. Re-estimate the absorption coefficient, $a^{(1)} (z)$ : update the estimate of $E_{d}^{m} (z) (= Q (z) L_{u}^{m} (z))$ and use the mean cosine, $\bar{μ}$ , in place of $μ_{0}$ $a^{(1)} (z) = \bar{μ} (z) K_{E}^{m} (z) = \frac{E_{d}^{(1)} (z) - E_{u}^{(1)} (z)}{E_{0}^{(1)} (z)} [- \frac{d}{d z} \ln (E_{d}^{m} - E_{u}^{m})]$
8. Re-estimate the backscattering coefficient, $b_{b}^{(1)} (z)$ : first, calculate the error in $R$ , $Δ R = R^{m} - R^{(1)}$ , and the error in $X$ , $Δ X (z) = 3 {Δ R (z) - \frac{d Δ R (z)}{d z} \int_{z}^{z_{m a x}} {[\frac{E_{d}^{m} (z^{'})}{E_{d}^{m} (z)}]}^{2} d z^{'}}$
Then calculate the error in $b_{b}^{(0)} (z)$ , $Δ b_{b} (z) = Δ X (z) a^{(1)} (z)$ . Finally, estimate $b_{b}^{(1)} (z)$ using $b_{b}^{(1)} (z) = b_{b}^{(0)} (z) + f Δ b_{b} (z)$ , where $f$ is a constant governing the speed of convergence and equals 0.2.
9. Repeat steps 4 - 8 until the residual error, $δ^{(n)}$ , after $n$ iterations reaches a minimum, or until 30 iterations have been reached: $δ^{n} = \frac{1}{2 N} \sum_{i = 1}^{N} | \ln (E_{d}^{(n)} (z_{i})) - \ln (E_{d}^{m} (z_{i})) | + \frac{1}{2 N} \sum_{i = 1}^{N} | \ln (E_{u}^{(n)} (z_{i})) - \ln (E_{u}^{m} (z_{i})) |$
A limit of 30 iterations was selected for the algorithm as, for most of our data, we found the inversion results converged after 15 - 20 iterations.
10. The inverted IOPs are the $a$ and $b_{b}$ estimates which were used to solve the radiative transfer equation on the n-th iteration. Note the resultant IOPs are total IOPs i.e. they include a component due to water.

The major updates were (1) replacing the radiative transfer calculations with EcoLight-S and (2) changing the scattering phase function to a Fournier-Forand phase function. Gordon and Boynton [31] used a Monte Carlo radiative transfer code to solve the radiative transfer equation which is computationally expensive. We reduced the computational time by using EcoLight-S, which was developed with the aim of achieving fast run times and has been shown to accurately estimate IOPs [34]. In steps 4 and 8, we replaced this formula with one based on Zaneveld’s method of treating radiative transfer in stratified media [37] to investigate if this improved the inversion; however, nearly identical retrievals were obtained.

In the present inversions we have ignored Raman scattering. For the waters we are investigating, the influence of Raman scattering is small, and decreases rapidly with depth, because the irradiance attenuation coefficient at the excitation wavelength (approximately 450 nm) is larger than that at the wavelength at which we will be performing inversions (532 nm). This will be the case in the Gulf of Maine for inversions in the blue as well. (For low-chlorophyll, clear ocean water Raman scattering does affect the retrievals at 532 nm [32] and will have to be corrected for.)

2.2 Modeled data

To test the new implementation of the inversion algorithm, we generated stratified profiles following Gordon and Boynton [31], with IOPs defined by

a (z) = a_{0} + a_{1} \exp [- \frac{{(z - z_{a})}^{2}}{2 σ_{a}^{2}}] b_{b} (z) = b_{0} + b_{1} \exp [- \frac{{(z - z_{b})}^{2}}{2 σ_{b}^{2}}]

where

z_{a}

and

z_{b}

are the depths of the maximal values of

a

and

b_{b}

respectively, and the parameters

a_{1}

,

b_{1}

,

σ_{a}

and

σ_{b}

control the magnitude of the maximal values above

a_{0}

and

b_{0}

and the depth range of the maxima. Two different IOP profiles were simulated, denoted aub2 and a2b2 (notation consistent with Gordon and Boynton [31], where the letter-number combination indicates the depth of the maximal value and “u” indicates a uniform profile e.g. aub2 indicates a uniform

a (z)

profile and a

b_{b} (z)

profile with a maximum value at 2 m). The values for the profile parameters were taken from Gordon and Boynton [31] (for a direct comparison of the inversion performance) and are given in Table 1. The water types modeled are moderately turbid, with a maximum

a

= 0.2 m

^{- 1}

and

b_{b}

= 0.014 m

^{- 1}

. Note that in neither of the two simulated profiles do absorption and scattering covary with depth. This makes the inversions more difficult. A Fournier-Forand phase function was selected with a backscattering fraction of 0.036 (in agreement with the Henyey-Greenstien function with g = 0.85 used by Gordon and Boynton [31]). In addition to varying the IOPs of the water column, the solar zenith angle,

θ_{0}

, was set to 0° and 60°, giving a total of 4 different profiles. The vertical light distributions for these profiles were generated with EcoLight (Sequoia Scientific) at 20 equally spaced depths between 0 and 4.75 m, using the “Measured IOPs” run and with all other parameters set to the default.

Table 1. Modeled IOP profile parameters

View Table | View all tables in this article

Initially, the inversion was run on these simulated profiles using the particle phase function which was used to generate the data i.e. Fournier-Forand with ${\tilde{b}}_{b} =$ 0.036. In addition, the inversion was run using a Fournier-Forand phase function with a backscattering fraction of 0.011 to compare with Gordon and Boynton [31] for cases where an incorrect phase function was used in the retrieval of the IOPs.

2.3 Field data

As part of the Gulf of Maine North Atlantic Time Series (GNATS [38],), two Slocum Gliders equipped with a suite of optical and biogeochemical sensors have been used to observe the vertical structure of the water column across the Gulf of Maine. The two gliders, called Henry and Grampus, have similar configurations: both have CTDs, a WETLabs EcoPuck which measures $b_{b} (532)$ , chlorophyll fluorescence (excitation/emission at 470/695 nm) and CDOM fluorescence (excitation/emission at 370/460 nm), and Satlantic spectral downwelling irradiance and upwelling radiance sensors for ocean color satellite validation. Grampus also has an oxygen sensor. The wavelengths of the radiometers are slightly different: Henry samples at 412 nm, 443 nm, 490 nm, 510 nm, 532 nm, 555 nm & 670 nm, and Grampus samples at 380 nm, 443 nm, 490 nm and 532 nm.

For this study, we considered the 532 nm data for all the glider missions as we could validate the backscattering coefficient estimated from the inversion method directly with the $b_{b} (532)$ sensor. In order to validate the inverted absorption coefficient, the inversion was performed at 443 nm for one of the glider missions. In order to ensure only data sampled under optimal light conditions were analyzed, we selected data which occurred within 3 hours of local apparent noon.

The only assumption needed for this inversion method is the particle phase function. Surface IOP measurements have been made within the Gulf of Maine since 1998 as part of GNATS. Therefore, we can use a Fournier-Forand phase function based on the median backscattering fraction observed within the Gulf of Maine. The backscattering coefficient at 532 nm was determined from the volume scattering function measured by Wyatt Technologies DAWN and EOS sensors. The scattering coefficient was calculated as the difference between the attenuation and absorption coefficients as measured by a WETLabs ac-9. Based on over 10,000 surface IOP measurements [39], the median backscattering fraction was calculated to be 0.0074 $\pm$ 0.0033. It is worth noting that there is a slight mismatch in the wavelengths used in this calculation, as the closest ac-9 wavelength to 532 nm is 550 nm.

2.3.1 Data processing protocol

The glider data are split up into individual profiles, or yos. The sensors are oriented at complimentary angles to place them parallel to the sea surface when the glider is ascending, thus only upward yos are used for radiometric sampling. Between the two gliders, there are a total of 22 individual missions over a period of 10 years (2008 – 2018), providing over 1,300 yos for analysis. For each individual yo, the data are binned to 1 m depths. Each yo is then quality-controlled and classified based on its overall quality with the details for these procedures given in the next two sections.

2.3.1.1 Quality control

The gliders are continually sampling under varying sky conditions hence the quality control procedure outlined below was developed to detect and reject data which have been subjected to illumination variations, mainly due to changing cloud cover during the glider ascent.

Any data with $E_{d} (532) < 0.01$ $μ$ Wcm $^{- 2}$ nm $^{- 1}$ or $L_{u} (532) < 0.0002$ $μ$ Wcm $^{- 2}$ nm $^{- 1}$ sr $^{- 1}$ were rejected as dark values. To minimize wave focusing affects, data shallower than 10 m were rejected. After these values had been removed, only yos which had more than 10 depth measurements were kept. At this point, the $E_{d} (z)$ and $L_{u} (z)$ data were log transformed. Each individual measurement and full yo were assigned a classification flag and type, respectively (see details in Sections 2.3.1.2 and 2.3.1.3).

Passing clouds can cause sharp decreases in radiometric profile data hence these artifacts were identified by checking if $\log (E_{d} (z))$ is monotonically decreasing with increasing depth. If it is not, $\log (E_{d} (z))$ was checked against $\log (E_{d} (z - n))$ , where $n$ was varied from 2 to 10 m, to see if any points in the shallower 10 m are also smaller than $\log (E_{d} (z))$ , and if so, those points are removed from both the $\log (E_{d} (z))$ and $\log (L_{u} (z))$ profiles.

The final part of the quality-control procedure was to iteratively fit polynomials for $\log (E_{d} (z))$ vs $z$ and $\log (L_{u} (z))$ vs $z$ , removing outliers at each iteration. First, a linear (first-order polynomial) was fitted to the data via least squares regression. The square of the residuals, or the deviation, $σ$ , was determined for each data point, and we rejected any points with $σ > 3 \bar{σ}$ , were $\bar{σ}$ is the mean deviation. This fitting and rejection of outliers was repeated twice more but using a third-order polynomial fit and finally a fourth-order polynomial fit. We used the data points from the final fourth-order polynomial fit as our quality-controlled $E_{d}$ and $L_{u}$ profiles.

2.3.1.2 Classification of each measurement

Following Organelli et al. [40], each measurement of every glider yo was assigned a classification flag, indicating the quality of that data point. The classification is done after rejecting the dark values (i.e. the low values) in the raw $E_{d}$ and $L_{u}$ profiles, but before the removal of cloud artifacts. We summarize the procedure below but see Organelli et al. [40] for full details.

A fourth-order polynomial is fit between $\log (E_{d})$ and $z$ via least squares. If $R^{2} < 0.996$ , every data point within the yo is flagged as 3 (or “probably bad”), and the classification procedure stops. However, for $R^{2} \geq 0.996$ , any residual ( $δ$ ) outside two times the standard deviation (SD) of the mean of the residuals ( $\bar{δ}$ ) is flagged as 3. Any points with $δ > \bar{δ} \pm 2 S D$ are removed, and the fourth-order polynomial fit is redone. As before, if $R^{2} < 0.996$ every data point within the yo is flagged as 3, and the classification procedure stops. If $0.998 > R^{2} \geq 0.996$ and $δ > \bar{δ} \pm 2 S D$ , assign flag 3 to the data point, but if $δ < \bar{δ} \pm 2 S D$ , assign flag 2 (or “marginal”) to the data point. Finally, for $R^{2} \geq 0.998$ , if $δ > \bar{δ} \pm 2 S D$ assign flag 3, if $δ > \bar{δ} \pm S D$ assign flag 2, and $δ < \bar{δ} \pm S D$ assign flag 1 (or “good”).

2.3.1.3 Classification of full yo

The full yos were classified based on the coefficient of determination of the fourth-order polynomial fits, as in Organelli et al. [40]. Any yos which had $R^{2} < 0.996$ after the first fit in Section 2.3.1.2 were assigned as Type 3 (or “probably bad”). After the second polynomial fit, any yos still with $R^{2} < 0.996$ were classified as Type 3, $0.996 \leq R^{2} < 0.998$ were classified as Type 2 (or “marginal”), and $R^{2} \geq 0.998$ were classified as Type 1 (or “good”). Note in Organelli et al. [40], different $R^{2}$ values used as the classification limits depending on the wavelength of interest. Here, we have used their values for 490 nm, as the closest wavelength to 532 nm.

2.4 Inversion performance metrics

To evaluate the performance of the inversion on the field data, the following performance metrics were used: the slope and coefficient of determination $(R^{2})$ of a least squares regression, the mean absolute error (MAE) and the bias. The absolute error is defined as $∣ y - x ∣$ , and the bias is defined as $(Σ_{i = 1}^{n} (y_{i} - x_{i})) / n$ , for $y =$ the inversion result and $x =$ the measured value. Statistical analysis was done on logarithmically transformed data, hence the resulting MAE and bias had to be back-transformed from $\log_{10}$ space. Explicitly [41],

b i a s = 10^(\frac{Σ_{i = 1}^{n} \log_{10} (y) - \log_{10} (x)}{n}), M A E = 10^(\frac{Σ_{i = 1}^{n} ∣ \log_{10} (y) - \log_{10} (x) ∣}{n})

When determined in this way, the bias and MAE are multiplicative metrics, that is, a bias of 1.5 indicates a positive bias of 50% (greater than the measured value), whereas a bias of 0.5 indicates a negative bias of 50%. Similarly for the MAE, but those values always exceed unity.

For consistency with Gordon and Boynton [31], the performance of the inversion method on the synthetic data set was evaluated on the relative error (RE) and the mean absolute relative error (MARE), where the relative error is defined as $(y - x) / x$ for $x =$ the modeled value and $y =$ the inversion result.

3. Results and discussion

3.1 Comparison with the modeled data

The inversion results for the eight different modeled data runs are shown in Fig. 1, with the MAREs presented in Table 2. Also included in Table 2 are the results of the original model by Gordon and Boynton [31] for comparison (reproduced from their Table 2 and Table 3).

Fig. 1 Relative error (%) of the estimated values of (a) $a$ and (b) $b_{b}$ for profile aub2 and of (c) $a$ and (d) $b_{b}$ for profile a2b2. Note the details of differences between aub2 and a2b2 are described in the methods and Table 1.

Download Full Size | PDF

Table 2. Mean absolute relative error (%) for the retrieved $a$ and $b_{b}$ using modeled data

View Table | View all tables in this article

The absorption coefficient was recovered with a lower MARE than the backscattering coefficient, irrespective of solar zenith angle or whether the correct phase function was used. Using the incorrect phase function in the inversion led to an increase in the MARE of approximately 3.5% for $a$ and between 2.5% and 6% for $b_{b}$ . Increasing the solar zenith angle had a minimal effect of the estimation of $a$ (MARE < 0.5% for all cases). Similarly, when the correct phase function was used, the difference in the MARE for estimating $b_{b}$ with differing solar angles was < 0.5%, however, when the incorrect phase function was used, this difference increased to approximately 3%.

The MAREs from this study were smaller than observed by Gordon and Boynton [31] when the correct phase function was used in the inversion. The only exception to this was for $b_{b}$ , for the aub2 profile. In comparison, when the incorrect phase function was used, the MAREs were smaller for Gordon and Boynton [31]. There are three components at play here: (1) the use of Fournier-Forand phase functions in place of Henyey-Greenstein, (2) replacing the Monte Carlo method for solving the RTE with an invariant imbedding method (EcoLight-S) and (3) Gordon and Boynton [31] used a coupled ocean-atmosphere radiative transfer code, i.e., the atmosphere (with molecular and aerosol scattering) was a part of each iteration of the Monte Carlo simulation, while in ELS a separate model is used to provide the radiance incident on the water surface from the sun and sky. In terms of the phase function, Mobley et al. [42] found changing the overall shape altered the underwater light field by only a few percent. The most important factor was to have the correct backscatter fraction. While it is hard to fully separate out the effects of changing the phase function and the RT model, our results are generally in agreement with Mobley et al. [42]: the difference between the MARE is <1.5% for our study vs Gordon & Boynton [31] when the correct phase function is used, but when the backscatter fraction was changed there was a larger change in the MARE.

A major advantage of the inversion method presented here is a 98% reduction in computational time. The method of Gordon and Boynton [31] (with a Monte Carlo approach to solving the radiative transfer equation) takes approximately 4 hours to invert each profile. Whereas, the EcoLight-S based approach takes approximately 5 minutes to invert each profile, which makes it a lot more feasible to invert 1,300 glider yos. This reduction in processing time does not affect the performance of the inversion; overall, we are estimating profiles of $a$ and $b_{b}$ with a MARE of less than 10%. These results are similar to the performance of Gordon and Boynton [31] and a different version of the inversion method as modified by Fan et al. [20] and tested on modeled data.

3.2 Application to the field data

The inversion results for the backscattering coefficient for all the glider missions, irrespective of how they were classified (see Section 2.3.1.3), are shown in Fig. 2. Linear least squares regression was performed on the data, with a resultant slope of 0.567 on linear scales and $R^{2} = 0.63$ (see Fig. 2a). There are two populations visible in the Fig. 2a which correspond to the two different gliders. For the logarithmically transformed data, the slope was 0.667 and $R^{2} = 0.42$ .

Fig. 2 Inversion results for all glider missions (Henry in orange, Grampus in green), shown as $b_{b}$ regressions on (a) linear scale, and (b) log scale, with the dashed line indicating the 1:1 line, the solid line the least squares regression and histograms on the top and right axes indicating the distribution of the data points; and (c) the frequency distribution of the ratio of $b_{b}$ inverted to $b_{b}$ measured, with the dashed lines indicating $\pm 25$ %.

Download Full Size | PDF

The regression results appear to show a lot of variability between the estimated and measured $b_{b}$ values (wide spread around the 1:1 line in Fig. 2a and b, and the low $R^{2}$ values). There is a total of 157,299 points within this analysis, hence the outliers look significant in the regression plots. However, if we consider the ratio between recovered $(b_{b}^{e s t})$ and measured $(b_{b}^{m e a})$ backscattering coefficients, it was found 73% of $b_{b}^{e s t}$ is within 25% of $b_{b}^{m e a}$ (Fig. 2c).

3.2.1 Depth offset analysis

When considering the full profiles, depth offsets were sometimes observed between the $b_{b}^{e s t}$ and $b_{b}^{m e a}$ profiles (see Fig. 3a). To identify how often this depth offset occurred, we first determined the similarity of the profiles by calculating the correlation coefficient between $b_{b}^{e s t}$ and $b_{b}^{m e a}$ . If the correlation coefficient was >0.75, a cross-correlation was performed between the two profiles. Note profiles with a correlation coefficient of <0.75 were categorized as too different to successfully compare. The depth offset was taken as the location of the peak value in the cross-correlation (see Fig. 3b).

Fig. 3 (a) the measured and estimated backscattering profiles for one yo, with a depth offset visible between the two profiles, (b) cross correlation between $b_{b}^{e s t}$ and $b_{b}^{m e a}$ (see text for details).

Download Full Size | PDF

For all of the glider yos, 47% had a correlation coefficient >0.75. Of those yos, 20% had a depth offset of 1m or more, which translated to, 9.5% of the total number of glider yos with a measurable depth offset of 1m or more.

We investigated the relationship between the depth offset and different variables as measured by the gliders at the top of the yo (including the time of day, $L_{u} (532)$ , $E_{d} (532)$ , pitch and roll). Unfortunately, no significant relationships were found; the depth offset occurred for both the gliders, and often neighboring yos did not have a depth offset. Despite not being able to determine the cause of the offset, we can correct for it. For the relevant yos, the offset identified by the cross-correlation analysis was removed from $b_{b}^{e s t}$ .

3.2.2. Depth-corrected $b_{b}$ inversion results

The regression analysis was re-done for the depth-corrected backscattering coefficient data for (1) all the data, (2) data with a measurement classification flag = 1, and (3) data with a full yo classification type = 1. Figure 4 shows the slope, coefficient of determination, bias and MAE for the all the glider data together and separated by individual glider missions, with statistics performed on log transformed data.

Fig. 4 Performance of the inversion method for estimating $b_{b}$ for individual glider missions. On the x-axis, each glider mission is labeled by a letter (g = Grampus, h = Henry) and a mission number, and “all” refers to all the missions together. (a) slope of a least squares regression between $b_{b}^{e s t}$ and $b_{b}^{m e a}$ , (b) $R^{2}$ for the least squares regression, (c) bias between $b_{b}^{e s t}$ and $b_{b}^{m e a}$ , and (d) MAE between $b_{b}^{e s t}$ and $b_{b}^{m e a}$ . All performance metrics were calculated on log-transformed data. See Section 2.4 for details of performance metrics and 2.3.1 for details of classification flags and types.

Download Full Size | PDF

The regression results are similar for all the missions regardless of which glider was used. For all the data, the minimum slope was 0.65 and maximum slope was 0.89. The slopes are similar within a given mission irrespective of which flag classifications or yo-types are considered. The $R^{2}$ values vary between 0.50 and 0.86 for all the data, but these increase to a mean $R^{2} = 0.74$ and $R^{2} = 0.75$ for data with flag = 1 or of yo-type = 1 respectively. The $R^{2}$ statistic is sensitive to outliers and does not reveal any biases there may be in the data. Additionally, the regression slope is sensitive to outliers and biases at a particular point in the full data range (which ultimately can produce meaningless slopes). Hence, Seegers et al. [41] suggest using the bias and MAE as algorithm performance metrics which are less sensitive to outliers and accurately reflect error magnitudes and direction.

For the glider inversion results, the MAE and bias revealed a difference in performance between gliders Henry and Grampus. Grampus missions had a positive bias and MAE of 30-60%, whereas the Henry missions had both positive and negative biases of less than 12%, and a MAE of < 20%. Both the bias and MAE were relatively unchanged when only the flag = 1 and type = 1 data were considered. Despite the larger errors in the Grampus retrievals, all the glider missions taken together have a bias of 0.01% and a MAE of 21%. These results are in rough agreement with Fan et al. [20] who applied their modified inversion method to profiles of radiometric data within the Ligurian Sea.

There are two major differences between gliders Grampus and Henry. First, the glider casing where the optical sensors are located on Grampus is black, but yellow on Henry. Second, the backscattering sensor on Henry is on the bottom of the glider, whereas on Grampus it is approximately 45° from the bottom towards the right hand side of the glider. In the summer of 2016, both Henry and Grampus were deployed at the same time within the Gulf of Maine. The data from that experiment (not shown) revealed when the gliders were within 2 km of each other and 5m vertically within the water column, the radiometric sensors had a MAE of approximately 15% and the backscattering sensors had a MAE of approximately 28%. Hence, the difference observed between the backscattering sensors could explain the higher bias and MAE in the inversion results for Grampus. We plan to do a more controlled comparison between the sensors within the lab to better understand this difference.

3.2.3 Depth-corrected $a$ inversion results

Before discussing the retrieval of the absorption coefficient it is useful to consider Gershun’s law [43], which can be written as

a = μ_{d} K_{d} \frac{1 - R (K_{u} / K_{d})}{1 + R (K_{u} / K_{d})}

where the K’s are the upward (u) and downward (d) irradiance attenuation coefficients, R is the irradiance ratio (upward/downward), and the

μ

’s are the average cosines of the upward and downward radiance distributions. Because R is small (in our case R is of the order of 0.01), this can be approximated by

a \approx μ_{d} K_{d}

. Furthermore, a good approximation to

μ_{d}

is the cosine of the solar zenith angle beneath the water surface:

μ_{0}

. This equation can usually yield the absorption coefficient with an uncertainty within 10-15%. In our case, we have measurements or direct retrievals of all of the quantities in the exact equation for a and therefore a better approximation than the somewhat crude estimate

a = μ_{0} K_{d}

, especially when reasonable agreement between the measured and retrieved backscattering coefficients has been obtained in the inversion. Thus, we expect our estimate of a is at least as good as the crude one and likely much better.

Unfortunately, there is no sensor on either of the gliders measuring the absorption coefficient, therefore we cannot do a direct comparison to investigate the accuracy of the estimation of $a$ from the inversion method. Within the Gulf of Maine, $a$ is expected to co-vary with respect to phytoplankton chlorophyll-a absorption, hence we can compare our estimated $a$ with the chlorophyll-a (Chl) fluorescence to get a qualitative measure of the performance. In order to put a quantitative measure on this, we used an empirical relationship to estimate the Chl concentration from the Chl fluorescence (this relationship was derived between 26 samples of extracted Chl taken at launch and recovery of the gliders and the Chl fluorescence measured by the glider fluorometers, with $R^{2} = 0.65$ ). The absorption due to phytoplankton at 443 nm, $a_{p h} (443)$ , can then be estimated from Chl concentration [44] (see their Eq. (6).

Figure 5 shows the relationship between the inverted $a (443)$ and measured Chl fluorescence and estimated $a_{p h} (443)$ for Henry mission 4, which took place 27^th April - 20^th May 2009 and sampled the spring bloom. It can be seen that Chl fluorescence (and hence $a_{p h} (443)$ ) co-varies with $a (443)$ as expected in the phytoplankton dominated waters of the Gulf of Maine [45]. There are other materials affecting $a (443)$ , as shown by the position of the data cloud in comparison to the reference 1:1 line between $a (443)$ and $a_{p h} (443)$ which indicates the contribution of $a_{p h} (443)$ to the total absorption. A component of the total absorption will be due to CDOM which will co-vary with Chl concentration, however in coastal shelf seas there may be additional CDOM from terrestrial environments which does not co-vary with Chl [46]. This is true within the Gulf of Maine which has been found to have component of $a$ due to CDOM which is not explained by Chl [39] and which seen increased levels of CDOM over the last century that are related to increased freshwater inputs [38]. It is worth noting that our estimate of $a_{p h} (443)$ is from two empirical relationships which have their own errors, hence the quantitative agreement observed between $a (443)$ and $a_{p h} (443)$ in Fig. 5 is encouraging.

Fig. 5 Relationship between $a (443)$ , the measured Chl fluorescence and the corresponding estimated $a_{p h} (443)$ . The 1:1 line between $a (443)$ and $a_{p h} (443)$ is shown as reference.

Download Full Size | PDF

4. Conclusions

In this study, we have presented a modified version of the inversion method of Gordon and Boynton [31] to estimate profiles of absorption and backscattering coefficients from spectral downwelling irradiance and upwelling radiance measurements as measured from a single platform such as a glider or profiling float. Our version of the inversion method performed similarly with the original on modeled data, estimating the absorption and backscattering coefficients with a mean MARE (on linear data) of 2.26% and 4.61%, respectively.

Part of this method involved an automated QC procedure for $E_{d}$ and $L_{u}$ glider-measured profiles which involves no pre-screening of the data for varying illumination due to sky conditions. The quality of each measurement and the full yo containing a series of measurements at different depths were classified to give an indication of sub-optimal conditions [40]. We applied the combined QC procedure and modified inversion method to over 1,300 glider measurements within the optically complex water of the Gulf of Maine. A comparison of the estimated and measured $b_{b}$ revealed variability between missions, but overall, we found a bias of 0.01% and a MAE of 21% on log-transformed data. The original study of Gordon and Boynton [31] did not evaluate the inversion method on field data, but the method was later validated in clear ocean [32]. Our results are comparable with a study which applied a different modified version of Gordon and Boynton [31] within the Ligurian Sea [20].

While we are not able to directly validate the inverted absorption coefficients for the glider data, our modelling work and Gordon and Boynton [31] found it was possible to estimate $a$ more accurately than $b_{b}$ . This inversion method, in particular, relies on numerically solving the RTE, hence if we are recovering $b_{b}$ to a reasonable degree of accuracy and are confident in our method for solving the RTE, we can assume the error in our recovered $a$ is, at least, no worse than the $b_{b}$ and likely smaller. The comparison with fluorometry derived phytoplankton absorption is encouraging, especially given the strong covariance between $a$ and chlorophyll-a fluorescence in the Gulf of Maine.

These results are promising. They show the feasibility of estimating IOPs from $E_{d}$ and $L_{u}$ under all sky conditions with an automated approach including QC of the data. This procedure can be used as a building block in the study of marine ecosystems, such as estimating different biogeochemical parameters (e.g. particulate organic carbon) or studying phytoplankton dynamics (e.g. growth). Combining these methods with the temporal and spatial coverage achieved by gliders will allow detailed studies of marine ecosystems throughout the euphotic zone.

Funding

National Aeronautics and Space Administration (NASA) (Grant Nos. NNX14AM77G, NNX14AQ43A, NNX14AQ41G, NNX17AI77G).

Acknowledgments

We would like to acknowledge Teledyne Webb Research’s technical support.

References

1. IOCCG, “Remote Sensing of Inherent Optical Properties: Fundamentals, Tests of Algorithms and Applications,” in Reports of the International Ocean-Colour Coordinating Group, No. 5, Z. P. Lee, ed. (2006), Vol. 5, p. 126.

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

3. P. J. Werdell, L. I. W. Mckinna, E. Boss, S. G. Ackleson, S. E. Craig, W. W. Gregg, Z. Lee, S. Maritorena, C. S. Roesler, C. S. Rousseaux, D. Stramski, J. M. Sullivan, M. S. Twardowski, M. Tzortziou, and X. Zhang, “An overview of approaches and challenges for retrieving marine inherent optical properties from ocean color remote sensing,” Prog. Oceanogr. 160, 186–212 (2018). [CrossRef]

4. M. Ramírez-Pérez, M. Twardowski, C. Trees, J. Piera, and D. McKee, “Inversion of In Situ Light Absorption and Attenuation Measurements to Estimate Constituent Concentrations in Optically Complex Shelf Seas,” J. Geophys. Res. Oceans 123(1), 720–737 (2018). [CrossRef]

5. O. Schofield, T. Bergmann, M. J. Oliver, A. Irwin, G. Kirkpatrick, W. P. Bissett, M. A. Moline, and C. Orrico, “Inversion of spectral absorption in the optically complex coastal waters of the Mid-Atlantic Bight,” J. Geophys. Res. 109(C12), C12S04 (2004). [CrossRef]

6. G. Zheng and D. Stramski, “A model for partitioning the light absorption coefficient of suspended marine particles into phytoplankton and nonalgal components,” J. Geophys. Res. Oceans 118(6), 2977–2991 (2013). [CrossRef]

7. Z. Lee, V. P. Lance, S. Shang, R. Vaillancourt, S. Freeman, B. Lubac, B. R. Hargreaves, C. Del Castillo, R. Miller, M. Twardowski, and G. Wei, “An assessment of optical properties and primary production derived from remote sensing in the Southern Ocean (SO GasEx),” J. Geophys. Res. Oceans 116, 1–15 (2011).

8. M. K. Barnes, G. H. Tilstone, T. J. Smyth, D. J. Suggett, R. Astoreca, C. Lancelot, and J. C. Kromkamp, “Absorption-based algorithm of primary production for total and size-fractionated phytoplankton in coastal waters,” Mar. Ecol. Prog. Ser. 504, 73–89 (2014). [CrossRef]

9. D. Stramski, R. A. Reynolds, M. Kahru, and B. G. Mitchell, “Estimation of particulate organic carbon in the ocean from satellite remote sensing,” Science 285(5425), 239–242 (1999). [CrossRef] [PubMed]

10. A. Cunningham, L. Ramage, and D. McKee, “Relationships between the inherent optical properties and the depth of penetration of solar radiation in optically complex waters,” J. Geophys. Res. 118, 2310–2317 (2013).

11. C. Mitchell and A. Cunningham, “Remote sensing of spatio-temporal relationships between the partitioned absorption coefficients of phytoplankton cells and mineral particles and euphotic zone depths in a partially mixed shelf sea,” Remote Sens. Environ. 160, 193–205 (2015). [CrossRef]

12. C. Mitchell, A. Cunningham, and D. McKee, “Derivation of the specific optical properties of suspended mineral particles and their contribution to the attenuation of solar irradiance in offshore waters by ocean color remote sensing,” J. Geophys. Res. Oceans 121(1), 104–117 (2016). [CrossRef]

13. H. R. Gordon, “Spectral Variations of the Volume Scattering Function at Large Angles in Natural Waters,” J. Opt. Soc. Am. 64(6), 773–775 (1974). [CrossRef]

14. H. R. Gordon, “Absorption and Scattering Estimates from Irradiance Measurements: Monte Carlo Simulations,” Limnol. Oceanogr. 36(4), 769–777 (1991). [CrossRef]

15. R. A. Leathers and N. J. McCormick, “Ocean inherent optical property estimation from irradiances,” Appl. Opt. 36(33), 8685–8698 (1997). [CrossRef] [PubMed]

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

17. 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(3), 628–641 (2000). [CrossRef]

18. E. Rehm and C. D. Mobley, “Estimation of hyperspectral inherent optical properties from in-water radiometry: error analysis and application to in situ data,” Appl. Opt. 52(4), 795–817 (2013). [CrossRef] [PubMed]

19. D. McKee, A. Cunningham, and S. Craig, “Estimation of absorption and backscattering coefficients from in situ radiometric measurements: theory and validation in case II waters,” Appl. Opt. 42(15), 2804–2810 (2003). [CrossRef] [PubMed]

20. Y. Fan, W. Li, V. S. Calzado, C. Trees, S. Stamnes, G. Fournier, D. McKee, and K. Stamnes, “Inferring inherent optical properties and water constituent profiles from apparent optical properties,” Opt. Express 23(15), A987–A1009 (2015). [CrossRef] [PubMed]

21. K. S. Johnson, W. M. Berelson, E. S. Boss, Z. Chase, H. Claustre, S. R. Emerson, N. Gruber, A. Kortzinger, M. J. Perry, and S. C. Riser, “Observing Biogeochemical Cycles At Global Scales With Profiling Floats and Gliders: Prospects for a global array,” Oceanography (Wash. D.C.) 22(3), 216–225 (2009). [CrossRef]

22. K. Niewiadomska, H. Claustre, L. Prieur, and F. d’Ortenzio, “Submesoscale physical-biogeochemical coupling across the Ligurian current (northwestern Mediterranean) using a bio-optical glider,” Limnol. Oceanogr. 53(5part2), 2210–2225 (2008). [CrossRef]

23. O. Schofield, R. Chant, B. Cahill, R. Castelao, D. Gong, A. Kahl, J. Kohut, M. Montes-Hugo, R. Ramadurai, P. Ramey, X. Yi, and S. Glenn, “The Decadal View of the Mid-Atlantic Bight from the COOLroom: Is Our Coastal System Changing?” Oceanography (Wash. D.C.) 21(4), 108–117 (2008). [CrossRef]

24. J. Bouffard, A. Pascual, S. Ruiz, Y. Faugère, and J. Tintoré, “Coastal and mesoscale dynamics characterization using altimetry and gliders: A case study in the Balearic Sea,” J. Geophys. Res. 115(C10), C10029 (2010). [CrossRef]

25. L. Merckelbach, D. Smeed, and G. Griffiths, “Vertical water velocities from underwater gliders,” J. Atmos. Ocean. Technol. 27(3), 547–563 (2010). [CrossRef]

26. F. Bourrin, G. Many, X. Durrieu de Madron, J. Martín, P. Puig, L. Houpert, P. Testor, S. Kunesch, K. Mahiouz, and L. Béguery, “Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions,” Cont. Shelf Res. 109, 135–149 (2015). [CrossRef]

27. V. S. Hemsley, T. J. Smyth, A. P. Martin, E. Frajka-Williams, A. F. Thompson, G. Damerell, and S. C. Painter, “Estimating Oceanic Primary Production Using Vertical Irradiance and Chlorophyll Profiles from Ocean Gliders in the North Atlantic,” Environ. Sci. Technol. 49(19), 11612–11621 (2015). [CrossRef] [PubMed]

28. N. Briggs, M. J. Perry, I. Cetinić, C. Lee, E. D’Asaro, A. M. Gray, and E. Rehm, “High-resolution observations of aggregate flux during a sub-polar North Atlantic spring bloom,” Deep. Res. Part I Oceanogr. Res. Pap. 58, 1031–1039 (2011).

29. C. A. Brown, Y. Huot, M. J. Purcell, J. J. Cullen, and M. R. Lewis, “Mapping coastal optical and biogeochemical variability using an autonomous underwater vehicle and a new bio-optical inversion algorithm,” Limnol. Oceanogr. Methods 2(8), 262–281 (2004). [CrossRef]

30. M. A. Moline, I. Robbins, B. Zelenke, W. S. Pegau, and H. Wijesekera, “Evaluation of bio-optical inversion of spectral irradiance measured from an autonomous underwater vehicle,” J. Geophys. Res. Oceans 117, 1–12 (2012).

31. 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(18), 3886–3896 (1998). [CrossRef] [PubMed]

32. H. R. Gordon, M. R. Lewis, S. D. McLean, M. S. Twardowski, S. A. Freeman, K. J. Voss, and G. C. Boynton, “Spectra of particulate backscattering in natural waters,” Opt. Express 17(18), 16192–16208 (2009). [CrossRef] [PubMed]

33. M. Jonasz and G. R. Fournier, Light Scattering by Particles in Water (Academic Press, 2007).

34. C. D. Mobley, “Fast light calculations for ocean ecosystem and inverse models,” Opt. Express 19(20), 18927–18944 (2011). [CrossRef] [PubMed]

35. 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(2), 417–427 (1975). [CrossRef] [PubMed]

36. C. D. Mobley, B. Gentili, H. R. Gordon, Z. Jin, G. W. Kattawar, A. Morel, P. Reinersman, K. Stamnes, and R. H. Stavn, “Comparison of numerical models for computing underwater light fields,” Appl. Opt. 32(36), 7484–7504 (1993). [CrossRef] [PubMed]

37. J. Ronald and V. Zaneveld, “Remotely sensed reflectance and its dependence on vertical structure: a theoretical derivation,” Appl. Opt. 21(22), 4146–4150 (1982). [CrossRef] [PubMed]

38. W. M. Balch, T. Huntington, G. Aiken, D. T. Drapeau, B. C. Bowler, L. C. Lubelczyk, and K. Butler, “Toward a quantitative and empirical dissolved organic carbon budget for the Gulf of Maine, a semienclosed shelf sea,” Global Biogeochem. Cycles 30(2), 268–292 (2016). [CrossRef]

39. W. M. Balch, D. T. Drapeau, B. C. Bowler, E. S. Booth, J. I. Goes, A. Ashe, and J. M. Frye, “A multi-year record of hydrographic and bio-optical properties in the Gulf of Maine: I. Spatial and temporal variability,” Prog. Oceanogr. 63(1-2), 57–98 (2004). [CrossRef]

40. E. Organelli, H. Claustre, A. Bricaud, C. Schmechtig, A. Poteau, X. Xing, L. Prieur, F. D’Ortenzio, G. Dall’Olmo, and V. Vellucci, “A novel near-real-time quality-control procedure for radiometric profiles measured by bio-argo floats: Protocols and performances,” J. Atmos. Ocean. Technol. 33(5), 937–951 (2016). [CrossRef]

41. B. N. Seegers, R. P. Stumpf, B. A. Schaeffer, K. A. Loftin, and P. J. Werdell, “Performance metrics for the assessment of satellite data products: An ocean color case study,” Opt. Express 26(6), 7404–7422 (2018). [CrossRef] [PubMed]

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

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

44. A. Bricaud, A. Morel, M. Babin, K. Allali, and H. Claustre, “Variations of light absorption by suspended particles with chlorophyll a concentration in oceanic (case 1) waters: Analysis and implications for bio-optical models,” J. Geophys. Res. Oceans 103(C13), 31033–31044 (1998). [CrossRef]

45. W. M. Balch, D. T. Drapeau, B. C. Bowler, and T. Huntington, “Step-changes in the physical, chemical and biological characteristics of the Gulf of Maine, as documented by the GNATS time series,” Mar. Ecol. Prog. Ser. 450, 11–35 (2012). [CrossRef]

46. H. R. Gordon and A. Morel, Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery. A Review (Springer-Verlag, 1983).

IOP Profiles	$a_{0}$ (m $^{- 1}$ )	$a_{1}$ (m $^{- 1}$ )	$b_{0}$ (m $^{- 1}$ )	$b_{1}$ (m $^{- 1}$ )	$σ_{a}$ (m)	$σ_{b}$ (m)	$z_{a}$ (m)	$z_{b}$ (m)
aub2	0.1	0	0.0072	0.0072	1	1	-	2
a2b2	0.1	0.1	0.0072	0.0072	1	1	3	1

IOP Profile	Solar Zenith Angle	Correct vs incorrect phase function	MARE $a$ (%) from this study		MARE $b_{b}$ (%) from this study	MARE $b_{b}$ (%) from [31]
aub2	0	Correct	0.62	0.68	3.48	2.81
		Incorrect	4.07	1.65	6.55	5.99
	60	Correct	0.37	1.20	3.55	3.16
		Incorrect	3.91	2.30	9.25	8.42
a2b2	0	Correct	0.80	1.72	1.65	2.32
		Incorrect	4.02	1.45	4.13	4.88
	60	Correct	0.65	1.92	1.10	1.79
		Incorrect	3.65	1.59	7.13	5.87

IOP Profiles	$a_{0}$ (m $^{- 1}$ )	$a_{1}$ (m $^{- 1}$ )	$b_{0}$ (m $^{- 1}$ )	$b_{1}$ (m $^{- 1}$ )	$σ_{a}$ (m)	$σ_{b}$ (m)	$z_{a}$ (m)	$z_{b}$ (m)
aub2	0.1	0	0.0072	0.0072	1	1	-	2
a2b2	0.1	0.1	0.0072	0.0072	1	1	3	1

IOP Profile	Solar Zenith Angle	Correct vs incorrect phase function	MARE $a$ (%) from this study		MARE $b_{b}$ (%) from this study	MARE $b_{b}$ (%) from [31]
aub2	0	Correct	0.62	0.68	3.48	2.81
		Incorrect	4.07	1.65	6.55	5.99
	60	Correct	0.37	1.20	3.55	3.16
		Incorrect	3.91	2.30	9.25	8.42
a2b2	0	Correct	0.80	1.72	1.65	2.32
		Incorrect	4.02	1.45	4.13	4.88
	60	Correct	0.65	1.92	1.10	1.79
		Incorrect	3.65	1.59	7.13	5.87

Optical inversions of the water column based on glider measurements

Abstract

1. Introduction