Sensitivity analysis of volume scattering phase functions

: To solve the radiative transfer equation and relate inherent optical properties (IOPs) to apparent optical properties (AOPs), knowledge of the volume scattering phase function is required. Due to the difficulty of measuring the phase function, it is frequently approximated. We explore the sensitivity of derived AOPs to the phase function parameterization, and compare measured and modeled values of both the AOPs and estimated phase functions using data from Monterey Bay, California during an extreme “red tide” bloom event. Using in situ measurements of absorption and attenuation coefficients, as well as two sets of measurements of the volume scattering function (VSF), we compared output from the Hydrolight radiative transfer model to direct measurements. We found that several common assumptions used in parameterizing the radiative transfer model consistently introduced overestimates of modeled versus measured remote-sensing reflectance values. Phase functions from VSF data derived from measurements at multiple wavelengths and a single scattering single angle significantly overestimated reflectances when using the manufacturer-supplied corrections, but were substantially improved using newly published corrections; phase functions calculated from VSF measurements using three angles and three wavelengths and processed using manufacture-supplied corrections were comparable, demonstrating that reasonable predictions can be made using two commercially available instruments. While other studies have reached similar conclusions, our work extends the analysis to coastal waters dominated by an extreme algal bloom with surface chlorophyll concentrations in excess of 100 mg m − 3 .


Introduction
In the field of ocean optics, volume scattering phase functions play a key role in the solution of the radiative transfer equation.Derived from basic energy conservation, the radiative transfer equation (RTE) is of principle importance in ocean optics because it describes how light propagates through optically thick media, given a set of inherent optical properties (IOPs).With the radiances and irradiances found from the RTE, one can determine the apparent optical properties (AOPs) of a medium [1].The apparent optical properties depend on the in situ light conditions, and include ratios of radiative quantities and reflectances, properties which minimize the effects of ambient light and are in general readily measurable.IOPs, on the other hand, are properties of a medium that are independent of ambient light conditions, and include absorption (a), scattering (b), and attenuation (c) coefficients.IOPs are generally much harder to measure than AOPs, often requiring specialized instruments [2].A fundamental goal of ocean optics is to achieve optical closure, wherein IOPs and AOPs converge to a single solution of the radiative transfer equation, allowing estimation of these parameters from remote sensing data [3].Optical closure gives us verification of the integrity of measurements of optical quantities and is paramount for evaluating the accuracy of numerical solutions to the RTE.However, to solve the radiative transfer equation, one must have knowledge of the volume scattering phase function.A common form of the radiative transfer equation used for ocean optics applications is shown below [1]: where L is the radiance, ζ is the optical depth, ω o is the single scattering albedo, μ is the cosine of the nadir angle, φ is the azimuthal angle, c is the attenuation coefficient, and S is the source function.Note that this solution assumes a plane-parallel geometry, and ignores the effects of polarization.
The phase function term β  is inconveniently located inside an integral in the radiative transfer equation (Eq.( 1)), and is required to obtain a numerical solution.A physical understanding of the phase function is best derived from another quantity, the volume scattering function, divided by the scattering coefficient ( / b ).The volume scattering function (VSF) is defined as the scattering cross section per unit volume, and can be thought of as the angular distribution of light scattering off a given volume [1].Integrating the volume scattering function over 4π steradians gives the scattering coefficient, b, while integrating over the backwards hemisphere (relative to the incident light field) gives the backscattering coefficient, b b .The phase function, with units of sr −1 , can then best be described as a normalized version of the volume scattering function, describing the angular distribution of light scattering off of a medium, independent of the intensity of incident light sources.
Due to the difficulty of measuring the phase function and the lack of commercial instruments available that can directly measure it, the phase function is typically estimated indirectly, or is approximated based on typical conditions [4,5].Examples of such approximations include the Petzold phase function, based on measurements from San Diego Harbor and often considered to be "typical" of coastal marine waters [6], and the Fournier-Forand (FF) phase function, a continuous phase function derived from Mie scattering that takes the index of refraction and hyperbolic slope of a Junge fit of the particle size distribution as arguments [7].A Fournier-Forand particulate phase function may also be constructed with knowledge of the backscatter fraction, the ratio of particulate backscatter and scattering coefficients, for a given water column [8].These approximations of phase functions are frequently used and are included as options in the Hydrolight radiative transfer model, but while convenient, these phase functions may not be appropriate for a given water mass.Here we aim to determine how modifying the phase function term in the RTE, based on IOPs measured with commercially available instruments, affects the numerical calculations of remote sensing reflectances.Using phase function measurements derived from multiple instruments, we compare the calculated reflectances to see which combination of measurements best approximates the measured reflectances, or achieves optical closure.Since phase functions are often significantly dependent on the composition of ocean waters, we also seek to determine the errors introduced when using estimates rather than direct or inferred measurements of the phase function.Finally we compare phase functions derived from measurements made with two types of commercially available instruments: HOBI Labs backscatter meters and a WET Labs VSF meter, to identify which instrument or combination of instruments most accurately retrieve the phase function for our study site.
We focused primarily on assessing the sensitivity of radiative transfer models to phase function measurements in a complex coastal environment.Coastal environments, such as our site in Monterey Bay, California, experience large variations in chlorophyll concentrations and frequent harmful algal blooms [9].Environments such as this are primarily driven by upwelling and relaxation cycles and may exhibit large seasonable variability in phytoplankton abundances [10].The optical properties of these waters differ significantly from the standard values assumed in the formulation of phase functions.Previous authors have directly compared backscatter meters in both open-ocean and coastal sites [11][12][13]; here we extend these comparisons to include red tide conditions with optical properties dominated by the dinoflagellate Akashiwo sanguinea in extreme Case 1 waters.

