## 2 Forward Model Development
### Average Cosine of the Asymptotic Light Field \(\overline{\mu}_{\infty}\)  
The \(\overline{\mu}_{\infty}\) term in Equation (7) must be expressed in terms of IOPs to invert. Using a fit to theoretical calculations of radiance fields by [PERSON] and [PERSON] [42], [PERSON] [28] recommended \(1/\overline{\mu}_{\infty}\) be modeled empirically with respect to the single scattering albedo \(\omega\) (=\(b/c\)) using a quadratic fit. Additionally, through more detailed RT computations, [PERSON] et al. [43] found a 4 th order dependency of \(\overline{\mu}_{\infty}\) on the albedo.  
In the study by [PERSON] and [PERSON] [40], \(\overline{\mu}_{\infty}\) was parameterized in terms of \(b_{b}/a\) instead of \(\omega\), 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 \(\overline{\mu}_{\infty}\) 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_{\text{rs}}\) without additional information, whereas \(b_{b}\) and \(a\) are (e.g., [PERSON] 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 \(\eta_{bb}\), the fraction of \(b_{b}\) attributable to molecular scattering, \(\eta_{bb}=b_{bw}/(b_{bp}+b_{bw})\)[6]. The natural range for \(\eta_{bb}\) is from \(\
eg 0\) to \(\
eg 0.98\)[40; 51], and the range used for \(b_{b}/a\) was \(10^{-4}\) to \(10^{-1}\). Since the water components of \(\eta_{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 \(\eta_{bb}\), the resulting fit was obtained for \(\overline{\mu}_{\infty}\):  
\[\begin{array}{rl}\overline{\mu}_{\infty}\left(\frac{b_{b}}{a},\eta_{bb} \right)\approx&\left[m_{1}(\log\eta_{bb})^{3}+m_{2}(\log\eta_{bb})^{2}+m_{3} \log\eta_{bb}+m_{4}\right]\left(\log\frac{b_{b}}{a}\right)^{3}\\ &+\left[m_{5}(\log\eta_{bb})^{3}+m_{6}(\log\eta_{bb})^{2}+m_{7}\log\eta_{bb}+ m_{8}\right]\left(\log\frac{b_{b}}{a}\right)^{2}\\ &+\left[m_{9}(\log\eta_{bb})^{3}+m_{10}(\log\eta_{bb})^{2}+m_{11}\log\eta_{ bb}+m_{12}\right]\log\frac{b_{b}}{a}\\ &\qquad+m_{13}(\log\eta_{bb})^{3}+m_{14}(\log\eta_{bb})^{2}+m_{15}\log\eta_{ bb}+m_{16}.\end{array} \tag{8}\]  
Coefficients \(m\) are provided in Appendix A, Table A2. Fits to simulated data were again made with the polyfit function from MATLAB. Absolute errors %\(\delta_{abs}\) for this fit vary from 0.19 to 3.5 for \(\eta_{bb}\) ranging from 0.98 to 0.0098, respectively. Since both \(b_{b}/a\) and \(\eta_{bb}\) are spectrally dependent, \(\overline{\mu}_{\infty}\) will be as well (not shown for clarity).