Ocean Color Analytical Model Explicitly Dependent on the Volume Scattering Function

An analytical radiative transfer (RT) model for remote sensing reflectance that includes the bidirectional reflectance distribution function (BRDF) is described. The model, called ZTT (Zaneveld-Twardowski-Tonizzo), is based on the restatement of the RT equation by Zaneveld (1995) in terms of light field shape factors. Besides remote sensing geometry considerations (solar zenith angle, viewing angle, and relative azimuth), the inputs are Inherent Optical Properties (IOPs) absorption a and backscattering bb coefficients, the shape of the particulate volume scattering function (VSF) in the backward direction, and the particulate backscattering ratio. Model performance (absolute error) is equivalent to full RT simulations for available high quality validation data sets, indicating almost all residual errors are inherent to the data sets themselves, i.e., from the measurements of IOPs and radiometry used as model input and in match up assessments, respectively. Best performance was observed when a constant backward phase function shape based on the findings of Sullivan and Twardowski (2009) was assumed in the model. Critically, using a constant phase function in the backward direction eliminates a key unknown, providing a path toward inversion to solve for a and bb. Performance degraded when using other phase function shapes. With available data sets, the model shows stronger performance than current state-of-the-art look-up table (LUT) based BRDF models used to normalize reflectance data, formulated on simpler first order RT approximations between rrs and bb/a or bb/(a + bb) (Morel et al., 2002; Lee et al., 2011). Stronger performance of ZTT relative to LUT-based models is attributed to using a more representative phase function shape, as well as the additional degrees of freedom achieved with several physically meaningful terms in the model. Since the model is fully described with analytical expressions, errors for terms can be individually assessed, and refinements can be readily made without carrying out the gamut of full RT computations required for LUT-based models. The ZTT model is invertible to solve for a and bb from remote sensing reflectance, and inversion approaches are being pursued in ongoing work. The focus here is with development and testing of the in-water forward model, but current ocean color remote sensing approaches to cope with an air-sea interface and atmospheric effects would appear to be transferable. In summary, this new analytical model shows good potential for future ocean color inversion with low bias, well-constrained uncertainties (including the VSF), and explicit terms that can be readily tuned. Emphasis is put on application to the future NASA Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) mission.


Introduction
Radiative transfer (RT) approximations linking inherent optical properties (IOPs), such as spectral absorption a(λ) (m −1 ) and spectral backscattering b b (λ) (m −1 ) to ocean color remote sensing reflectance R rs (λ) are vital to interpreting R rs because it is not possible to analytically invert the full RT equation [1,2]. Once forward RT approximations are developed, inversions can then be explored to solve for IOPs from R rs and subsequently ocean biogeochemical properties [3]. Ocean color R rs (λ) (sr −1 ) here is defined as L w (λ)/E d (λ), or water leaving radiance (W m −2 sr −1 nm −1 ) normalized to above water downwelling irradiance (W m −2 nm −1 ) (see Reference [4] for complete definitions of all optical parameters and Appendix A, Table A1 for notation).
Ocean color expressions to date have almost exclusively relied on first order approximations of RT relating R rs to b b /a through a proportionality represented as f /Q [5][6][7][8], or to b b /(a + b b ) with multi-term polynomial expressions based originally on Gordon et al. [9] with coefficients represented as l or G ( [9][10][11][12]; see review by Werdell et al. [13]). The coefficients describing the relationship between R rs and IOPs are detailed in look-up tables (LUTs), or a neural network in the case of [8], with dependencies on geometry (i.e., solar zenith, viewing angle, relative azimuth) and in some cases wavelength, wind speed, atmospheric conditions, and/or chlorophyll concentration [Chl]. The LUTs are generated from full RT computations with so-called synthetic data sets where IOPs and their interrelationships are assumed and referenced to [Chl]. LUT coefficients are thus also implicitly dependent on IOPs.
These ocean color relationships have been tremendously useful to the ocean color community for decades. Morel [14] discusses the first order analytical relationship and the empiricism necessary to invert. Coefficients describing the relationship between R rs and IOPs are dependent on the bidirectional reflectance distribution function (BRDF), which describes the transformation of downwelling irradiance for different solar zenith angles into the distribution of upwelling radiance. Morel et al. [7] describe the current state-of-the-art BRDF model (herein referred to as M02), currently implemented operationally by the ocean color community (NASA Ocean Biology Processing Group (OBPG); [15]) to convert L w measured at any viewing geometry to a conceptual "exact normalized" L w , [L w ] ex N , with Sun at zenith and nadir viewing in a non-attenuating atmosphere. This allows any measurement at any Sun-viewing geometry to be directly intercompared. Inversion algorithms to derive IOPs are then typically applied to [L w ] ex N and are often based on the same type of first order R rs to IOP approximation [13,16]. The current semi-analytical algorithm (SAA) application to water-leaving radiances is thus a two-step process (after applying calibration coefficients, georeferencing, and atmospheric correction): (1) A BRDF correction to carry out the conversion to [L w ] ex N , followed by (2) application of an inversion algorithm to derive IOPs (see [13]).
A key potential source of uncertainty in current SAA approaches is associated with the volume scattering function (VSF; β(ψ) m −1 sr −1 , where ψ is scattering angle), and this uncertainty is ambiguously dispersed in both steps. The VSF dependency is of paramount importance; as Morel and Gentili [5] note, the BRDF " . . . is essentially controlled by the shape of the VSF . . . ." Historically, very few measurements of the VSF have been available to develop and test SAAs. The M02 BRDF tables were developed from extensive RT modeling using VSF shapes (also known as phase functions P(ψ), defined as the VSF normalized to total scattering, P(ψ) = β(ψ)/b, with units sr −1 ) tied to estimated chlorophyll concentrations [Chl]. To obtain P(ψ) for modeling, Morel et al. [7] first mixed two phase functions for populations dominated by large and small particles representing high and low [Chl] extremes, respectively, and then mixed that particulate phase function with the phase function for molecular seawater to a degree that was also linked to estimated [Chl]. The [Chl] estimate was initially derived from R rs with an empirical band ratio algorithm. Uncertainty is thus associated with how close the assumed P(ψ) used in the BRDF modeling matches the actual P(ψ) associated with any given L w measurement. While the M02 model has been validated in Case 1 waters with atmospheric conditions and IOPs that presumably agree with the underlying atmospheric and bio-optical models [17,18], there have been limited attempts to assess any embedded phase function uncertainties with actual VSF measurements in diverse water types [19,20]. No subsequent BRDF approach has demonstrated enhanced performance for Case 1 waters relative to the M02 approach while being feasible to implement for ocean color. However, Talone et al. [21] have recently shown the Lee et al. [11] LUT-based approach was more accurate for Case 2 waters sampled in the Adriatic, Baltic, and Black Seas.
Uncertainty also arises from exclusively using the b b term in the inversion algorithm, i.e., in step two of the SAA approach described above. Effort has been devoted to trying to correct this uncertainty by layering additional VSF dependence in the f /Q proportionality [22]. Interestingly, after the current M02 BRDF normalization is applied (i.e., measurement geometry transferred to Sun at zenith and nadir viewing), the scattering parameter that should be most closely linked to [L w ] ex N , and thus presumably should provide the lowest associated uncertainties in IOP retrievals, is β(π) if single scattering is assumed, which is typically a good approximation [23]. We are not aware however of any algorithm using β(π) in lieu of b b . The practical but still arbitrary choice of converting measured L w to [L w ] ex N with Sun at zenith and nadir viewing may thus not optimize uncertainties in inversion algorithms. The parameter β(π) is almost completely unknown in the ocean, as there are virtually no direct measurements of β(π) in the literature and typical models of particle scattering with simplified particle shapes do not account for possible particle-particle coherent scattering that may cause significant, but poorly understood, enhancement near β(π) [24]. For example, a recent estimate of b b /β(π) made with a combination of airborne lidar and in situ b b measurements [25] was 50% the value expected from extrapolation of available β measurements [26].
Although the M02 BRDF correction as currently applied can be considered "conceptual" at L w (0,π) without creating a problem from a geometry point of view (i.e., L w at any viewing geometry can be corrected to a conceptual standard L w at another geometry as long as it is consistent), it will be problematic if we try to develop algorithms based on β(π) to reduce uncertainties because the necessary data for algorithm development and validation are lacking. This would also suggest a possible benefit in using the BRDF correction to obtain [L w ] ex N at another unique geometry (solar zenith, viewing angle, and azimuth), one that was representative of single scattering at an angle that we can measure directly with available instrumentation. To minimize the magnitude of BRDF correction, this angle could be chosen as the centroid angle of the maximum in the frequency distribution of in-water single scattering angles observed for polar orbiting satellites, which is about 150 • (Figure 1). Since we can accurately measure β at or near 150 • with commercially available instrumentation (e.g., WET Labs ECO, In-situ Marine Optics IMO-SC6, and legacy HOBI Labs Hydroscat sensors), there would appear to be potential to reduce uncertainties with β(150 • )-based algorithms matched to a β(150 • )-based BRDF correction using M02 relative to the current b b -based algorithms applied to a β(π)-based BRDF correction. In the former approach, uncertainties associated with the phase function are thus mostly restricted to the BRDF correction step.
Herein we explore a different approach, working with a RT expression from Zaneveld [27,28] that explicitly incorporates a dependency on VSF shape, as well as specific viewing geometry. A path to inversion to IOPs is also presented where SAA steps 1 and 2 mentioned above are combined in a single relationship and uncertainties related to VSF shape and other parameters can be directly assessed. This approach has several potential benefits, including (1) a single, fully analytical and invertible expression describing the RT process for all remote sensing geometries, (2) optimal retention of native RT relationships with more degrees of freedom than the first order approximation, directly linked to physically meaningful terms, (3) all parameters including the VSF are explicit in the model with readily characterized uncertainties, (4) model can be readily enhanced by tuning one or more terms rather than developing new LUTs from complete recomputations of full RT, and (5) native viewing geometries produce scattering at angles that can be resolved with available instrumentation. A challenge in an explicitly VSF-dependent approach is amenability to inversion, as some information about the phase function is ostensibly required [1,20]. However, it should be pointed out this is also the case with any model, including the M02 BRDF correction, and subsequent inversion algorithms and their inherent assumptions about VSF shape. In fact, as is assessed herein, the same assumptions of M02 in linking changes in VSF shape with an independent estimate of chlorophyll can be directly applied to this RT model, with the benefit of being able to directly quantify associated errors. In this paper we focus on the performance of a forward implementation of the Zaneveld expression using a modified formulation in terms of IOPs that is now amenable to inversion. Performance of the inversion to IOPs is being fully assessed in ongoing work. A key advance promoting inversion is the finding by Sullivan and Twardowski [26] that the shape of the particulate VSF in the backward direction is relatively constant for a wide range of water types, and so may be represented by a constant function in an algorithm without introducing significant error. Importantly, we have also now overcome a limitation in practically assessing this approach by collecting a database of measured VSFs resolved over a large dynamic range concurrently with other high quality IOPs and radiometry under ideal cloud-free environmental conditions [29].

