Inherent Optical Properties-Reflectance Relationships Revisited

Remote sensing reflectance, Rrs(λ) (sr), has been listed by the Global Climate Observation System as an essential climate variable thanks to its fundamental role in enabling observation of biogeochemical processes in the upper ocean (GCOS, 2003, and following updates; Groom et al., 2019). Theoretical studies have shown that apparent optical properties (AOPs) of the water such as sub-surface and above-surface remote sensing reflectance, rrs(0, λ) and Rrs(0, λ) respectively, are related to inherent optical properties (IOPs) such as absorption, a(λ) (m), and backscattering, bb(λ) (m), coefficients which in turn are controlled by the types and concentrations of in-water optical constituents. Significant effort has been devoted to deriving both IOPs and constituent concentrations from ocean color signals (IOCCG, 2006, 2014, 2019). Abstract Understanding the relationship between remote sensing reflectance, Rrs(λ) and the inherent optical properties (IOPs) of natural waters is potentially a key to improving our ability to determine biogeochemical constituents from radiometric measurements. These relationships are usually described as a function of absorption, a(λ), and backscattering, bb(λ), coefficients, with the literature providing various forms of equation operating on either bb(λ)/a(λ) or bb(λ)/[a(λ)+bb(λ)] to represent the impact of variations in light field geometries and changes in sea-water composition. The performance of several IOP-reflectance relationships is assessed using HydroLight radiative transfer simulations covering a broad range of Case 1 and Case 2 water conditions. While early versions of IOP-reflectance relationships assigned variability to associated proportionality factors (e.g., f/Q) or low-order polynomial functions, recent studies have demonstrated relationships between Rrs(λ) and bb(λ)/[a(λ)+bb(λ)] are wellcharacterized by nonlinear (high-order polynomial), monotonic functions. This study demonstrates that this approach is also valid for relationships operating on bb(λ)/a(λ) and that there is no intrinsic benefit to functions operating on bb(λ)/[a(λ)+bb(λ)] compared to bb(λ)/a(λ) for Case 2 waters, contrary to recent suggestions in the literature. In all cases it is necessary to carefully consider the performance of best fit relationships across the full range of variability of IOPs and Rrs(λ), with higher order polynomials required to enable equivalent performance across the range of natural variability. The analysis further demonstrates insignificant wavelength sensitivity across the visible region, limited sensitivity to changes in solar zenith angle and extends to relationships for below surface remote sensing reflectance, rrs(λ).

