As mentioned in the introduction, the “Four-term” IOP model manages all components in Eqs. (1) and (2), separately. The IOPs of pure water are known and are not the subject of further modeling. The fundamental concepts of constructing IOPs of gelbstoff, detritus, and phytoplankton – so that it fits to our observational data – are illustrated in the following. The corresponding parameterizations, which we determined based on analyses of measured data, are given in section 3.2.

The absorption of gelbstoff in surface water, as a proxy of CDOM, can be expressed as:

where the reference wavelength, λ 0 , is usually at 440 nm, and a g + ( λ ) is the normalized a g ( λ ) at λ 0 , which has typically a nearly exponential shape in the visible spectral range, which is most relevant in the ocean color context. To substitute the exponential function into Eq. (5), a g ( λ ) can be written as

where S g is the spectral slope estimated by nonlinear regression for a specific wavelength range, such as 350 to 500 nm (Babin et al., 2003b). The form of Eq. (6) has been widely accepted to compare the CDOM absorption coefficients across different systems but may have a risk when extrapolate to longer wavelengths, e.g., > 500 nm.

The absorption coefficient of detritus, a d ( λ ) , is assumed to be controlled by both phytoplankton [Chl] and inorganic suspended matter, [ISM]. Following the concept of Ramírez-Pérez et al. (2018) and Lo Prejato et al. (2020), the absorption coefficient is constructed as:

where a b d * ( λ ) is the Chl-specific biogenic detritus absorption coefficient in the unit (m²/mg) and a m d * ( λ ) is the mass-specific minerogenic detritus in (m²/g). Note that the units of [Chl] and [ISM] are (mg/m³) and (g/m³), respectively.

The scattering of detritus is determined by the subtraction of its absorption from its attenuation:

The attenuation coefficient of detritus, c d ( λ ) , is described as a power law function to avoid any negative values:

where the reference wavelength, λ 0 , is often at 550 nm for its lower susceptibility to phytoplankton absorption, and γ d is the power law exponent. Using Eq. (9), we can determine the spectral shape of c d . To ensure a positive scattering value obtained from Eq. (8), we calculate the magnitude of its attenuation based on the single-scattering albedo at λ 0 , ω d ( λ 0 ) , which is a relative IOP parameter and can be empirically obtained. The magnitude of the spectrum, represented by c d ( λ 0 ) , can be then calculated by:

The absorption coefficient of phytoplankton is expressed as

where a p h * in (m²/mg) is the Chl-specific absorption coefficient, and [Chl] is commonly regarded as a proxy of the total phytoplankton concentration.

To capture the natural variability in phytoplankton IOPs due to the occurrence of different taxonomic groups, several phytoplankton groups with different optical properties and pigment compositions are considered. The Chl-specific absorption coefficient for phytoplankton is then expressed as:

where f p h , i is the fraction from 0 to 1 of each group in sum equal to 1, and N p h is the number of different groups used. f p h , i represents the proportionate contribution of each group to the total chlorophyll a concentration.

We define a p h + ( λ ) as the normalized spectrum of phytoplankton absorption coefficient spectrum with respect to its value at the reference wavelength, which can be expressed as:

where λ 0 is the reference wavelength. Eq. (13) means one can obtain the spectral shape by any given phytoplankton absorption or specific absorption spectrum. Because we will afterwards obtain a p h ( λ ) once a p h ( λ 0 ) is provided. This value is deemed as the spectral magnitude and can be obtained based on its relationship with [Chl]. The relationship is often described using a power law function:

where A is the scale factor, and the power exponent, E, is often observed to be slightly below one in global data sets (Bricaud et al., 1998; Woźniak and Dera, 2007). E represents the effects of pigment packaging and its interaction with phytoplankton cell size. Combining Eqs. (13) and (14), a p h ( λ ) of specific phytoplankton groups can be spectrally determined.

The reference wavelength λ 0 is often set at 443 nm (Bricaud et al., 2004). However, in this study, it was set at 676 nm, a wavelength where chlorophyll a is the dominating pigment absorption, and where a p h * is assumed to be very similar between different taxonomic groups, while variations at 443 nm would be larger due to a different pigment composition in each group.