Data set
For our sensitivity analysis of phase functions, we used in situ measurements from the 2006 COAST project in Monterey Bay, California [14] aboard the R/V John Martin.This data set was collected in a region exhibiting sharp spatial gradients between a dinoflagellatedominated community (Akashiwo sanguinea) and more typical coastal waters.Optical properties for this site are extreme compared to standard case I or case II waters, with surface chlorophyll exceeding 100 mg m −3 .Despite the high biomass at this time of year Monterey Bay exhibits Case I properties with minimal input of sediments and other optical constituents from upwelling or fluvial inputs, leading to dominance by particulate organic material [15].
Data were collected on September 7, 2006, and include multiple measurements of the volume scattering function (Table 1).Measurements included vertical profiles of pressure, temperature and salinity from a SeaBird SBE49 CTD, equipped with a WETLabs WETStar fluorometer.A WETLabs 0.25 m pathlength ac-s was deployed concurrently to obtain total absorption (a) and attenuation (c) at 83 distinct wavelengths between 401.5 and 752.5nm.A WETLabs 0.25 m pathlength ac-9 was deployed with a Gelman 0.2 µm filter inline to estimate colored dissolved organic material (CDOM) absorption at 9 wavelengths.The ac-s and ac-9 data were corrected for temperature, salinity and pure water [16] with pure-water calibrations at the beginning and end of the day, and the ac-s was also corrected for particulate absorption using the 713.9 nm channel as reference.A Satlantic Hyperpro II freefalling optical profiler was deployed to obtain vertical profiles of in water radiances and reflectances, with one cast conducted using a floatation collar to obtain near-surface data.The data were processed using the Satlantic ProSoft software to estimate remote sensing reflectances (R rs ) following standard protocols.Briefly, three profiles were conducted from approximately 0-10 m with "shutter darks" enabled.Post-processing of level 2 (L2) data to level 3 (L3) and then level 4 (L4) extrapolated water-leaving radiances across the air-sea interface using calculated attenuation coefficients over the first 3-5 m of the cast, and R rs was calculated using water-leaving radiance (L u ) and the deck-reference (E s ) sensor.Data were filtered for excessive tilt, and casts were considered acceptable if the extrapolated abovewater E d values were within 10% of the E s values for the wavelength range 400-700 nm.Data from the three casts and the surface floating acquisition were combined to generate a single R rs data set.
Data for estimating the volume scattering function were provided by two separate devices.The first, a WETLabs ECO-VSF3, measured the volume scattering function in the backwards direction at 3 angles (100°, 125°, and 150°) and 3 wavelengths (450, 530, and 650 nm).From the ECO-VSF-3 measurements of the volume scattering function (β), numerical integration was used to obtain the backscattering coefficient (b b ).A HOBI-Labs HydroScat-6 was also deployed.The HydroScat-6 measures the VSF at 6 wavelengths (442, 470, 550, 589, 620, and 671nm) and at a single angle (140°).The HydroScat was deployed with 2 single-wavelength HOBI-Labs a-Betas (420 and 510 nm), adding an additional 2 wavelengths of VSF measurements.Data were collected from the same station within a few minutes of each other as three separate profiles.The CTD, ac-s, ac-9, and ECO-VSF were deployed together, as were the HOBI instruments, both using a hydrowire.The HyperPro II was profiled independently multiple times, spatially and temporally concurrent with the IOP measurements to the extent possible.Total chlorophyll and CDOM absorption were measured from discrete samples collected by bucket (surface) and CTD (subsurface).Chlorophyll was extracted in 90% acetone and analyzed spectrophotometrically using a 10-cm cell on a Cary-Varian 50 spectrophotometer.CDOM was measured on the same spectrophotometer after filtration with a 0.2 µm filter.The extracted chlorophyll data were used to calibrate the WetStar fluorometer, without a separate correction for non-photochemical quenching.This correction was assumed to have a negligible impact on modeled R rs since chlorophyll was only used to estimate inelastic scattering contributions of fluorescence.
The vendor-supplied software for the HydroScat instrument obtains b b using the Oishi (1990) proportionality argument for β and b b at angles around 120°, adjusted to the 140° angle of the HS6 instrument [17,18].The software also provides a sigma-correction, but we elected to apply independent estimates of the sigma-correction term, since recent literature has demonstrated that the manufacture-supplied sigma-correction is inadequate [19].We first sigma-corrected our data following the vendor-specified method as modified by [19] by multiplying the raw data values by σ, where σ was given by the approximation used in the HydroScat-6 manual that 1 exp( ) where k 1 and k exp are the instrument calibration constants and K bb is the attenuation function of the medium, an AOP which can be approximated as 0.
. The absorption and scattering coefficients (a and b) were known for each of our b b values.The magnitude of the sigma corrections using this approach varied significantly with the wavelength and depth for which the b b measurements were derived.At deeper, lower chlorophyll concentration waters the correction was barely greater than 9%, while in the highest chlorophyll concentrations near the surface, the correction exceeded 260% for the shorter wavelengths.We subsequently corrected the data using the approached described by Doxaran et al. 2016 [13] with a modified K bb function: where nw a is wavelength-specific non-water absorption (from the ac-s).This led to corrections of between ~2 and 158% compared to the uncorrected data.