Forward Model Development
Zaneveld [27,28] derived an exact restatement of the RT equation assuming a plane-parallel, optically deep water column in terms of upwelling radiance L u in viewing direction θ v : where E od is the scalar downwelling irradiance, K Lu is the attenuation coefficient for upwelling radiance, b f is forward scattering, and c is the attenuation coefficient. Parameters f b and f L are light field shape factors representing the path radiance term of the RT equation [27,28]. Factor f b describes the redirection of downwelling radiance into the upwelling viewing angle for a given VSF and is normalized to the redirection that would be observed if the VSF was isotropic. Factor f b is thus directly linked to the shape of the VSF in the backward direction. The factor f L is defined similarly but describes the redirection of upwelling radiance into the upwelling viewing angle and is thus directly linked to the shape of the near forward VSF. These shape factors are expected to vary within a relatively narrow range near unity. All terms are a function of depth z and wavelength λ. The Zaneveld expression does not account for any inelastic processes, such as molecular Raman scattering or fluorescence [30][31][32][33]. Using physically reasoned approximations for f b including the assumption of single scattering, the following was obtained by Zaneveld ([28]; his Equation (14)): where β is the VSF (including water) and ψ is the in-water scattering angle formed between peak incident sunlight at solar zenith angle θ s and scattered light traveling in viewing direction (θ v ,Φ), where cos(ψ) = cos(θ s )cos(θ v ) − sin(θ s )sin(θ v )cos(φ). Azimuth φ is relative to the Sun's direction. Zenith angles θ are in-water, refracted through the air-sea interface and determined from a vertically downward direction. For nadir viewing, ψ = π − θ s and −cos(θ v ) = 1 in the denominator. Note approximations of Equation (2) provided in Zaneveld [28] assumed only nadir viewing while we retain the full BRDF functionality here. In theoretical analyses, Weidemann et al. [34] showed b b retrievals based on the Zaneveld expression had errors as large as −20% and +40%. However, an extraordinarily wide range of VSF shapes was applied in their simulated data, including VSF shapes for specific particulate subcomponents, such as bacteria, minerals, and phytoplankton. These populations were considered separately as quasi-monodispersions with phase function shapes computed from Lorenz-Mie theory, i.e., assuming homogeneous spheres. These phase function shapes thus had large oscillations with respect to angle, structure that we now know is extremely unrealistic for VSFs representative of bulk in situ particle populations (e.g., [26,[35][36][37]). Weidemann et al. also showed with these phase functions that the Zaneveld approximation of the shape factor at nadir viewing, i.e., f b ≈ 2πβ(π − θ s )/b b , had an average error of only 5% (their Figure 11) and this included overcast conditions that would not occur in remote sensing. To our knowledge, the only study attempting to test the Zaneveld model with directly measured data was He et al. [20] with a subset of the NOMAD data set, which included no VSF measurements, where performance for the BRDF component of the model was comparable to the current state-of-the-art [7,11].
Equation (2) demonstrates the direct link between L u (θ s , θ v , φ) and β(ψ). However, for the Zaneveld expression to be a practical tool for ocean color, the terms K Lu and f L must be expressed in terms of IOPs. Furthermore, the IOPs in the model should ideally be coefficients that are closely linked to reflectance and directly measurable with good accuracy using existing sensor technology. For example, the term in the denominator [c − f L (θ s , θ v , φ)b f ] has two IOPs that are difficult to determine directly because of acceptance angle issues with standard transmissometer designs [38,39]. However, since f L is close to unity, this term is also closely related to a + b b , immediately recognizable from commonly used first order approximations. Goals of the next several sections are: 1) To represent Equation (2) entirely in terms of such IOPs, 2) to rework in terms of the commonly used remote sensing reflectance L u /E d , since these are the measurements currently available in validation data sets, and 3) to include the inelastic effects associated with water Raman scattering.