Having c p h ( λ ) = [ Chl ] c p h * ( λ ) , the Chl-specific attenuation coefficient of phytoplankton, c p h * ( λ ) , is expressed as the sum of those of the different groups:

The phytoplankton Chl-specific scattering spectrum, b p h * ( λ ) , is then taken from the difference between those for attenuation and absorption:

The backscattering of phytoplankton, b b p h , is calculated by multiplying b p h ( λ ) with the weighted sum of the backscattering probabilities of each phytoplankton group, b ˜ b p h , which is expressed as:

For phytoplankton groups for which c p h * ( λ ) in Eq. (16) is unavailable, we will use ω p h ( λ 0 ) and re-normalized c p h ( λ ) at the reference wavelength, i.e., c p h + ( λ ) / c p h + ( λ 0 ) , to determine c p h * ( λ ) :

where λ₀ is set at 676 nm likewise and ω p h ( λ 0 ) can be empirically obtained from the measurements or literature for specific phytoplankton groups, which controls the magnitude of c p h * ( λ ) , and c p h + ( λ ) only serves the spectral shape and can be assumed safely as a power law function ( λ 0 / λ ) γ p h where the exponent γ p h determines the spectral slope. The power law function was also used to extrapolate c p h * ( λ ) beyond the measured wavelength range.

While the “Four-term” IOP model is highly adaptable, it may not be the most efficient option for simulating Case-1 waters, which make up the majority of the Earth’s oceans and are influenced by only two principal components: pure water and phytoplankton. Although the “Four-term” model can simulate these waters, it often requires more computational time and the gain in accuracy is small. To improve efficiency, we added a “Two-term” IOP model, as a supplement, by following the data synthesis process outlined in the IOCCG (2006), but using the setups for pure water and phytoplankton from the “Four-term” model (section 2.1) to mimic different phytoplankton groups. In the “Two-term” model, the fraction of phytoplankton groups from pico-phytoplankton to micro-phytoplankton is constrained by varying limits based on [Chl] levels. At low [Chl] levels, oceanic groups dominate, while other groups become mixed in as [Chl] increases (Brewin et al., 2010; Losa et al., 2017). However, in the “Four-term” model, there are no bounded limits as phytoplankton communities may be not correlated with [Chl] in optically complex water such as coastal and inland waters. The fraction setting is a preliminary biological limit that prevents unrealistic scenarios, but can be further refined with additional knowledge inputs.

This study evaluated the performance of an IOP model by comparing its ability to reproduce co-measured concentrations and IOP results using component concentrations as input. The aim was to demonstrate the capability to accurately describe the bio-optical process. To evaluate the performance, several statistical metrics including linear regression were utilized. The squared correlation coefficient, R ², measures how well a model fits the observed data. The slope, S, represents the change in the measured value for each unit increase in the estimated values, while the intercept, I, represents the systematic offset when comparing measurements and estimations. Better fits are indicated by R ² and S values closer to one and I value closer to zero. Additionally, the unbiased median percentage difference, D, was used to evaluate the accuracy of IOP models, with lower values indicating better accuracy. D was calculated as the median of ( X − Y ) / ( X + Y ) × 200 % , where X and Y are measurement and estimation values. The measurement and estimation values are scaled using a logarithmic base 10 transformation prior to the calculation of these statistical metrics.

Recall that the two-term model takes only [Chl] as input, while the four-term model takes [Chl], [ISM], and a g ( 440 ) as inputs. The R18 model uses [Chl], [TSM], and a g ( 440 ) , while the L23M model takes only a p h ( λ ) . All of these models output spectral IOPs of the components. In this section, we compared these benchmark models in terms of a p h , a d , a g , a g p , b p , and b b p , focusing specifically on the results at 440 nm. Scatter plots for the Hereon data set and the external data sets are in Figures 6 , 7 , respectively. Additional scatter plots at other wavelengths are provided in the Supplementary Material ( Figures A3 , A4 : 560 and 670 nm). Figure 8 shows spectral analysis to evaluate model performance from 400 to 700 nm. Multispectral samples were excluded in the spectral analysis to avoid outliers, primarily from OC-CCI v3, resulting in different valid estimation numbers (N) in scatter plots and spectra plots, but this does not significantly affect the spectral trend.