Constructing phase functions
Due to the innate variability of measured IOPs of given water column as a function of depth and wavelength, we opted to construct different phase functions for each discrete depth and wavelength.To obtain continuous phase functions from our discrete measurements, we used the computed backscattering fractions, / , to construct Fournier-Forand particulate phase functions for each depth and wavelength.The Fournier-Forand function is dependent on both the index of refraction, n, and the hyperbolic slope of the particle size distribution, μ, and since B is dependent on both these quantities, we can approximate a Fournier-Forand phase function with a single B value [20].
The backscattering fraction, B, was obtained separately from both the ECO-VSF3 and HydroScat-6 data.We also calculated Fournier-Forand phase functions using B = 0.015, which is a value of the backscattering fraction that falls near the middle of our range of measured B values, and is representative of a reasonable estimate for the phase function if direct measurements were not available.We also included the Petzold Average Particle phase function, which is frequently used in the solution of the radiative transfer equation when the phase function is unknown.

Solving the radiative transfer equation
We used the Hydrolight © 4.2 and 5.3 software from Sequoia Scientific for surface reflectance calculations, using measured IOP values and phase functions as inputs.Hydrolight uses a numerical model to solve the radiative transfer equation and compute the spectral radiance distribution and derived quantities of the medium [21,22].For each calculation, all user supplied quantities were kept constant except for the phase function.For each run, we used absorption and attenuation data from the ac-s measurements as inputs for the IOPs of the ocean water.We used Pope and Fry [23] values with saltwater scattering for the inherent optical properties of pure water and the built-in pure water phase function to specify pure water IOPs for each Hydrolight calculation.For inputs of chlorophyll biomass, we entered the measured values of chlorophyll concentrations from the WETStar fluorometer data, calibrated with discrete chlorophyll measurements.Although chlorophyll was not used directly in Hydrolight (since we supplied the measured IOPs), it was used for one of the inelastic scattering terms (chlorophyll fluorescence, CDOM fluorescence, and Raman scatter).We opted to use the "Chl model for case 1 water" that was built-in to Hydrolight.To include the effects of CDOM fluorescence, we entered our measured values.
We also entered information for a variety of corrections included in Hydrolight including atmospheric and air-water surface conditions.We entered the recorded weather conditions in Monterey Bay on September 7, 2006 [24], including atmospheric pressure, wind speed, horizontal visibility, relative humidity, and cloud cover.The incident solar irradiance was estimated from the time, latitude and longitude of the location for the 2006 COAST data set.We used the semi-empirical RADTRAN-X model to compute sky irradiances, and used the default HCNRAD model as the angular pattern of the sky radiance.We assumed an infinitely deep water column, since we found that varying our depth and bottom reflectances didn't significantly affect our calculated AOP values.Comparison of RADTRAN-X output with the HyperPro II E s sensor for 400-700 nm demonstrated 3% deviation between the measurements with no obvious bias (r 2 = 0.92; slope = 1.00 for RADTRAN-X E d versus Hydrolight E s values).
We used the phase functions described in the previous section, the ECO-VSF3, HydroScat-6, Fournier-Forand (B = 0.015), and Petzold average particle phase functions, to estimate R rs between 400 and 700 nm. .Calculations were performed at depth intervals of 0.5 m until a maximum depth of 16 m, for wavelengths between 400 nm and 700 nm with a bandwidth of 30 nm.For the HydroScat-6 data using the improved sigma-correction, wavelength resolution was increased to 10 nm bandwidth.The modeled remote-sensing reflectance was used for comparisons between model runs and observations from the HyperPro II data.