• The mathematical form and performance of relationships between above and below surface remote sensing reflectance and inherent optical properties (IOPs) of natural waters are examined for both forward and inverse applications • Two forms of IOPs relationship that have previously been presented in the literature are shown to be equivalent in performance across a wide range of optically complex water conditions • Performance of new and previously published relationships varies across the range of natural optical variability and higher order polynomial functions are needed for accuracy in some cases The quasi-single scattering approximation (QSSA) (Gordon, 1973;Gordon et al., 1975) models are widely employed in standard ocean color processing and applications because they provide an explicit relationship between IOPs and r rs (λ), as formulated by Gordon et al. (1988), coupled with a simple relationship for converting R rs (λ) to r rs (λ) (Lee et al., 2002). The weakness of QSSA models is their inability to account for multiple scattering effects and their lower accuracy, compared to radiative transfer (RT) codes (Werdell et al., 2018). On the other hand, forward RT models such as HydroLight (Sequoia Scientific, Inc.) can provide the full radiance distribution below and above the water, as a direct solution to the RT equation (Mobley, 1994). Therefore, given a set of input IOPs and appropriate boundary conditions, it is possible to generate precise estimates of a full set of AOPs, including r rs (0 − , λ) and R rs (0 + , λ). However, the inverse problem of determining IOPs from AOPs is not straightforward and relies on empirical and/or semi-analytical relationships.
It is intuitively obvious that reflectance signals will increase with backscattering and decrease with absorption. Early work in this area focused on irradiance reflectance below the sea surface, R (0 − ,λ) (nondim), defined as the ratio of the upwards to downwards planar irradiances, E u (0 − ,λ) (Wm −2 nm −1 )/E d (0 − ,λ) (Wm −2 nm −1 ), as this was both practically measureable with available instrumentation and computationally convenient (e.g., Kirk, 1984). The results of these early studies (Gordon et al., 1975(Gordon et al., , 1988Gordon & Morel, 1983;Morel & Prieur, 1977) confirmed that R (0 − ,λ) is a function of the absorption and backscattering coefficients. However, in the literature, at least two different expressions are found to approximate this relationship. Morel and Prieur (1977), modeling the results from RT calculations, found that while Gordon et al. (1975) showed that where F represents a polynomial function with factors for up to third order given in the original paper. The factor f in Equation 1 has been found to be variable in natural conditions (Morel & Gentili 1991, 1996 accounting for most of the directional effects due to changes in the light field, degrees of multiple scattering and in water bio-optical characteristics. The exact form of F in Equation 2 (second or third order polynomial) and the associated polynomial coefficients have not been unambiguously established for all situations, though in several cases a second order polynomial has been adopted (e.g., Feng et al., 2005).
With the advent of ocean color satellites measuring top of atmosphere radiance, L TOA (λ) (Wm −2 sr −1 nm −1 ), the focus of effort shifted toward remote sensing reflectance, R rs (0 + , λ), defined as the ratio of the water-leaving radiance, L w (0 + , λ), to downwards irradiance, E d (0 + , λ), just above the air-sea interface (Mobley, 1994). R rs (0 + ,λ) has commonly been used in satellite oceanography as input to inversion methods for deriving IOPs and water constituent concentrations (IOCCG, 2006) as it minimizes the effects of environmental conditions such as the influence of sun angle. R rs (0 + , λ) measurements are converted to below surface r rs (0 − , λ) using Equation 3.
where T is a transmission factor incorporating information including the Fresnel transmittance from water to air, the refractive index of seawater (invoking the n 2 law of radiances) and several additional factors to deal with propagation of downwards irradiance (see Mobley, 1999 and IOCCG, 2019 for more details). The resulting conversion factor T exhibits only limited variability (0.50 < T < 0.57) and is usually assumed to have a value of 0.54. Morel and Prieur (1977) adapted Equation 1 to express subsurface remote sensing reflectance r rs (0 − ,λ) as a function of the ratio b b (λ)/a(λ) (w M from here onward): where Q(sr) is the bidirectional function defined by Morel and Gentili (1996) as the ratio of upwards irradiance E u (0 − ,λ) to upwards radiance, L u (0 − ,λ), and expresses the nonisotropic character of the radiance distribution.
In contemporary ocean color processing (Werdell et al., 2018), r rs (0 − ,λ) is more often expressed as a function of IOPs following the approach of Gordon et al. (1988) Taking Equations 3-5 together, it is clear that R rs (0 + ,λ) can be related to IOPs through either w G or w M . However, the apparently simple form of Equations 4 and 5 is deceptive. From the outset it was shown that f and Q were variable. In a series of papers, Morel and Gentili (1991, 1993, 1996 explored variability in f and Q, with the f/Q factor found to be variable in the range 0.075-0.12 in oceanic (Case 1) waters and affected by solar zenith angle, sensor viewing geometry, in-water constituent concentrations and wavelength. Furthermore, in coastal (Case 2) waters, additional concentrations of CDOM and mineral particles that do not co-vary with the chlorophyll concentration were expected to influence the variability of f/Q, though this variability was found to be minimal in the nadir viewing direction and the f/Q factor was almost insensitive to different wavelengths when the Sun is at the zenith (Loisel & Morel, 2001). The net result of this work was establishment of look up tables (LUTs) for Case 1 waters (Morel & Gentili, 1996; to provide estimates of f/Q which, while certainly addressing the complexity of the underpinning RT physics, leaves Equation 4 difficult to use in practice. In addition, attempts have been made to provide suitable LUTs for Case 2 waters (e.g., Hlaing et al., 2012;Park & Ruddick, 2005;Salem et al., 2017), but none of these approaches have been adopted unanimously.
One possible reason for the relative popularity of IOP-reflectance relationships operating on w G is that the quadratic form expressed in Equation 5 captures some of the nonlinear behavior in the relationship with IOPs that is less obviously elucidated in the LUT approach for f/Q. The limiting factor for Gordon style versions of the relationship with IOPs stems from failure to achieve consensus on a single form or set of coefficients that performs equally well across the known range of variability for natural waters. Gordon et al. (1988) suggested a quadratic form with g 0 = 0.0949 and g 1 = 0.0794 for Case 1 waters, while for highly scattering coastal waters Lee et al. (1999) suggested g 0 = 0.084 and g 1 = 0.17. Later Lee et al. (2002), aiming at applying the forward model to both coastal and open-ocean waters, proposed average values of the coefficients from previous studies, g 0 = 0.0895 and g 1 = 0.1247. Subsequently, focusing on above surface R rs , Albert and Mobley (2003) and Park and Ruddick (2005) developed fourth order polynomial relationships in deep and shallow waters, while more recently Lee et al. (2011) and Hlaing et al. (2012) presented second and third order polynomial variants respectively. In all cases, the effect of these efforts was to directly address apparent variability in the relationship between reflectance and IOPs using a relatively simple equation that was easier for end-users to implement than working with LUTs. Of course it is also true that each version of the relationship is potentially restricted by the bounds of the modeling parameters used to establish that particular relationship, so it is not surprising that a single optimum relationship has not been forthcoming.
As time has passed, the community has gradually become accustomed to working with IOP-reflectance models operating on w G . In recent years, there have been claims that w G is generally applicable to all waters, including optically complex coastal waters (D'Sa & Miller, 2005) while w M has been suggested to be inappropriate for waters where b b (λ) cannot be considered negligible compared to a(λ) (e.g., Loisel & Morel, 2001;Morel & Muller, 2002). This has been interpreted by some as suggesting that the w M approach is not suitable for use in coastal waters (Marrable, 2014). Anecdotally, we are aware of concerns being raised about the end case scenario of zero absorption causing w M to go to infinity. However, it is difficult to find strong evidence in the literature to support the idea that w M models are less effective than w G for coastal waters or elsewhere, and the zero absorption argument is a red herring as the absorption simply doesn't go to zero in any natural water. On the other hand, while w M is a priori unbounded, w G lies within the interval [0, 1]. However, unless dealing with extremely turbid waters (total suspended solid >> 100 gm −3 ) the maximum theoretical value of 1 for w G is barely approached. In this study, for example, the maximum value of w G was 0.6, having considered concentrations of mineral suspended particles up to 100 gm −3 , which is consistent with values found for moderate to very turbid waters.
While there has been considerable improvement in our understanding of IOP-reflectance relationships since the pioneering studies by Morel and Gordon, there remains confusion in the field about the relative merits of w M and w G approaches. Moreover, there is an outstanding requirement to develop easily implemented relationships that permit end users to relate IOPs and reflectance values for both above and below surface reflectances corresponding to remote sensing and in situ applications. Most importantly, it is essential that the performance of any such relationship is equivalent across the wide range of concentrations of optical constituents found in natural waters and can be applied equally well in both clear Case 1 waters and optically complex Case 2 waters.
The aim of this study is to exploit the strength of forward RT modeling to generate a consistent set of IOPs and corresponding simulated r rs and R rs values which enables investigation of IOP-reflectance relationships across a wide variety of optically complex water conditions. The majority of the work considers the nadir-viewing direction with the Sun at zenith, reflecting fully normalized conditions for remote sensing applications, but also briefly considers the impact of other solar zenith angles to better represent in situ measuring conditions. Both the sub-surface remote sensing reflectance, r rs , and the remote sensing reflectance in air, R rs are considered, again to better represent different practical scenarios. This study broadly follows the approach of Wong et al. (2019), but the analysis is extended to include w M , and both above (R rs ) and below (r rs ) surface reflectances. In all cases, additional care is taken to test the performance of IOP-reflectance relationships over each decade of variability of w G and w M .

Field Data
All in situ data used in this study were collected in March 2009 during the BP09 cruise in the Ligurian Sea ( Figure 1). Located in the northwest part of the Mediterranean Sea, this area includes both deep clear oceanic waters (considered Case 1) and shallow turbid coastal waters (considered Case 2). Among 60 sampled stations, 11 offshore stations and 23 onshore stations returned sufficient data set for the purpose of this work. A detailed description of the sampling location, data and methods can be found in McKee et al.

In Situ Optical Measurements
In situ nonwater absorption a n (λ) and attenuation c n (λ) (m −1 ) measurements were made with a WETLabs 25 cm path length AC-9 instrument, operating at nine wavelengths across the visible-NIR (412,440,488,510,532,555,650,676, and 715 nm, 10 nm FWHM) and corrected for salinity and temperature using data from a Seabird SBE 19 plus CTD. The proportional method by Zaneveld et al. (1994) was applied to correct AC-9 absorption data. In situ backscattering b b (λ) measurements were made with a WETLabs BB-9 backscattering meter. BB-9 data were corrected for path length absorption using corrected AC-9 absorption data and the backscattering of pure water (Smith & Baker, 1981) was subtracted to obtain the particulate backscattering coefficient, b bp (λ), as detailed in Lefering et al. (2016). BB-9 data were linearly interpolated where necessary to match AC-9 wavelengths.

Laboratory Measurements
The absorption of all dissolved and suspended components minus water was measured using a point source integrating cavity absorption meter (PSICAM), with an estimated accuracy of ±2% even in highly LO PREJATO ET AL. 4 of 20 10.1029/2020JC016661 turbid waters (Röttgers et al., 2005;Röttgers & Doerffer, 2007). Measurements were made using fresh Milli-Q ultrapure water references and then corrected for salinity and temperature effects on water absorption (Röttgers & Doerffer, 2007;Röttgers et al., 2014). Measurements of absorption by colored dissolved organic matter, a cdom (λ), were performed using a 1 m liquid waveguide capillary cell (LWCC) with an Ocean Optics USB2000 mini spectrometer, with a noise range of ±0.0001 m −1 at 532 nm (Lefering et al., 2017). In this study, the concentration of colored dissolved organic material is represented by its absorption at 440 nm, CDOM = a cdom (440). Total particulate absorption, a p (λ) was derived by subtracting a cdom (λ) from PSICAM nonwater absorption. Total particulate absorption was also obtained by measuring the optical density of particles, OD p (λ) (nondim), retained on a 25 mm GF/F filter using a Shimadzu UV-2501 PC benchtop spectrophotometer. After applying dilute bleach and rinsing with filtered seawater to remove algal pigments, the measurement was repeated for each sample to provide the detrital optical density of nonalgal particles, OD det (λ). The particulate and detrital absorption coefficients, a p (λ) and a det (λ) respectively, were obtained from: where a x (λ) is the required absorption coefficient, A fp is the exposed area of the filter pad, V f is the volume of sample filtered (m 3 ), β is the path length amplification factor (nondim) and OD x (λ) is the required optical density. A linear regression approach (Lefering et al., 2016;McKee et al., 2014) using a combination of PSICAM and LWCC data was employed for the determination of path length amplification factors (β) and scattering offset corrections for both filter pad absorption measurements. Absorption by phytoplankton, a ph (λ), was derived by subtracting absorption of detrital particles from the total particulate absorption, a ph (λ) = a p (λ)−a det (λ).
Concentrations of chlorophyll-a (CHL) and total suspended solids (TSS) were measured in triplicate by colleagues from Management unit of the North Sea mathematical models (MUMM) with averaged results reported here. CHL concentrations were determined by standard HPLC measurements: triplicate samples were analyzed using a reversed phase, acetone-based method with a C18 column and a Jasco FP-1520 fluorescence detector. TSS concentrations were estimated using preashed, rinsed and preweighed 47 mm GF/F filters, careful rinsing each sample with ultrapure water. Filters were stored frozen and returned to the lab where they were dried and reweighed. A simple relationship between TSS and CHL for Case 1 waters, obtained by least squares regression forced through zero, gives biogenic TSS (TSS bio ) per unit CHL (TSS bio = 0.234 CHL). This is used to estimate the concentration of biogenic TSS for Case 2 waters, with mineral (i.e., nonbiogenic) suspended solids (MSS) obtained for Case 2 waters by subtracting the estimate of TSS bio from total TSS .

Bio-Optical Model Development
A linear bio-optical model, based on spectral material-specific IOPs (SIOPs) and populated with a wide range of log-spaced constituent concentrations (CHL, MSS, and CDOM) was employed for describing a broad range of optical water conditions. In this study, the SIOPs have been determined by simple linear regression of partitioned IOPs against related optically dominant constituent concentrations, following McKee and Cunningham (2006). This methodology minimizes the effect of the propagation of measurement uncertainties in both optical and biogeochemical measurements .
Following the methodology reported in Ramírez-Pérez et al. (2018) for Case 2 waters, each total IOP can be considered as the sum of partial IOPs, and each partial IOP can be expressed as the product of SIOPs, and associated constituent concentrations: where the subscripts represent the following five bio-optical constituents: phytoplankton (ph), biogenic detritus (bdet), nonbiogenic detritus (ndet), colored dissolved organic material (cdom) and pure water (w). The constituent concentrations are: chlorophyll, CHL, absorption of colored dissolved organic material at 440 nm, CDOM, and mineral suspended solids, MSS (the nonbiogenic detrital component of TSS). Here the detrital particulate absorption has been considered as the sum of two separate biogenic and nonbiogenic components. Unfortunately, it is not possible to experimentally partition scattering and backscattering measurements so the level of discrimination possible for these parameters is reduced. Biogenic partial IOPs are assumed to co-vary with CHL, while the nonbiogenic partial IOPs are assumed to co-vary with MSS.
A set of material-specific absorption and scattering coefficients was previously prepared for Ramírez-Pérez et al. (2018). However, in the present study, only values at each standard AC-9 wavelength have been considered (Figures 2a-2f, black lines). These SIOPs have been determined by simple linear regressions forced through zero of partitioned IOPs against associated constituent concentrations from surface water samples. The optimal SIOP value has been estimated using the linear regression slope with the 95% confidence interval (CI) representing the associated uncertainty ( Figure 2, red dashed lines). In addition to the previously provided SIOPs for partitioned absorption and scattering coefficients, the material-specific backscattering coefficients are presented here following the same methodology. The chlorophyll-specific phytoplankton backscattering coefficient, , shown in Figure 2g, was determined using data from offshore stations (considered Case 1 waters), where the particle population is assumed to be biogenic in origin. However, it is assumed to be representative of algal and biogenic detrital backscattering for all stations, offshore and onshore. For onshore stations (considered Case 2 waters) it is assumed that there will be an additional nonbiogenic contribution, , that can be found after subtraction of the water and biogenic contribution ( . The MSS-specific nonbiogenic detrital backscattering coefficient, shown in Figure 2h was obtained by regressing

Radiative Transfer Simulations
Input IOPs for the simulations were generated by populating the bio-optical model, described in the previous section at 9 standard AC-9 wavelengths (412, 440, 488, 510, 532, 555, 650, 676, and 715 nm). A total of 1,690 unique combinations of constituent concentrations and associated IOPs were calculated from logspaced distributions of constituent concentrations in specific ranges (0.01< CHL <100 mg/m 3 , n = 13; 0.01< MSS <100 g/m 3 , n = 13; 0.01< CDOM <10 m −1 , n = 10). The data set of modeled IOPs was arranged in the form of 1690 virtual AC-9 and BB-9 type files in the MATLAB® environment (MathWorks Inc.). Nonwater absorption and beam attenuation coefficients were given as functions of depth and wavelength via the user-supplied AC-9 type data files. The RT numerical model HydroLight 5.2 (Sequoia Scientific Inc.) was used to process the input IOP data and generate a synthetic data set of remote sensing reflectance spectra, both below and above the sea-air interface, r rs (0 − , λ) and R rs (0 + , λ), respectively. As all combinations of the log-spaced biogeochemical concentrations were included in the RT runs, the resulting input data set broadly represents the conditions of a substantial portion of natural coastal waters. Based on the considered constituent concentration ranges, the total absorption at 412 nm, a (412), was in the range 0.0218-27.1346 m −1 , while the total backscattering coefficient, b b (412), varied from 0.0036 to 2.3240 m −1 , where the minimum (maximum) values represent the total IOPs when all the constituent concentrations are at their lowest (highest). Pure water absorption and scattering values were taken from Pope and Fry (1997) and Smith and Baker (1981). All the input files were processed in HydroLight using the ABACBB model which provides sample dependent phase functions, dynamically generated through the Mobley et al. (2002) parameterization of the Fournier-Forand scattering phase function. This method allows estimation of wavelength-dependent scattering phase functions based on the particle backscattering ratio, b bp (λ)/b p (λ) which showed significant variability at each waveband (e.g., between 0.01 and 0.03 at 412 nm and between 0.007 and 0.025 at 650 nm). The subsurface remote sensing reflectance, r rs (0 − ,λ) calculated by HydroLight as the ratio of upwards radiance L u (0 − ,λ) to downwards irradiance E d (0 − ,λ) and the remote sensing reflectance in air, R rs (0 + ,λ), calculated as the ratio of water-leaving radiance L w (0 + ,λ) to downwards irradiance LO PREJATO ET AL. 8 of 20 10.1029/2020JC016661 . These runs were repeated for solar zenith angles of 30 and 60° in order to investigate how the solar geometry affects the results. For the sake of simplicity, the water body was modeled as a homogeneous and infinitely deep medium. Sky radiance was modeled using RADTRAN-X (Gregg & Carder, 1990) with no wind and the mean Earth-Sun distance was employed. Other atmospheric conditions such as sea-level pressure, relative humidity, horizontal visibility, and ozone concentration were set at default values, which are described in the HydroLight 5 technical documentation (Mobley & Sundman, 2008). Raman (inelastic) scattering by water itself is ubiquitous and was considered since it is easy to model without any extra information, while fluorescence emissions due to chlorophyll and CDOM were not included, due to lack of information about fluorescence quantum yields.

Establishing the Nature of IOP-Reflectance Relationships
Variations in the bidirectional properties of the radiant flux leaving the water and returned toward the atmosphere have been extensively studied (Morel & Gentili, 1991, 1993, 1996. The term f/Q in Equation 4 implicitly describes the magnitude of these variations, mainly due to changes in the light field geometry and bio-optical characteristics of the water body. Figure 3 shows distributions of f/Q as obtained in this study from RT simulations, calculated from Equation 4. These results are broadly consistent with Morel and Gentili (1993)  noted to be affected by variations in Sun angles and viewing geometries. This is confirmed in Figures 4c and  4d where it can be seen that changing solar zenith angle has a small but defined impact on relationships between R rs and both w M and w G . In the present study, unless otherwise stated, the results refer only to a zenith Sun and a vertical viewing direction, meaning that the variability of the f/Q factor here is due mainly to differences in water constituent concentrations.
The RT simulations can also be used to assess variability in the transmission factor T by rearranging Equation 3. Mobley (1999) suggested a range between 0.50 and 0.57 for typical ocean waters and solar and sensor geometries. In this study, T was found to vary in the range 0.52-0.63 considering the Sun at the zenith and the sensor in the vertical direction, with the differences in ranges between the two studies presumably associated with the range of water conditions (and therefore IOPs) considered. Overall, it is reasonable to suggest that the variability described by f/Q and T factors observed in this study are broadly consistent with previous studies. This preliminary check is encouraging because it confirms the capability of the synthetic data set to realistically describe the variability of both f/Q and T, allowing a systematic study of the IOP reflectance relationships in optically complex waters. No angular dependency has been considered in this part of the study, therefore the variability of f/Q and T mentioned here are determined by the bio-optical properties of the water body only. We would anticipate further variability in both f/Q and T in the case of extending the study to other sun-sensor angle combinations.
Empirical relationships between r rs (λ) and IOPs (primarily operating on w G ) have already been expressed in several studies through quadratic or higher polynomial functions (Albert & Mobley, 2003;Gordon et al., 1988;Lee et al., 2002Lee et al., , 2004Wong et al., 2019) with several others focusing on relationships between above surface R rs (λ) and IOPs ( 2005). This study investigates relationships between IOPs and both r rs (λ) and R rs (λ) considering both w M and w G forms, and also considers both forward and inverse directions.
Plotting the HydroLight output r rs (λ) values against the input IOP terms, w M and w G , (Figures 4a and 4b), the relationships between r rs (λ) and either w M or w G are found to be nonlinear but, crucially, monotonic. It is evident that a significant amount of the apparent variability in f/Q is associated with the nonlinear nature of the relationship with w M . Whereas the frequency distributions shown in Figure 3 do not offer any means of predicting a particular value of f/Q, Figure 4 suggests the possibility of establishing relationships between r rs and either w M or w G with strong predictive power in both the forward and reverse directions. A first visual comparison of the two plots in Figure 4 suggests that both approaches (w M and w G ) are equivalently valid and confirms that the dependence of the subsurface remote sensing reflectance on wavelength (400 nm<λ < 700 nm) is minimal in Case 2 waters, when considering the nadir-viewing direction only and the Sun at the zenith (Loisel & Morel, 2001). Changing the solar zenith angle (θs = 0°; 30°; 60°) generates similar but nonidentical nonlinear, monotonic relationships (Figures 4c and 4d). The main focus of this study is to establish IOP-reflectance relationships for the fully normalized case of both sun and sensor at zenith, with LO PREJATO ET AL. 11 of 20 10.1029/2020JC016661 Figure 6. Performances of second to seventh order polynomials (colored lines) in fitting HydroLight simulated data (gray points). In the upper two plots the ratios of (a) w M and (b) w G are considered as independent variables while in the lower two plots (c and d) r rs is considered the independent variable. All plots are on log-log scales to better show the behavior of polynomials across the range of values. Note that some of the best fits produced negative values that cannot be represented using this scaling (see truncated lines in the upper two plots).
application to fully processed ocean color data in mind. We note that other combinations of sun-sensor angles require further additional modeling, consistent with Hlaing et al. (2012).
Plotting the HydroLight output R rs (λ) values against the input IOP terms, w M or w G , (Figures 5a and 5b), the relationships between R rs (λ) and either w M or w G have an analogous behavior compared to the corresponding cases of subsurface reflectance, with slightly less sensitivity to changing solar angles (Figures 5c and 5d).

Determination of Optimal IOP-R rs Relationships
This study considers a wide input data range and applies a systematic curve fitting procedure in order to find second-to-seventh order polynomial models for the simulated set of data. The polynomial curve fitting has been performed in the Matlab environment using the polyfit operator. This built-in function returns the coefficients of a polynomial of the chosen degree, which is a best fit for the dependent variable in a least squares sense, along with the 95% confidence bounds for each coefficient. Different polynomial forms have been investigated for both forward and inverse applications, considering r rs (λ) versus IOPs (Figures 6a  and 6b) and vice versa (Figures 6c and 6d), focusing on the case of a solar zenith angle,   θ 0 S , and nadir viewing sensor direction,    0 . Both w M and w G have been used as representations for IOPs and the performance of each polynomial has been assessed across the full range of data (Figure 7). Because high-order polynomials cannot be easily inverted algebraically (González-Cardel & Díaz-Uribe, 2006) and numerical inversions are computationally demanding, a practical choice was to apply the same fitting procedure for estimating IOPs from r rs (λ) instead of mathematically inverting the forward model. It is clear from Figure  Comparison between second to seventh orders of polynomial models, in terms of mean absolute error: (a and b). The inherent optical properties terms are considered as independent variables in the fitting process; (c and d) r rs is considered as an independent variable. Data have been binned by order of magnitude therefore mean absolute error is referred to the order of magnitude of the independent variable.
6 that there is significant divergence in the performance of each best fit polynomial for low values of both w M and w G , and similarly when these parameters are derived from r rs (λ). In contrast, performance for high values of r rs (λ), w M and w G is less sensitive to choice of polynomial. It is therefore important to consider performance over all relevant ranges of signal strength when attempting to establish a set of optimal empirical relationships.
The performance of each best fit polynomial is primarily assessed using the mean absolute error (MAE), a commonly used forecasting metric that describes the typical magnitude of the residuals, giving the average error in units of the variable of interest (Figures 7a-7d). When estimating r rs (λ) from w M or w G (Figures 7a  and 7b) it is possible to directly compare the two approaches because the derived term in both cases is r rs (λ). However, when comparing retrievals of b b (λ)/a(λ) and b b (λ)/[a(λ)+b b (λ)], it is important to bear in mind that the two variables are not numerically equivalent (Figures 7c and 7d). As both IOPs and r rs (λ) values vary over several orders of magnitude, the data set has been split into subranges that permit evaluation of performance for each order of magnitude. The results of this systematic exercise show that higher order polynomial models provide better performance than second order models to ensure similar goodness of fit across all decades of signal variation. Based on analysis of Figure 7, the recommended relationships in retrieving r rs (λ) from w M and vice versa are sixth and fifth order polynomials, respectively (Figures 7a and  7c). A third order polynomial is the suggested choice for retrieving r rs (λ) from w G , (Figure 7b), while a sixth order polynomial is needed for estimating w G from r rs (λ) (Figure 7d). The coefficients for these four polynomial relationships, with related 95%CI, are reported in Note. The polynomial relationships are in the form y = p 0 + p 1 x + p 2 x 2 + … + p n x n , where n is the order of the polynomial. Each polynomial coefficient is supplied with associated 95% confidence bounds. Abbreviation: IOP, inherent optical properties.

Table 1 Coefficients of Best Fit Polynomials for IOP-R rs (λ) Relationships and Vice Versa, With Sun at zenith
example; in terms of computational speed, but overall there is no significant accuracy benefit of one form of IOP expression over the other. The results of this exercise show that with suitably careful selection of empirical relationship, using w M or w G is broadly equivalent. This is particularly true for high reflectance values typical of Case 2 waters, strongly contradicting previous suggestions that the w M form was inappropriate for these waters (e.g., Marrable, 2014).

Determination of Optimal IOP-R rs Relationships
The same systematic curve fitting procedure applied in Section 3.2 for r rs (λ) has been employed in order to find second-to-seventh order polynomial models for the simulated set of data considering R rs (λ) versus IOPs (Figures 8a and 8b) and vice versa (Figures 8c and 8d). The assessment of the fitting performance of each polynomial is shown in Figure 9, with MAE supplied. Similarly to the previous case, higher order polynomials are suggested to ensure similar goodness of fit across all decades of signal variation, with the greatest impact occurring at small data values. From a closer analysis of Figure 9, the recommended relationships in retrieving R rs (λ) from w M and vice versa are the fifth order polynomials in both cases (Figures 9a and 9c, respectively). On the other hand, for retrieving R rs (λ) from w G , a fourth order polynomial is the suggested LO PREJATO ET AL. 14 of 20 10.1029/2020JC016661 choice (Figure 9b), while a seventh order polynomial is needed for estimating w G from R rs (λ) (Figure 9d). The coefficients for these four polynomial relationships, with related 95%CI, are reported in Table 2.
Using the suggested models, R rs (λ) can be expressed as a function of the absorption and backscattering coefficients of the water and it is equally valid to represent the IOPs using w M or w G . The limiting factor is not choice of either w M or w G , rather it is in selecting an appropriate formulation to represent variation the nonlinear relationship with IOPs.

Comparison With Previous Studies
Many studies (Albert & Mobley, 2003;Chami et al., 2006;Lee et al., 2004;Mobley et al., 2002;Morel & Gentili, 1991, 1993Park & Ruddick, 2005;Wong et al., 2019) have shown that the nonlinear increase of reflectances with IOPs (in the form of either w M or w G ) can be attributed to multiple scattering. However, the existing established models for representing r rs (λ) (Gordon et al., 1988;Lee et al., 2002Lee et al., , 2004 or R rs (λ), based on low-order polynomial relationships, are unable to fully account for the effect of multiple scattering which is particularly relevant in Case 2 waters (Wong et al., 2019). Interestingly, the results presented above suggest that the performance of low-order polynomials tends to be worst for clear waters. It is clear from this analysis that higher order polynomials are generally necessary in order to achieve optimal performance across the full range of natural variability.
LO PREJATO ET AL. 15 of 20 10.1029/2020JC016661 Figure 9. Comparison between second to seventh orders of polynomial models, in terms of mean absolute error, as in Figure 7 but considering above surface R rs (λ).
Previous studies have typically reported best fit IOP-reflectance relationships for conditions with the sun placed away from zenith. For the purpose of comparison with these previous studies, we have assembled best fit polynomial parameters for the case of the Sun at θ s = 30° and the sensor in nadir viewing direction (Table 3). Figure 10a shows comparison with a number of earlier studies, while Figure 10b shows the MAE on r rs (λ) by order of magnitude of w G . The performance of our proposed Poly5 model is very close to the version of Gordon's model with coefficients suggested by Lee et al. (2002) which was intended to be applicable to both Case 1 and Case 2 waters. Wong et al. (2019) presented an updated version of a model by Lee et al. (2004) that operates on molecular and particulate components (u w and u p ) separately. Applying the Wong et al. (2019) approach to our data set presents smaller MAEs for lower values of w G , but higher MAE for higher values of w G . It is important to note that the family of curve models shown in Figure 10a falls within the 95% confidence bounds of the proposed Poly5 relationship at the high end of the range (Figure 10c) while the 95% confidence bounds at the low end ( Figure 10d) are an order of magnitude greater than the MAE differences between the various proposed relationships for the equivalent range in Figure 10b. This comparison ( Figure 10) does not identify a clear cut "winner" but does illustrate clearly the need to consider performance across the full range of values of interest.

Discussion and Conclusions
In this study, a series of RT simulations across a very broad range of optical conditions has been conducted using the HydroLight RT code. Input IOPs have been supplied through a linear spectral bio-optical model, varying the concentrations of in-water optical constituents across a wide range of concentrations to simulate natural variability. The resulting set of synthetic above and below surface reflectances have As in Table 1 but for IOP-R rs (λ) relationships. Abbreviation: IOP, inherent optical properties. enabled investigation of relationships with IOPs expressed as either w M or w G . Consistent with previous studies, the results demonstrate that the relationship between IOPs and either r rs or R rs are highly predictable and can be well-modeled by a nonlinear but monotonic curve, which is not significantly wavelength dependent for 400 nm<λ<700 nm. However, it has been shown here that in order to achieve consistent levels of prediction performance across the broad range of optical coastal water conditions considered, it is generally necessary to consider higher order polynomial relationships. Comparison of w M and w G variants demonstrate effectively equivalent performance: there is no immediate advantage or disadvantage to use of either form in coastal waters. This contradicts previous claims that w M was not valid for Case 2 waters (e.g., Loisel & Morel, 2001;Marrable, 2014;Morel & Mueller, 2002) and highlights the general proposition that either combination of absorption and backscattering contains sufficient relevant information to both predict and interpret remote sensing signals. There may be situations where future end users may have reason to prefer use of one form over another and this study suggests there is no reason to be concerned over either.
We note that the relationships presented here are empirical and only valid for the range of conditions from which they have been derived and due to their high order are likely to deviate rapidly from reality beyond these conditions. Nonetheless, the range covered is particularly broad so the results are expected to be generally useful for a wide range of natural waters. The details of the bio-optical model, such as the use of linear relationships between IOPs and constituents, used to set up the RT simulations are potentially open to criticism. However, it should be noted that it would have been possible (and in some senses possibly preferable) to have produced a comprehensive suite of values of b b (λ) and a(λ) with no reference to a bio-optical model at all. The purpose of the bio-optical model in this case is simply to assure readers that the range of IOPs considered genuinely covers a broad range of optical water conditions. We also note that The polynomial relationships are in the form y = p 0 + p 1 x + p 2 x 2 + … + p n x n , where n is the order of the polynomial. Each polynomial coefficient is supplied with associated 95% confidence bounds. Abbreviation: IOP, inherent optical properties. the results are restricted to the wavelength range (400-700 nm) and sun-sensor angle conditions used in the simulations. Relationships derived for both sun and sensor in the zenith position (Tables 1 and 2) are offered for consideration for ocean color remote sensing applications where data have been fully corrected for BRDF effects. The relationships provided in Table 3, for the sensor at zenith and Sun at 30° from zenith, are offered for applications involving vertically profiling in situ radiometers; the differences in practice are small but potentially significant. Figure 10. (a) Scatter plot of r rs (λ) on w G from HydroLight simulations (gray points) compared to previous models (Gordon et al., 1988;Lee et al., 1999Lee et al., , 2002Wong et al., 2019) and the model Poly5 from the present study.   θ 30 S and nadir viewing sensor direction have been considered for this comparison, with Poly5 considered as the best performing polynomial model among second-seventh-order polynomials at  30 solar angle. (b) Mean absolute error as in Figure  7b has been calculated for each model. The 95% confidence bounds of the Poly5 relationship are shown at both the high end (c) and the low end (d) of the data distribution, using a linear and log-log scale, respectively.