Figure 6 shows the reproduced IOPs for the Hereon data from our model in comparison with two other models, whereby for our model, the “Two-term” setup was applied to the oceanic data in Hereon and the “Four-term” setup to coastal and river data. The comparison between the two- and four-term results of the “oceanic” samples indicates no significant deviation, which is crucial for the seamless modeling of IOPs during the transition from Case-1 to Case-2 waters. Figure 6A demonstrated accurate a p h ( 440 ) estimates with a slope close to one, which verifies the reliability of the a p h * ( λ ) determination and the phytoplankton group spectral shapes, as the result at 440 nm can be regarded as being extrapolated from 676 nm. The percentage difference in a p h from 400 to 700 nm, shown in Figure 8 (A7), was within approximately ±5%, with some results within ±10%. The evaluation of a d ( 440 ) in Figure 6D revealed good performance across a range of four orders of magnitude, especially when a d ( 440 ) was greater than 0.1 m^-1, where the particulate matter was dominated by the inorganic component. However, the a d ( 440 ) estimation was more scattered for lower values, possibly due to increased variability in the a b d * ( λ ) shape due to increased organic components and uncertainties in concentration inputs caused by general measurement errors for organically dominated [TSM], hence [ISM]. The spectral statistical results of a d in Figure 8 (A8) exhibited relatively consistent values across the wavelength range, with percentage differences between ±5%. The a g ( 440 ) results, shown in Figure 6G , were in line with the “Four-term” input, while the “Two-term” results were dependent on the random parameter values. Despite the wide variation range, the error lines covered the measured values. The statistical results of a g in Figure 8 (A9) showed percentage differences within ±5% up to 600 nm, but degradation was observed beyond 600 nm as the values approached zero. Our IOP model performed well in the evaluation of a g p ( 440 ) , the sum of a p h ( 440 ) , a d ( 440 ) , and a g ( 440 ) , across different water types and over wavelengths, with no outliers observed. The b p ( 440 ) estimates, shown in Figure 6M , were accurate across three orders of magnitude, but with a few outliers due to higher measurement variation of [TSM]. The b b p ( 440 ) results were close to the 1-to-1 line for values greater than 0.01 m^-1, but with some deviations in the low-value region, especially for the “Two-term” model in the oceanic data. The spectral distribution of statistical values for b p and b b p in Figure 8A remained consistent, with percentage differences within ±5%.

The R18 model generally yielded comparable results at 440 nm to our model, albeit with some differences for certain parameters. Specifically, R18 model underestimated a d at longer wavelengths, as shown in Figure 8 (A8), with a slope of 0.78 lower than that of our model (0.93) at 560 nm, and with 0.71 lower than 0.99 at 670 nm. This difference could be due to the imperfect separation of biogenic and minerogenic parts of detritus in R18, as the pigment residuals in longer wavelengths, e.g., ≥ 560 nm, were still evident in their a b d * ( λ ) , leading to the underestimation of a m d * ( λ ) in these wavelengths. Additionally, consistent discrepancies between measurements and estimations of b p ( λ ) and b b p ( λ ) were observed across wavelengths, which were also found in other data sets, as shown in Figure 8B . Based on the points of Hereon – Ocean and some points of Hereon – Coast resembling oceanic waters, as shown in Figure 6 , the L23M model demonstrates favorable results for Case-1 waters when a p h ( λ ) spectra are provided. Conversely, the model generally underestimates all IOPs for Case-2 waters across the wavelength range concerned, as demonstrated in Figure 8B . This outcome is anticipated since the model is primarily designed for global oceanic waters. The superior performance of the L23M model for Hereon – Ocean is primarily due to the fact that a p h ( λ ) is its primary input, which mainly governs the total IOPs.