Diffuse Attenuation of Upwelling Radiance K Lu
For the term K Lu , Zaneveld [28] suggested an assumption of equivalency to K ∞ , the diffuse attenuation coefficient in the asymptotic regime. Asymptotic theory is based on the principle that the shape of the light field with depth gradually transforms from being dependent on the incident surface light field to being constant, azimuthally symmetric (so L is only a function of θ), and dependent only on IOPs. Attenuation coefficients for all aspects of the light field, i.e., for all radiances and therefore all irradiances, are equivalent in the asymptotic regime and are also IOPs. This is a critical assumption for the purposes of model development since K can then be described only in terms of IOPs. Zaneveld justified the assumption, even though the omnidirectional light field in surface waters is far from asymptotic, due to the decoupling of upwelling radiances to downwelling radiance distributions. Measurements of upwelling radiance fields from the 1960s and 1970s showed a near constant shape and attenuation rate with depth. Recently, Twardowski and Tonizzo [40] confirmed this assumption in RT simulations with no more than 3% error when the sun was at solar zenith and the following relationship held for above water solar zeniths θ s up to 75 • : which included a full range of possible natural water types. Coefficients f are provided in Appendix A, Table A2. Angle θ s is related to θ s by Snell's Law, i.e., θ s = sin −1 (1.34 sin(θ s )). Only nadir viewing was considered in the simulations for K Lu in Equation (3), and we make the assumption that K Lu for other viewing angles that define a specific scattering angle ψ can be approximated by Equation (3) with nadir viewing geometry after assigning a θ s that provides an equivalent in-water ψ. For example, for θ s = 60 • , θ s will be 40.3 • and ψ will be 139.7 • for nadir viewing; we assume the resulting F(60 • ) will also be applicable to any off-nadir viewing angle with ψ = 139.7 • . This assumption has been verified using Hydrolight RT simulations (methodology addressed in Section 3.2) to no worse than 2% in the solar plane and no worse than 5% within the upwelling hemisphere for in-water scattering angles consistent with remote sensing (i.e., Figure 1). This assumption implies a rotational reference frame, where the first order determinant of radiance field shape in the model, i.e., ψ, is preserved. Two potential drawbacks of this assumption are (1) the influence of skylight may be skewed in the rotated reference frame, and (2) the range of viewing angles is restricted since the smallest ψ is~134 • for underwater nadir viewing. The range ψ > 134 • , however, comprises >95% of the expected scattering angles that will be measured by the PACE imager ( Figure 1). Further work with field measurements is needed to verify this assumption. Reformulating Equation (3) in terms of dependency on the in-water scattering angle, we obtain: Coefficients f A are provided in Appendix A, Table A2. From Gershun's Law we can set K ∞ = a/µ ∞ , where µ ∞ is the average cosine of the asymptotic light field. After inserting Equation (4) into Equation (2) and allowing for c = a + b, the following is obtained: where The parameter Ψ K Lu is assumed spectrally independent with errors over the full range of possible water types estimated at <2% [40]. Dividing numerator and denominator of Equation (5) by b b , we obtain: where b b = b b /b is the backscattering ratio. For IOPs we now have the backward phase function in the numerator; in the denominator we have the recognizable b b /a, as well as b b , an IOP that incorporates information on bulk particle composition (see Section 2.4; [41]), and is not typically associated with ocean color remote sensing vis-à-vis the common first order approximation of r rs assumed proportional to b b /a.

Average Cosine of the Asymptotic Light Field µ ∞
The µ ∞ term in Equation (7) must be expressed in terms of IOPs to invert. Using a fit to theoretical calculations of radiance fields by Prieur and Morel [42], Zaneveld [28] recommended 1/µ ∞ be modeled empirically with respect to the single scattering albedo ω (=b/c) using a quadratic fit. Additionally, through more detailed RT computations, Berwald et al. [43] found a 4th order dependency of µ ∞ on the albedo.
In the study by Twardowski and Tonizzo [40], µ ∞ was parameterized in terms of b b /a instead of ω, including an assessment across a full range of environmentally representative phase function shapes using the Fournier-Forand analytical model [44,45]. The data set used in the assessment was not a representative synthetic data set (i.e., [46]), as it included a full range of possible b b /a values, possible phase functions, and permutations thereof. Representing µ ∞ in terms of b b /a has two distinct advantages. First, b b and a can be measured with commercially available in situ instrumentation with accuracies of a few percent [47][48][49] to enable performance assessment for algorithms. Parameters c and b cannot be measured without significantly larger errors, typically >25-50%, because of the acceptance of near forward scattered light in conventional transmissometer designs [38,39]. Secondly, b and c are not parameters that are closely linked to r rs without additional information, whereas b b and a are (e.g., Gordon et al. [50]), and a goal here is to represent the entire model in terms of b b and a to enable inversion.
The Twardowski and Tonizzo [40] parameterization also explicitly depended on η bb , the fraction of b b attributable to molecular scattering, η bb = b bw /(b bp + b bw ) [6]. The natural range for η bb is from 0 to~0.98 [40,51], and the range used for b b /a was 10 −4 to 10 −1 . Since the water components of η bb may be assumed known [52], the effective unknown here is b bp . After extending the analysis from Reference [40] to include near zero b b /a and increased resolution in η bb , the resulting fit was obtained for µ ∞ : Coefficients m are provided in Appendix A, Table A2. Fits to simulated data were again made with the polyfit function from MATLAB. Absolute errors %δ abs for this fit vary from 0.19 to 3.5 for η bb ranging from 0.98 to 0.0098, respectively. Since both b b /a and η bb are spectrally dependent, µ ∞ will be as well (not shown for clarity).

Backward Phase Function β(ψ)/b b
The scattering parameters in Equation (7) must be expanded into water and particulate components. Expanding β(ψ)/b b gives: where the p and w subscripts represent particles and molecular seawater, respectively. In Equation (9), the pure seawater terms, which are temperature and salinity specific, can be directly computed with an estimated error of no more than 2% [52]. Introducing a term for the particulate phase function in the backward direction, P bb (ψ) = β p (ψ)/b bp , Equation (9) then becomes: Note b bp and b bw both have spectral dependencies as does β/b b , and the unknowns are b bp (λ) and P bb (ψ). Equation (10) can be easily rewritten in terms of η bb , with the same unknowns. In coastal waters where η bb is near zero, β p >> β w , b bp >> b bw , and β(ψ)/b b will be approximated by P bb (ψ). For clear ocean waters, phase functions are represented by a mixture of both particles and pure seawater with VSF shapes dependent on η bb .