Results
The measured IOPs of the water column varied significantly as a function of depth and wavelength (Fig. 1).Whereas the absorption coefficients appeared to vary strongly with both wavelength and depth, the total attenuation seemed to primarily depend only on depth.This variability in absorption and attenuation was also observed with our measured volume scattering functions and backscattering coefficients [25] and was consistent with the measured chlorophyll concentrations.At the very near surface, chlorophyll was up to 370 mg m −3 decreasing to 81 mg m −3 when depth-averaged over the first meter, and declining to 13 and 3.5 mg m −3 at 6 and 12 m depth.We constructed Fournier-Forand phase functions using particulate backscattering fractions, and the variability of these IOPs caused our phase functions to vary as well.Example calculated phase functions are provided in Fig. 2. The measured remote-sensing reflectance spectra from the HyperPro II data were typical of "red tide" conditions, with maxima between 550 and 600 nm and 680-750 nm (Fig. 3).There was considerable spatial and temporal heterogeneity, but the overall shape of the reflectance spectra remained uniform, based on the standard deviation of the multiple HyperPro II casts.
After entering our IOP measurements as inputs in Hydrolight simulations, the modeled R rs spectra were similar in shape (as expected, since only the phase-function was varied) and generally captured the overall shape of the reflectance spectrum, with a pronounced peak between 550 and 600 nm (Fig. 3).The pronounced peak in the near infrared was located outside of the range of wavelengths used in the Hydrolight models, so it is not represented in the calculated values.
To assess the goodness of fit of the models computed with different phase functions compared to the measured R rs values, we used a reduced chi-squared test.While none of our models perfectly matched the measured values, it is clear that they are capturing the general R rs shape.The reflectances computed with HydroScat-6 data using the standard sigmacorrection are considerably larger than the measured R rs measurements, and have the largest chi-squared value (Fig. 3).This is significantly improved when using the sigma-correction based on [13], with a final χ 2 red = 0.255.Reduced chi-squared values using the standard corrections from the ECO-VSF3 were of χ 2 red = 0.778, comparable to the HydroScat-6 results.
The R rs values computed using the two estimates of the phase function (Fournier-Forand (B = 0.015) and Petzold) were significantly larger than the measured values.Representative b bp values for 1 and 6 dbar are provided in Fig. 4. Previous analyses have demonstrated that the ECO-VSF3 and HydroScat instruments provide comparable estimates of b bp when deployed concurrently and when data are processed using comparable methods [11,12].Consistent with those findings, our data exhibit <15% variability for the full vertical profile when comparing b bp at 530 nm (directly measured by the ECO-VSF3 and estimated at 530 nm using a power-law function from the corresponding HydroScat wavelengths) after applying the sigma-correction from [13].In contrast, using the manufacturer-suggested sigma-correction, there was reasonable agreement at low b bp values (low chlorophyll values) with deviations between the two instruments increasing rapidly with chlorophyll concentrations above approximately 50 mg m −3 .This is consistent with the findings of [13] who demonstrated that the ECO instruments, which use short optical pathlengths, require negligible pathlength scattering corrections and are less sensitive to absorption corrections.In contrast, the HydroScat-series of instruments, and in particular the physically large HydroScat-6 with its correspondingly longer optical pathlength, must be corrected for both absorption and scattering.As demonstrated by [13] the standard correction overestimates the scattering contribution, resulting in increasingly large deviations from the ECO-VSF3 with increasing chlorophyll, but instrument closure can be achieved when the proper sigmacorrection is applied (Fig. 5).Here we show that the correction scheme of [13] works well in both highly scattering waters and also in waters strongly dominated by phytoplankton, which, relative to inorganic particles, exhibit strong wavelength-specific absorption and weaker scattering properties [15].