Figure 7 displays an independent evaluation of IOP models based on in situ data from external sources. Our model performs better in estimating a p h ( 440 ) for various water types, including CDOM-dominated lakes from M17, eutrophic lakes from C22, and oceanic waters from OC-CCI v3, as evidenced by the lower percentage difference. In addition, our model produces a larger number of estimated a p h ( 440 ) values, which may interfere with the reliability of other statistical metrics if other models are unable to estimate these points. However, the a d ( 440 ) results are slightly underestimated in C22 and overestimated M17 in both our model and R18. The higher [ISM] values observed in M17 compared to similar samples in the Hereon data may explain the overestimation of a d ( 440 ) , which can also be observed from the overestimation of b p ( λ ) in Figures 7M, N . The overestimation of a d ( λ ) was slightly reduced by L23M in Figure 7F , when only [Chl] was used to model a d ( λ ) . One should consider the uncertainty associated with component concentration when using IOP models (IOCCG, 2019). Despite this, the overall performance of a g p ( 440 ) remains strong by our model, with an R² of 0.83, even in the presence of more variability in OC-CCI v3. The spectral percentage difference of a g p mostly fell within ±10%, as shown in Figure 8 (B10). The b p ( 440 ) estimates by our model are slightly underestimated for extremely turbid water (mainly from HYPERMAQ). For such highly turbid waters, the uncertainty of in situ measurements should be considered (IOCCG, 2019). Only a few of points of b b p ( λ ) are observed due to the limited data from these external data sets. The different performance between our model and R18 in estimating b b p ( 440 ) may result from the use of similar measurement devices in the studies of R18 and M17. The reason for the absence of results for Hereon – Ocean and OC-CCI v3 in R18 is the unavailability of the input parameter [TSM] in that specific part of the data. Note that the smaller data size and larger measurement uncertainties in the external data sets compared to the Hereon data can impact statistical results in Figure 8B . Nonetheless, our model still demonstrates consistent and robust performance across the spectrum, with most components falling within a ±10% difference range.

To summarize, our IOP model shows comparable performance as other the two models in terms of reproducing IOPs from concentrations in specific water environments and in terms of absolute values and spectral shapes. However, our model provides a higher amount of reproducibility across diverse optical water types, encompassing oceanic, coastal, and inland waters. Our spectral statistical metrics demonstrate that the estimations of model have a percentage difference of within ±5% for most component IOPs throughout the visible spectrum, with some falling within ±10%, which is essential for meeting the accuracy requirements outline in GCOS (2011) for developing ocean color algorithms. Another advantage of our IOP model is its flexibility, allowing for replacement of specific model parameters with more optimal values in special scenarios. Although these parameters defined so far have proven to be representative, this enables continuous improvement and refinement of the model over time.

In the previous section, we established that our IOP model can accurately reproduces water optical properties based on component concentrations. However, we also need to determine whether the model can capture most of the optical variability, which depends on the concentration range and the dynamic range of model parameters like γ and ω . We admit that our current ranges may not encompass some extreme scenarios. For instance, to address the underestimation depicted in Figure 7M , we increased the original value of ω d ( 550 ) and reduced γ d to simulate a highly turbid water type. By adjusting model parameters such as γ and ω , the resulting ranges offer valuable insights into the appropriate direction for fine-tuning the model. Although these modifications improved that particular case, they lay outside 96% of the training data distribution. Despite this, we used the original model parameters established in section 3.2 and the concentration range provided in Table 4 to do the simulation as a showcase, which is used to test the data coverage with other data sets.