Backscattering Ratio b b
Expanding the backscattering ratio where b bp is the particulate backscattering ratio. This b bp is the "true" b bp , distinct from the b bp typically derived from measurements that include c data from transmissometers with significant acceptance angle errors (e.g., [38,39]) (we note models linking particle biogeochemical properties and measured b bp should account for the acceptance angle of c measurements [41,53]). Bootstrapping exercises using Equation (7) readily show an impact on reflectance of up to several percent when b bp is varied over the full~0.003 to~0.03 dynamic range observed in the oceanic environment [41,[53][54][55][56]. All terms have spectral dependencies.

Shape Factor f L
Zaneveld [28] recommended the term f L , i.e., the dimensionless upwelling radiance shape factor, could be set to 1.05 with small error. The natural range was estimated at 1 to 1.12 [27]. In the Results, we develop a new model for this term.

Remote Sensing Reflectance Formulation
Nearly all ocean RT algorithm work over the last several decades has used reflectance defined as L u /E d instead of L u /E od . Zaneveld [28], however, pointed out, as is evident from Equation (1), L u /E od is most closely aligned with RT theory. The irradiance parameter E od is also less dependent on solar zenith angle than E d . Furthermore, sensor technology has been available to measure E od , historically from companies Biospherical (www.biospherical.com), Satlantic (www.satlantic.com), and Trios (www.trio.de). Nonetheless, since nearly all field data over the last several decades has focused on E d , any testing and validation of the RT model described here requires modification in terms of E d . Substituting E d /E od = µ d (the average cosine of downwelling radiance, just below the air-water interface) into Equation (7) gives: where r rs is the classically known remote sensing reflectance just below the water surface. Full dependencies of variables not shown for clarity.

Average Cosine of the Downwelling Light Field µ d
A model for µ d is now required in the expression from Equation (12). If the sky is ignored and we assume a negligible fraction of the incident solar beam is scattered in the near-surface, a reasonable first-order approximation of µ d should be µ w ≡ cos(θ s ). Adding in a cardioidal radiance distribution for skylight, Morel and Prieur [57] obtained: and noted that for sun angles between 8 • and 62 • , µ d varied only from 0.79 to 0.94. Thus, even without knowledge of solar zenith angle, a median value could be used with an accuracy better than 10%. Numerator values of 0.6 and 0.4 in Equation (13) represent the fractions of direct (E dd /E d ) and diffuse (H = E ds /E d ) downwelling light, respectively, which together equal unity. These values primarily depend on θ s ' and horizontal visibility V, the latter of which depends on aerosol optical thickness (AOT). For 20 • ≤ θ s ' ≤ 60 • and for H between 0.2 and 0.5, the variability of µ d is~7% for a fixed θ s '.
The term µ d can be factorized in two parts, one part dependent on the IOPs, the other dependent on the atmospheric conditions and geometry: The atmospheric component can be represented as: As mentioned, Morel and Prieur [57] assumed H = 0.4. Gregg and Carder [58] later provided a relationship for H as a function of θ s and V that included skylight. Specifically, the results of ( [58]; their Figure 4) can be fit as follows: Coefficients e are provided in Appendix A, Table A2. Note the typical default value used in Hydrolight is V = 15 km.
The P 3 term in the M + d relationship is a 3rd order polynomial in cos(θ s ) to correct where Morel and Prieur's original approximation deviates from the Gregg and Carder relationship at large θ s . The P 3 term is: The IOP-dependent component M * d was modeled using the approach in [40] with the extended b b /a range and η bb resolution discussed in Section 2.2. The final fitted relationship is: Coefficients m * d are provided in Appendix A, Table A2. Absolute percent error with this relationship relative to computations using full radiative transfer (Hydrolight, see Section 3.2) varies from 0.06% to 0.4% across all θ s '. Error for the full µ d expression is <1%. Spectral dependency enters the expression through the V term, which is dependent on spectral AOT.

Including Inelastic Water Raman Effects
The molecular water Raman scattering contribution to r rs can be included in Equation (12) as an additive term, resulting in the expression: Full dependencies for all parameters in the model are shown except for λ; all parameters exhibit dependence on λ in the model except for Ψ K Lu (ψ). We note a similar approach has been taken in adding water Raman effects in other reflectance models (e.g., [59]; reviewed in [13]). The term r rs,Raman can be derived according to Westberry et al. ([60]; see their Equation (7)) with inputs of above water (z = 0 + ) downwelling irradiance E d (0 + ,θ s '), a, and b b . The NASA Generalized IOP (GIOP) inversion model implementation [16] also currently uses this Raman formulation. Terms for other inelastic effects, such as fluorescence from dissolved organic matter and pigments may also be added if representative models are available.

ZTT Model Summary
Equation (18) is the final model for ocean color reflectance defined as L u /E d , called the ZTT model hereafter. In Equation (18) the µ d term is approximated by Equation (14), the β(ψ)/b b term approximated by Equation (10), the Ψ K Lu term described by Equations (4) and (6), the µ ∞ term described by Equation (8), and b b represented by Equation (11). The ultimate assignment of the f L term is addressed in the Results. The geometry variables, E d (0 + ,θ s '), V, and molecular water scattering parameters in the above can be considered knowns. The model is fully spectral. Key unknowns are b b and a (or b b /a and η bb ), and since pure seawater absorption a w in the visible is considered known with good accuracy [61] the effective unknowns are b bp and absorption by non-water constituents a pg . The two additional unknowns are P bb (ψ) and b bp . In the forward implementation here, these four parameters must be provided from direct measurements or through some assumptions. In the inversion implementation, these are the parameters that may be solved through techniques to minimize errors in the expression if there are enough spectral bands (i.e., degrees of freedom) in r rs , although a priori assumptions may be required for P bb (ψ) and b bp .
For ocean color reflectance defined as L u /E od , the simpler relationship in Equation (7) can be used. Water Raman effects can still be considered by applying an algorithm, such as in Reference [60].

Synthetic Dataset
A synthetic data set referenced to [Chl] was developed to test the ZTT model and develop an expression for the f L term (see Section 4.1). Twenty values of [Chl] were assumed, logarithmically spaced between 0.01 and 30 mg m −3 . Total absorption coefficient (a) was represented by a sum of four components: Pure water absorption a w was taken from Reference [61]. Phytoplankton absorption a ph was calculated from chlorophyll concentration [Chl] and from spectrally averaged absorption coefficients of micro-and pico-plankton (a micro (λ) and a pico (λ), respectively) [62]: where S f = [0.25, 0.5, 0.75] is the shape mixing factor. Non-algal particulate absorption a d was given by [63]: where a d (440) is equal to R d a ph (440), with R d = [0.05, 0.1, 0.5, 1, 2]. Similarly, we used an exponentially decaying expression for the absorption of chromophoric dissolved organic matter, a g [64]: where a g (440) is equal to R g a ph (440), with R g = [0.1, 0.3, 0.5, 1, 2].
Total backscattering b b was represented by the sum of water and particulate components: To derive particulate backscattering b bp , total particulate scattering b p was first empirically estimated at 550 nm [65]: and then extrapolated spectrally [66]: where Pure seawater backscattering b bw was calculated according to Zhang et al. [52].

Radiative Transfer Simulations
RT simulations of r rs were performed with Hydrolight (Sequoia Scientific, Bellevue, WA), following the procedure in Tonizzo et al. [29]. Fournier-Forand analytical phase functions were derived for each [Chl] iteration following the method of Mobley et al. [67].
Inelastic water Raman scattering was not included in the simulations for the synthetic data set, as this is separately addressed in the model (Equation (18)); fluorescence from any seawater constituents was also not considered. Note water Raman effects were included in full RT simulations for all field data (see [29]). Output wavelengths ranged from 350 to 800 nm at 5 nm resolution. Hydrolight default atmospheric parameters were used for all runs. Computations were run for θ s of 10, 30 and 60 • . Altogether, 1500 different IOP permutations based on [Chl] were simulated for each θ s '.

Field Data Sets
Model performance was assessed using two aggregate data sets. The first is a high quality data set of 23 stations from the Ligurian Sea, waters around the Marine Optics BuoY (MOBY) west of Lanai, the southern California coast, and the New York Bight, collected in 2008 and 2009 as part of NASA Spectral Ocean Radiance Transfer Investigation and Experiment (SORTIE) and NASA Ocean Color Validation (OCVAL) exercises. Radiometric closure was assessed in detail for these data by Tonizzo et al. [29].
[Chl] ranged from 0.24 to 23.85 mg m −3 . The VSF was directly measured using the custom Multi-Angle SCattering Optical Tool (MASCOT) [26,36] and radiometric uncertainties were rigorously assessed by Voss et al. [68]. Here we use WET Labs ac-9 absorption measurements corrected for scattering using independent VSF measurements (i.e., the VSF98P correction [29]; also see Stockley et al. [49] for detailed evaluation of this correction). The second data set is the NASA NOMAD database, where 80 data records were identified containing the parameters needed for performance assessments here, i.e., a, b b , b bp , [Chl], and r rs .
[Chl] ranged from 0.23 to 10.68 mg m −3 in this data set.

Depth Weighting IOPs
Depth-weighted contributions of IOPs to water-leaving radiance were derived after Zaneveld et al. [69] using a 2-stream, first derivative approximation. The approach finds the IOP for a conceptual homogeneous ocean that reproduces the water-leaving radiance observed in a specific stratified ocean case. The remote sensing depth weighting for generic IOP X is: The exponential term can also be found in Gordon and Clark [70], where it is mentioned the 2K(z) term should actually be K u + K d , but approximating with K d alone was not expected to lead to appreciable errors. Several approximations of K d in terms of IOPs are available from the literature (e.g., [71,72]). The approximation by Lee et al. [73] for averaged K d within the euphotic zone is used here and conveniently represented in terms of a, b b and above water solar zenith θ s : For evaluating RT approximations, all IOP terms should ideally be weighted together. For example, Zaneveld et al. [69] considered the simple R ∝ b b /a approximation and demonstrated that <b b >/<a> was not equivalent to <b b /a> for an IOP stratified ocean. However, b b and a are considered independently in the model here and an ultimate objective is to invert the approximation to derive b b and a independently, so each IOP was individually depth weighted. For all stations sampled in the validation data set, the difference between surface b b and depth weighted b b was <1%; with a, this was also the case for most stations, but reached a difference of 5% for one station sampled.

Metrics for Error Assessment
Mean absolute percent error (MAPE), %δ abs , is a commonly used metric in assessing performance in r rs match ups: The MAPE metric takes into account the absolute magnitude of the residuals, giving them equal weight. Other commonly used metrics include root mean square error (RMSE), a measure of accuracy and potential forecasting errors in simulating r rs when the errors may be assumed to be unbiased and normally distributed [74]. RMSE gives greater weight to larger errors than δ abs . Since bias errors are often expected to be more significant than random, normally distributed errors in simulations, %δ abs is expected to be the most appropriate metric to assess match ups [29]. Other common metrics (e.g., [75]) are shown in the following plots as appropriate. Also see Werdell et al. [13] for a detailed discussion of performance metrics.

Developing an Expression for the f L Term
Dynamics of f L were explored in the synthetic data set. As found by Hoge et al. [76], dependencies on wavelength and θ s were most important, and an average spectral shape f L,ave (λ) scaled to solar zenith according to sin(θ s ') was ultimately found to be most representative ( Figure 2): Function f L,ave (λ) is provided in Appendix A, Table A3. Equation (30) was developed by minimizing residual errors between results from full radiative transfer and the ZTT model for each wavelength and each solar zenith angle for the full synthetic data set. Spectral shapes for f L were relatively consistent in the results, so the term f L,ave (λ) was derived from an average for these results after spectral normalization. The scaling factor in brackets was a suitable function to minimize errors with respect to solar zenith angle. Only nadir viewing was again considered in the simulations for f L , and we again make the assumption from Section 2.1 that f L for other viewing angles that define a specific in-water scattering angle ψ can be approximated by f L observed at the nadir viewing geometry with equivalent ψ. The potential caveats mentioned in Section 2.1 apply here as well. Restating Equation (30) in terms of ψ, the following is obtained: Similarities between f L,ave (λ) and a typical ocean absorption spectrum are noted (Figure 2), resulting from the effects of multiple scattering, where higher relative scattering (lower relative absorption) promotes a radiance field closer to that which would result from isotropic scattering, i.e., where f L approaches unity.    Figure 3A shows results deriving r rs using the ZTT model for the synthetic data set, applying the approximation from Equation (31) for f L and using all synthetic data IOPs as input. This match up is, thus, effectively assessing combined errors in the ZTT model from the approximations of f L , f b , K Lu , and µ d . We note f L was specifically optimized to this data set, so in this respect Figure 3A shows what may be considered a best-case matchup. MAPE %δ abs for the full data set was relatively small at 2.68%. The largest errors were observed in three r rs spectra, with the model showing overestimation bias centered around 570 nm. These spectra were associated with the highest [Chl] modeled in the data set and lowest b bp . Colors represent wavelengths, from 400 to 700 nm with gray used for 350 < λ < 400 nm and 700 < λ < 800 nm. [26] found a high degree of consistency in the shape of the particulate VSF in the backward direction, i.e., β p (ψ)/b bp , for a large data set covering a wide dynamic range in bulk particle composition. Deviations from this relationship were no more than 5% through the entire backward angular range. While it is well known the phase function shape over the full angular range varies substantially (e.g., [37,41,45,77]), only the shape in the backward direction that is important for remote sensing and the ZTT model is considered here.

Sullivan and Twardowski
If we assume β p (ψ)/b bp is a constant shape P bb,ST (ψ) after Sullivan and Twardowski [26] in the ZTT model, replacing P bb,FF (ψ, b bp ) from the synthetic data set, errors increase by only~0.3% ( Figure 3B). We assume P bb,ST (ψ) is constant spectrally after [41] and others. This result is significant, as it demonstrates the potential for eliminating one of the key unknowns in the ZTT model. We note that an averaged Fournier-Forand phase function in the backward direction could be used in the model, which would enhance the analytical character of the model, but the Sullivan and Twardowski [26] averaged phase function, derived from extensive measurements in a wide range of water types, may be more representative of the natural environment. This is assessed further in the next section. Figure 3C shows results of setting f L to a constant 1.05 after [28] and using a constant P bb,ST (ψ). MAPE increases significantly to 5.85%, showing the importance of using the f L model in the ZTT. Each of the three distinct data groupings along the 1:1 relationship are associated with the individual solar zeniths that were used, i.e., 10, 30, and 60 • , highlighting the influence of the f L model in removing the effects of solar zenith. All ZTT model runs hereafter use the f L model from Equation (31). The constant backward phase function P bb,ST (ψ) is used unless another phase function is specified. Figure 4A shows results for the ZTT model, with inputs of constant P bb,ST (ψ) and direct measurements of a pg , b bp , and b bp , compared to measured r rs for the validation data set [29]. MAPE %δ abs was 16% for all λ (summarized in Table 1; spectral %δ abs provided in Table 2). Interestingly, this result is slightly more favorable than r rs computed with the full radiative transfer (i.e., Hydrolight) using measured phase functions, where %δ abs was 17% for the same data (as reported in Tonizzo et al. [29]). Based on closure analyses for these data [29], this 17% error represented the aggregate inherent error from all sources within the data and computations, i.e., from the IOP and radiometric measurements, as well as any errors from the assumptions within the Hydrolight RT code. For the ZTT model evaluated using measured phase functions instead of the constant P bb,ST (ψ), i.e., the same approach followed with the full RT computations, %δ abs was also 17% ( Table 1). The different backward phase functions applied in the ZTT model are shown in Figure 5. Measured backward phase functions showed more variability at individual angles than the average P bb,ST (ψ) derived from a large data set of diverse Case 1 and Case 2 water types, likely the result of small scale hydrosol patchiness along the different light paths for individual measurements at each angle. This observation is addressed further in Section 5. Moreover, the strong agreement in absolute error between full radiative transfer simulations and the ZTT analytical approximation for this diverse data set is encouraging.

Assessment with High Quality Validation Data
Excluding the effects of water Raman from the ZTT model increased absolute error by 3%, demonstrating Raman had a significant effect, especially >580 nm (results not shown). Figure 4B shows ZTT model performance with a further constraint, setting b bp constant at 0.006, so the only IOP inputs were measured a pg and b bp . With respect to inversion, this approach has the same unknowns, i.e., a pg and b bp , as current algorithms based on the simple first order approximation [16]. Resulting MAPE %δ abs was 17%, comparable to full RT computations. ZTT runs with other fixed b bp are shown in Table 1. Some influence from b bp is apparent, which has not been fully appreciated in previous models based on the first order approximation of r rs to b b /a. ZTT model runs with Fournier-Forand phase functions for the particulate component, computed from measured b bp following the method of [67], resulted in %δ abs of 19%, 3% higher than results using the constant P bb,ST (ψ) ( Table 1) [29]. Note L11 does not account for Raman scattering. Chlorophyll input for M02 was derived from spectral absorption using the Nardelli and Twardowski [78] line height algorithm. MAPE %δ abs was 16%, 17%, 21%, and 21%, in A through D, respectively. defined as f /Q, which is a function of viewing geometry, wind speed, and [Chl].
[Chl] can be estimated from an empirical r rs band ratio algorithm. For the assessment here, [Chl] input was either quantified directly from discrete samples or derived from spectral particulate absorption measurements using the line height method with specific absorption a p *(676) of 0.0108 m 2 mg −1 after [78]. Note that bootstrapping exercises have shown relatively large errors in this estimated [Chl], i.e., up to +/−50%, affect %δ abs in match ups by typically no more than 1%. The other inputs were measured a pg and b bp . Match up %δ abs for M02 was 21% for the data set from [29] (Table 1). Notably, M02 performed well in the blue spectral region (Table 2). To derive phase functions, Morel et al. [7] used randomly-oriented spheroidal particles, computing the scattering phase function with the T-matrix method. Because the computations for spheroids were lengthy, a single refractive index n p of 1.06 was considered with sizes ranging from 0.02 to 14 µm. A Junge-type particle size distribution was assumed, i.e., a power law model, with exponents of 3.1 and 4.2 representing two end-member populations. Particulate scattering phase functions were then calculated from (see Figure 5): If the Morel et al. [7] approach for deriving the phase function is used as input into the ZTT model, replacing the constant shape P bb,ST (ψ), MAPE %δ abs increases significantly from 16% to 23% (Table 1).
Finally, match ups in r rs were assessed in a recently published BRDF model from Lee et al. ( [11]; called L11 hereafter) that, like M02, developed a LUT based on RT simulations with a [Chl]-referenced synthetic IOP data set ( Figure 4D; Tables 1 and 2). L11 uses a quadratic form of b b /(a + b b ) split into particulate and molecular seawater components with each term scaled by a G coefficient, with solutions provided in LUTs. In the synthetic data set, Lee et al. [11] used an averaged Petzold phase function and Fournier-Forand phase function assuming b bp = 1% as particulate phase function endmembers. L11 however does not include water Raman effects. Like M02, match up MAPE %δ abs was also 21% for all λ.    [26] with the assumption P bb,ST (ψ > 170 • ) = P bb,ST (170 • ); "Morel" is from Morel et al. [7] with endmember populations dominated by small and large particles also shown; "Measured" were directly measured; and "FF" are Fournier-Forand phase functions P bb,FF (ψ, b bp ) derived from measured b bp following [67]. Figure 6A shows results for the ZTT model applied to a subset of 80 measurements from the NASA NOMAD global data set. IOPs were averaged over the first optical depth. Match up MAPE %δ abs of 20% was observed for this global data set, comparing favorably to the error of 16% observed with the high quality data set from Tonizzo et al. [29] (Table 1; see Table 2 for spectral breakdown). This is especially true considering the different IOP and r rs processing approaches (e.g., there are several scattering error correction options for in-water absorption measurements) used for these data sets. If b bp is constrained to 0.006, MAPE %δ abs increases to 23% ( Figure 6B). Results from ZTT runs for other fixed b bp are shown in Table 1, where the dependency of performance on b bp is again apparent.

Assessment with Global NOMAD Data Set
If the Morel et al. [7] approach for deriving the phase function from [Chl] is used as input into the ZTT model, %δ abs increased from 20% to 26% for the NOMAD data (Table 1). If measured b bp are used to derive P bb,FF (ψ, b bp ) following [67], %δ abs increased to 25% for the ZTT model. Table 2. Summary MAPE %δ abs results for calculated r rs vs measured r rs for each λ, for SORTIE and OCVAL data sets (see Figure 4) and NOMAD data (see Figure 6). Match ups are also shown for implementation of the full M02 and L11 models with the NOMAD data set ( Figure 6C,D). MAPE statistics are summarized in Tables 1 and 2. There are several possibilities for the systematically lower measured r rs in the red spectral region for these data: (1) Suboptimal surface extrapolation of radiances due to water Raman effects [80,81], (2) residual scattering errors in reported in-water absorption measurements that may promote a bias to higher values of modeled r rs (e.g., [29,49]), and/or (3) suboptimal corrections for sensor self-shading, which can be strong in clear water at longer wavelengths due to water absorption [82,83].    Figure 4, but for the NASA NOMAD data set containing a, b b , r rs , and [Chl]. MAPE %δ abs were 20%, 23%, 25%, and 26% in (A-D), respectively.

Discussion
The ZTT analytical model for remote sensing reflectance is based on the restatement of the RT equation by Zaneveld [27,28] with the following refinements: (1) Full bidirectionality was retained in the final model (no assumption of nadir viewing), (2) the assumption of asymptotic diffuse attenuation for upwelling radiances at the surface was analytically described, (3) the entire model was formulated in terms of b b and a to enable inversion, (4) the VSF was explicitly included with β p (ψ)/b bp set to a constant shape defined by P bb,ST (ψ) [26], which showed the strongest performance of any phase function, including measured phase functions in disparate, diverse data sets, (5) water Raman effects were added, and (6) for r rs , an analytical model for µ d was developed that included dependency on atmospheric visibility. IOP inputs to the ZTT model are b bp , a pg , β p (ψ)/b bp , and b bp . Moreover, the ZTT model enables systematic assessment of the effect of IOPs on the BRDF for the first time.
ZTT showed strong performance in match ups with field data sets spanning a wide range of water types. In terms of MAPE, the model with constant phase function shape P bb,ST (ψ) performed as well as full RT computations. Stronger performance was observed for these available data sets relative to current state-of-the-art LUT-configured models based on the first order approximation between r rs and b b /a or b b /(a + b b ), i.e., M02 and L11.
MAPE for the ZTT model was slightly better using P bb,ST (ψ) (%δ abs of 16%) when compared to using directly measured backward phase functions for the validation data set (where %δ abs was 17%). It may be argued, under certain conditions, the broad average P bb,ST (ψ) may be more accurate than an isolated measurement because specific bias errors in any single VSF measurement at an individual angle may be avoided. This of course assumes such biases may be larger than the natural variability in the shape of the backward phase function. Another explanation is error assessments for the validation data set were made over a relatively few number of stations, 23. To roughly estimate uncertainty in this error metric, single stations were sequentially left out of the error determination for ZTT using P bb,ST (ψ), which gave a total range in %δ abs of 14.8 to 16.8% with standard deviation of 0.6%. A~1% difference in error may thus not be particularly meaningful.
The ZTT model is fully analytical in the sense it is described entirely by equations and "almost fully analytical" with respect to the inclusion of empiricism. Empiricism comes from optimizing the f L parameter to our synthetic data set based on [Chl]-based bio-optical models, with MAPE improving by a significant~3% compared to setting f L equal to a spectrally constant 1.05. Using the constant P bb,ST (ψ) backward phase function shape also introduced empiricism. If sample-specific Fournier-Forand analytical phase functions, based on b bp were used, errors in the validation data set increased by 3-5% (Table 1).
Particulate phase function shapes from Morel et al. [7] were modeled using the T-matrix approach, assuming homogeneous spheroidal shapes with constant n p of 1.06. Although the endmember phase function shapes were not consistent with P bb,ST (ψ) and P bb,FF (ψ, b bp ) shapes, many of the blended phase functions were. If P bb,ST (ψ) may be considered widely representative, increasing errors may be expected in applying the M02 model to the lowest and highest Chl conditions where the endmember phase functions would fully manifest. Phase functions in the Chl range 1 to 10 mg m −3 were in closest agreement to P bb,ST (ψ). Applying phase functions derived using the Morel et al. [7] approach in the ZTT model increased errors relative to using a constant P bb,ST (ψ) by a significant 7% for the validation data set of [29] (from 16% to 23%) ( Table 1). This error was also larger than the 21% error observed for M02 (addressed further below).
Since L11 used a Fournier-Forand phase function where b bp = 1% and the average Petzold phase function as endmembers [shown in 26, their Figure 3], many blended shapes from L11 should be relatively similar in the backward to the P bb,ST (ψ) used in ZTT (although the Petzold shape does exhibit systematic discrepancies). This particular P bb,FF with b bp = 1% is actually a very close match to P bb,ST . L11 was developed through least-squares fitting with respect to Hydrolight results, so all the blended phase functions are weighted in the results.
A "typical marine atmosphere" and associated incident sky radiance distribution were used in the Hydrolight RT computations for the synthetic data set, the effects of which enter the ZTT model through spectral optimization of the upwelling radiance shape factor f L (Section 4.1). The parameters Ψ K Lu and µ ∞ were also solved with Hydrolight RT simulations assuming the same atmosphere. The M02 and L11 models also assumed a single "typical marine atmosphere" in their simulations using Hydrolight based on Reference [58], so this aspect of the models should be consistent. The µ d term in the ZTT model has added flexibility in accommodating varying atmospheric visibility (closely related to AOT; [58]) inasmuch as V is an explicit term for µ d .
One question arising from this work is, if the phase function shape P bb,ST may indeed be considered highly representative, how would a LUT-based model following the approach of M02 but using P bb,ST perform? To test this, we developed a LUT based on the first order r rs proportionality to b b /a, using the synthetic data set described in Section 3.1. The particle phase function was constant, consisting of P bb,ST in the backward and a Fournier-Forand phase function with b bp = 0.006 for the forward direction (this had the best performance for the ZTT model in the validation data set [29], with MAPE = 17%). The LUT approach resulted in a MAPE of 22% for the validation data set and 26% for the NOMAD data. This performance is very close to M02 (21% and 25%, respectively) and L11 (21% and 26%, respectively). An insight we can take from this analysis is that a component, such as phase function shape in the different models cannot be considered in a vacuum, i.e., each model has a large set of assumptions that ultimately influence final performance. Furthermore, although not always explicitly stated in the literature, field validation efforts are likely to have had a role in updating assumptions through the course of development and testing of some models to ensure the model was well aligned with available data. For example, in our case, the f L expression (Section 4.1) was developed from an assumed synthetic data set and the b bp value of 0.006 was recommended as it ultimately provided the best results with the field validation data sets available.
Recent work by Zhang et al. [84] has suggested the possibility of substantial variability in phase function shape in the backward direction, including up to a 40% increase for β p (ψ)/b bp in the near backward relative to β p (120 • )/b bp . These results are in stark contrast to the comprehensive study of [26] and other published phase functions in the literature (e.g., [35,85]). Other recent works showing increases in the VSF at ψ > 150 • [86][87][88] were made in enclosed cuvettes and did not include any correction for internal reflections (e.g., [89]), which are typically a difficult problem in laboratory VSF devices, since the scattered signal in the backward direction is orders of magnitude smaller than the forward. The Zhang et al. [84] results are also inconsistent with the phase functions used in previous BRDF models, such as M02, which have been shown to have satisfactory performance in validation efforts [17,18], even in many Case II coastal waters [19,90]. The recent theoretical r rs modeling study of Xiong et al. [91] based on the phase function shapes presented in [84], where up to 50% disagreement with L11 was observed, consequently appear unrealistic.
Calibration and possible bias errors in the VSF device MASCOT used in [26] to derive P bb,ST (ψ) have been evaluated in detail [36,48]. The MASCOT is an open path, in situ device with 17 independently calibrated detectors resolving the VSF from 10 • to 170 • in 10 • increments. There are two calibration coefficients for each channel, a dark offset (obtained with the source occluded) and a scaling factor to convert digital counts to VSF units m −1 sr −1 . Moreover, the observation of a remarkably consistent shape in the backward phase function for VSFs spanning more than four orders of magnitude in dynamic range cannot be explained by any known bias error in these coefficients. In fact this observation argues against any systematic bias, as consistency in phase function shapes during serial particle suspension experiments is a useful method to assess the accuracy of calibrations.
Recently, He et al. [20] developed a BRDF model also based on Zaneveld [28] but it differs from our implementation is several ways: (1) only the BRDF aspect was considered; (2) one wavelength was considered; (3) the formulation was in terms of a LUT (similar to M02 and L11) for aggregated terms; (4) widely varying backward phase function shapes were assumed after [84] described above; (5) the µ d model depended only on "typical average" sky conditions; and (6) water Raman effects were excluded. Direct comparisons were thus not possible here.

Assessing Residual Bias in the Model
The assumption of single scattering in the simplification of f b (Section 2, Equation (2)) is not expected to invoke significant errors at turbidity levels found for most oceanic and coastal waters [5,13]. As shown by Gordon [23], secondary scattering events for light redirected toward the backward direction, i.e., downwelling radiance scattered into the upwelling radiance field, has a negligible effect on the shape of the upwelling radiance distribution because the vast majority of these scattering events are in the near forward direction. The full effects of multiple scattering were included in relationships for other terms in the model as they were determined using full RT modeling.
MASCOT VSF measurements used to determine P bb,ST (ψ) only extended to 170 • , so the application of this function for ψ > 170 • in the model is uncertain. There were no ψ > 170 • in the data sets tested here, although such scattering angles may be expected for some viewing geometries for the PACE imager ( Figure 1). There is a need for a better understanding of phase function shapes for ψ > 170 • and this extends to other applications, such as lidar [92].
Only nadir viewing was considered for the analytical relationships for K Lu and f L and we suggest these relationships may be transferable to other viewing geometries as long as ψ is preserved. This assumption has been verified in Hydrolight simulations to have small errors for representative cases, but cannot be assessed with the field data sets available here, which all have nadir viewing. This is being considered in future work, and, moreover, emphasizes the need for routine measurements of full upwelling radiance distributions in ocean color validation work as validation only with nadir radiance L u (θ s ,θ v = π) neglects the BRDF effects that are essential to any algorithm.
The optimization of f L was based on the [Chl]-based synthetic data set. A significant deviation from the bio-optical model assumptions of that synthetic data set could lead to a discrepancy in match ups, however it worth pointing out (1) f L has a relatively weak dependency on IOPs except for the VSF shape in the backward direction (which is effectively accounted for in the solar zenith dependency in Equation (30)), and (2) other state-of-the-art models (i.e., M02 and L11) are constructed entirely around such synthetic data sets. A benefit of the ZTT model is that an analytical relationship for a specific term of the model may be replaced later if a more optimal approach is demonstrated.
As shown in Zaneveld [27,28], L u normalized to E od is more closely linked to the RT Equation than L u /E d , where an additional µ d term must be added to the model (Section 2.6). Adding this term is expected to increase uncertainty in the RT model and inversion. Scalar irradiance has the additional advantage of being less dependent on solar zenith angle and sky radiance distributions. As discussed in Section 2.6, sensor technology is available to measure E od and inclusion of these measurements in field efforts focused on algorithm development and validation should be encouraged. When the µ d term must be included in association with E d measurements, AOT is also a useful validation measurement that may be made in the field [93] and used to derive visibility V [58] for input in Equation (16) or full RT computations.
Lastly, potential uncertainties in the full RT computations must be considered, as this was used to parameterize several ZTT terms. Polarization is not considered in the ZTT model or in Hydrolight, which is likely important [29,33]. There is a need in the research community for commercially available RT code with a focus on the ocean that includes full polarization, along with concomitant measurements of polarized radiance in the field. Neither the Hydrolight simulations or ZTT model included inelastic effects from DOM fluorescence, which could lead to bias errors as high as about 10% in the mid-visible [33]. Another potential source of uncertainty is the water Raman estimation [60], where bias errors as high as 10% can be typical.

Suitability for Inversion
The ZTT model can be used to directly solve for b b /a using least-squares minimization of a r rs measurement at any wavelength. This is shown in Figure 7 for the validation data set from Tonizzo et al. [29], assuming a constant b bp of 0.006. The typical value for V of 15 km was assumed. MAPE %δ abs in the inversion was 17%, which matches results from full RT simulations of closure in the forward direction. Residual errors thus appear almost entirely comprised of the inherent errors of the data set arising from IOP and radiometric measurements [29]. With respect to the use of P bb,ST (ψ), another potential approach to estimating VSF shape as it relates to the upcoming PACE mission is using planned multi-and hyper-angular polarimetry data from the Hyper-Angular Rainbow Polarimeter-2 (HARP2; contributed by University of Maryland, Baltimore County) and the Spectrometer for Planetary Exploration-one (SPEXone; contributed by the Netherlands Institute for Space Research) [94]. If information on the backward shape of the VSF could be gleaned with sufficient accuracy from angular polarimetry data (only a few angles would be necessary), this could be used directly in the ZTT model.
The ZTT model described here only addresses in-water RT, although the current approach to normalize water-leaving radiances measured by a satellite imager is expected to directly apply (see [15], Sections 3.1 and 6]). The BRDF is built-in to the Zaneveld model so an additional step to conceptually shift the geometry to a normalized geometry standard would not be required for IOP inversion. A normalized geometry (to nadir viewing or other geometry) could still be applied to intercompare reflectances in images. Translation to water-leaving radiance by accounting for all the reflection and refraction effects through the air-sea interface i.e., R (θ s ,θ v ) [6,15], can still be applied. Similarly, the effects of varying Earth-Sun distance and atmospheric attenuation can be removed following the current approach. Moreover, we would expect a practical implementation of this analytical model to have similar computational requirements as the current M02 LUT-based approach.    (Table 1) and the potential to invert r rs for b bp is also being investigated. The parameter b bp has not historically been associated with r rs and could provide new insights into particle composition in natural waters (e.g., [41]). Even a retrieval of b bp with relatively large errors could be useful for many research applications.
The ZTT model is fully spectral and can be applied over the full anticipated ocean color hyperspectral wavelength range for the future PACE mission (350-800 nm). A primary goal in model development was parameterization entirely in terms of a and b b . Inversion approaches, such as the General IOP (GIOP) inversion using a least-squares minimization to solve for IOP subcomponents [16] may consequently be readily applied. Spectral b bp and a pg can be solved using assumptions for subcomponent IOP spectra through least squares minimization of Equation (18). This is the step where empiricism prominently enters the problem as spectral shapes for IOP subcomponents must be assumed. This is being examined in ongoing work. Further investigation is needed to determine IOP subcomponent spectra with a high spectral resolution for error minimization over this full spectral range in inversion, which is currently being assessed by the PACE Science Team [13,95]. Hyperspectral (i.e., 5 nm) resolution for PACE is expected to improve fidelity in inverting for subcomponents (and potentially b bp ) as there will be higher degrees of freedom than for multi-spectral imagers [96,97].
Code for the ZTT model in MATLAB is available at ioccg.org.