Analysis and conclusions
It is well known that varying the phase function can significantly influence the calculated values of AOPs from Hydrolight runs.The differences in magnitude of calculated R rs demonstrate that model outputs for these waters are sensitive to the phase function term, as expected.Estimates using two commercially available instruments provided estimated R rs that varied by a factor of 2 in this complex coastal environment when using the manufacturersupplied correction schemes.The Hydrolight models that best predicted the Rrs spectrum were obtained using phase function measurements from the HydroScat instruments, when corrected using the methods developed by [13], with lower but reasonable agreement when using the ECO-VSF3 data.The other phase function measurements didn't lead to very accurate calculations of the remote sensing reflectances. The where b bw and β w are the backscattering coefficient and volume scattering function of pure water, and χ bb is the conversion factor for particulate backscattering [26].More recent studies suggest that there is the least variability in χ bb at angles around 120° [27].It is possible that the proportionality argument between the VSF and backscattering coefficient may not work as well for the angle used by the HydroScat-6.While this proportionality relation appears to work for typical ocean waters, it has not been applied to extreme algal bloom sites, where it may be expected to perform poorly given the optical design [13].In fact, in waters optically dominated by phytoplankton blooms, the conversion factor may be highly dependent on wavelength [28].Since this data set was taken at a dinoflagellate bloom site, and the HydroScat-6 assumes a constant value for the conversion factor, we expect this wavelengthdependence may introduce errors in the Fournier-Forand phase functions from HydroScat-6 b b values.This is consistent with what others have observed for high chlorophyll waters [19].
In contrast, the ECO-VSF3 is less sensitive to highly scattering and absorbing environments due to its smaller physical size and shorter optical pathlength, but there are fewer wavelengths to constrain the Fornier-Forand phase function.
In the absence of direct phase function measurements one can still somewhat accurately calculate remote sensing reflectances.Using the Fournier-Forand phase function with an estimate of the backscattering fraction yields a reasonable fit to the measured remote sensing reflectances.The Petzold average particle phase function overestimates the reflectance values, but as an approximation for the phase function it may still be ace ptable.These findings are consistent with recent optical closure studies using other radiative transfer models, where the Fournier-Forand phase function was shown to yield more accurate AOPs than the Petzold phase functions [29].
Based on Fig. 3 we failed to fully meet optical closure.The best modeled R rs spectra, based on the HydroScat-6, was close (within 1 standard deviation) of the R rs spectra from direct observations but still underestimated the maximum reflectance peak at ~555 nm.In contrast, the ECO-VSF3 exhibited lower overall goodness of fit (based on reduced chi-square values) but better captured the reflectance maximum.However, given the potential errors in both the observations and the assumptions made for the Hydrolight model, the results are quite good.Optical measurements within an extreme red tide event provide several challenges.There was substantial horizontal and vertical variability, especially in the uppermost layer of the water column where the Akashiwo were actively migrating to the surface [11].None of the optical instruments used in this study are particularly well-designed for these conditions.For example the ac-s and ac-9 pathlength of 0.25 m are optimized for clear-water conditions, while the HyperPro II cannot achieve the approximately 1 cm and 1% For the ECO-VSF3 and HydroScat6 phase functions, we entered the measured b b values and had Hydrolight construct a Fournier-Forand phase function from the backscatter fraction, changing phase functions for each step of / 0.0005 b b B b = =

Fig. 1 .
Fig. 1.IOP data (a p , b p , c p ) as a function of depth and wavelength.The offset at ~552 nm is caused by the change in holographic gratings on the ac-s, and was removed during postprocessing prior to use in Hydrolight.

Fig. 3 .
Fig. 3. Comparison of calculated remote sensing reflectances to measured remote sensing reflectances from Hyperpro data.Gray region represents one standard deviation from averaged R rs values using three vertical profiles and one "floating" data collection.

Fig. 4 .
Fig. 4. Comparison of b bp values as a function of wavelength for 1 and 6 dbar from the HydroScat-6 and ECO-VSF3.

Fig. 5 .
Fig. 5. Comparison of b bp values from the ECO-VSF3 and HydroScat-6, at 530 nm.For the HydroScat-6 data b bp was extrapolated from the 470 and 550 nm channels assuming a powerlaw spectrum.The dashed line indicates a 1:1 relationship.Variability between the two instruments was less than 15% and generally increased with increasing chlorophyll.
HydroScat-6 measures the volume scattering function at a single angle and multiple wavelengths, and uses a proportionality argument for the volume scattering function, β, and the backscattering coefficient, b b ,