The data distribution presented in Figure 9 compares in situ measurements and simulated data sets. Panels (A) and (C) show our simulated data and in situ measurements, while panels (B) and (D) compare other simulated data sets. Our selection of 560 nm was intended to showcase the data distribution, but it is worth noting that the results for other wavelengths are also comparable. The a p ( 560 ) and ω p ( 560 ) relationships demonstrate the ratio of b p to c p as particle absorption increases. The a p ( 560 ) values of Hereon vary up to four orders of magnitude, with C22 and HYPERMAQ having higher values, and M17 varying over a smaller range with lower values. Despite the variation in a p ( 560 ) , Figure 9A demonstrates that the variation in  ω p ( 560 ) is relatively stable, with measured values as low as about 0.9 and simulated values as low as 0.8. As anticipated, the “Four-term” model contributes to most of the value variation. Figure 9B indicates that other simulated data sets have similarities, with the L23 data set having a comparable distribution to our simulations, albeit with fewer points where a p ( 560 ) is greater than 0.1. The simulated CCRR data set shows a wider distribution in a p ( 560 ) but stable ω p ( 560 ) , primarily due to a constant γ d and a [Chl]-related γ p h . However, the simulated C2X data set displays more pronounced differences, with some points having significantly lower ω p ( 560 ) values for high a p ( 560 ) and phytoplankton-dominated water, suggesting that only 20% of the light is scattered for a given high particle concentration, resulting in low reflectance, which is inconsistent with real-world observations. The primary reason for such low ω p ( 560 ) values is the inappropriate parameterization of b p h ( λ ) in C2X. These points have lower b p h ( λ ) values with relatively higher a p h ( λ ) (not shown). Figures 9C, D show the relationship between a g p ( 560 ) and b b p ( 560 ) , with b b p ( 560 ) increasing as a g p ( 560 ) increases in the measured data, suggesting a stable proportion of a g ( 560 ) in a g p ( 560 ) that co-varies with [Chl], as reflected in CCRR and L23. C2X contains simulations of extremely absorbing waters, resulting in simulated points below the distribution of other data due to high CDOM absorption. Our simulated data, using a similar concentration limit from C2X, has a wider coverage compared to the C2X distribution, with most points originating from the “Four-term” model (CDOM is independent of [Chl]), while C2X has some data gaps.

For practical reasons, the study limited the number of IOP components to a manageable amount instead of including all particle species in water (Stramski et al., 2004). In the “Four-term” model, a g was modeled using a few spectra (N=12) with similar shapes ( S g = 0.0174 nm^-1 ± 0.1%). The resulting model was able to accurately reproduce the CDOM spectrum with precise a g ( 440 ) values. In the “Two-term” model, a g was correlated to [Chl] with a random value, resulting in relatively greater variability, but still covering the range of values observed in oceanic waters. The shape of specific detritus absorption, a d * , varies across different aquatic environments but is generally limited on a local scale (Woźniak and Dera, 2007). The slope of minerogenic detritus absorption, S m d * , in our study, 0.0104 nm^-1 in Figure 2B , was consistent with values observed in other natural water basins (Woźniak and Dera, 2007). The combination of biogenic and minerogenic detritus in our model enhances its applicability across different regions, which is important for simulating IOPs of complex waters on a global scale. The use of different quantiles (quasi-organic proportion of the total particle) to define a b d * aligns with the observed increase in biogenic detritus absorption when organic matter dominates. The slopes of our simulated a d values range from 0.0068 to 0.015 (as per the definition in Stramski et al. (2019), based on the natural logarithmic transformation at 440 and 550 nm), which are generally consistent with previous studies, e.g., 0.0095~0.0148 in Zheng et al. (2015) and 0.0031~0.0169 in Stramski et al. (2019), but lack the “flat” shapes (lower values) found in Stramski et al. (2019) due to limited data in the training set, which can be expanded in the future to incorporate such extreme scenarios. The IOP model proposed in this study includes seven phytoplankton groups that are used as general “color” groups, rather than specific phytoplankton species or functional types that have different absorption and scattering properties. By using these general groups, the model can capture the increased variability in phytoplankton IOPs that is visually seen in natural systems. In other words, the phytoplankton groups in our model are not based on specific biological characteristics, but rather on their broad spectral characteristics, allowing for more flexibility and accuracy in predicting phytoplankton IOPs. Variables a p h ( 412 ) / a p h ( 443 ) and a p h ( 555 ) / a p h ( 490 ) , used to describe the shape of a p h in Zheng et al. (2015), were found to be in good agreement with other in situ and simulated data sets, within the ranges of 0.8~1.0 and 0.2~1.0, respectively. Although some ranges may not be covered, it is not necessary to add more phytoplankton groups if no color differences are added. Our IOP model is also able to simulate coccolithophore blooms with a 490 nm reflectance peak by multiplying its absorption coefficient by 0.02, which is biologically based on the vanishing of the absorbing algae cells but remaining high scattering from detached coccoliths (Neukermans and Fournier, 2018; Cazzaniga et al., 2021).

Complete coverage of the optical diversity of as many natural waters as possible is also important for the exploitation of optical water type classification methods. Existing methods (e.g., Moore et al., 2014; Hieronymi et al., 2017; Jackson et al., 2017; Bi et al., 2021) can only partially cover the optical variability; many complex waters are out-of-scope for these methods, which is reflected in too low class memberships (not shown in detail here). This also applies to the classifiability of R r s from various ocean color atmospheric corrections (Hieronymi et al., 2023). Both aspects identify potential weaknesses of the OWT methods, which should be revised in the future, e.g., by means of our IOP model. It is noteworthy that the analysis here did not account for variability when modeling IOP to R r s , and we simplified backscattering as a spectral constant rather than a sophisticated phase function or a volume scattering function. This could impact the R r s angular distribution in many complex waters (Harmel et al., 2016). When analyzing OWT frameworks, environmental factors, and sun-viewing geometries should also be considered and varied (Bi et al., 2023).

Figure 10 demonstrates that the disparity in R r s between modeled b p h ( λ ) and b p h GM 83 ( λ ) was more noticeable at high [Chl] (>100 mg/m³), with b p h GM 83 ( λ ) showing lower values (up to 0.02 sr^-1 difference) due to the modeled b p h ( 550 ) discrepancy. At low [Chl] (< 1 mg/m³), b p h GM 83 ( λ ) displayed slightly higher R r s values (up to 0.003 sr^-1 or 40% difference). The most significant differences were observed in the green and NIR (540~590 nm and 700~730 nm), with differences of up to 90%, while the blue region (400~450 nm) was the least affected as non-phytoplankton signals prevail. Renormalizing b p h GM 83 ( λ ) to b p h Adj ( λ ) , which makes b p h Adj ( 560 ) equal to b p h ( 560 ) , significantly reduced differences in R r s except for higher [Chl] (>100 mg/m³). However, the values in NIR were still lower than those with b p h ( λ ) by up to 0.004 sr^-1 (or 35%). Figure 11 presents six intuitive samples of IOPs and R r s ordering by increasing [Chl] values. The use of b p h GM 83 ( λ ) underestimated b p ( λ ) from concentrations evaluated in the training and testing data sets (not shown) and did not fit the spectral features of phytoplankton observed in the measured b p data (Castagna et al., 2022; Röttgers et al., 2023), resulting in suppressed R r s values in the green and longer wavelengths for high [Chl]. This contradicts the high NIR reflectance commonly seen during algal blooms (Hieronymi et al., 2017; Bi et al., 2019; Castagna et al., 2022). The differences were significantly mitigated by using b p h Adj ( λ ) , although the shape difference still remains.

The sensitivity analysis results, on the forward modeling aspect, revealed that the assumption of b p h shape has a limited impact on R r s modeling, provided that the b p h magnitude is reasonable. This finding may explain why some inversion models perform well, despite using a simplified power law function for scattering shape. However, it is worth noting that the simplified scattering shape may not accurately represent the scattering spectra and may lead to variance transference between scattering and absorption in the non-linear optimizing approach, as discussed in Roesler and Boss (2003) from the perspective of inversion modeling. Furthermore, a power law function may not be sufficient for modeling optically complex water, particularly at high biomass, where scattering dominates over absorption in water color (McLeroy-Etheridge and Roesler, 1998).

It is important to note that adjusting b p h to b p h Adj , as presented in this section, is an ideal scenario that may not always be feasible due to limited knowledge of b p h magnitude in many practical cases. As discussed about the C2X in Figure 9B , modeling absorption and scattering of phytoplankton (detritus as well) separately can lead to significant uncertainties when the input concentrations are outside their training ranges. To address this issue, a relative IOP parameter, the single-scattering albedo, is used in our new model, which has been proven to improve accuracy. Additionally, attenuation, unlike scattering, is relatively simpler to mathematically describe (Twardowski et al., 2001; Roesler and Boss, 2003). Therefore, we suggest an improved parameterization of phytoplankton IOPs in terms of both their spectral shape and magnitude, which has been shown to be applicable for various natural water types.