A synthetic optical database generated by radiative transfer simulations in support of studies in ocean optics and optical remote sensing of the global ocean

Radiative transfer (RT) simulations have long been used to study the relationships between the inherent optical properties (IOPs) of seawater and light fields within and leaving the ocean, from which ocean apparent optical properties (AOPs) can be calculated. For example, inverse models used to estimate IOPs from ocean color radiometric measurements have been developed and validated using the results of RT simulations. Here we describe the development of a new synthetic optical database based on hyperspectral RT simulations across the spectral range of near-ultraviolet to near-infrared performed with the HydroLight radiative transfer code. The key component of this development is the generation of a synthetic dataset of seawater IOPs that serves as input to RT simulations. Compared to similar developments of optical databases in the past, the present dataset of IOPs is characterized by the probability distributions of IOPs that are consistent with global distributions representative of vast areas of open-ocean pelagic environments and coastal regions, covering a broad range of optical water types. The generation of synthetic data of IOPs associated with particulate and dissolved constituents of seawater was driven largely by an extensive set of field measurements of the phytoplankton absorption coefficient collected in diverse oceanic environments. Overall, the synthetic IOP dataset consists of 3320 combinations of IOPs. Additionally, the pure seawater IOPs were assumed following recent recommendations. The RT simulations were performed using 3320 combinations of input IOPs, assuming vertical homogeneity within an infinitely deep ocean. These input IOPs were used in three simulation scenarios associated with assumptions about inelastic radiative processes in the water column (not considered in previous synthetically generated optical databases) and three simulation scenarios associated with the sun zenith angle. Specifically, the simulations were made assuming no inelastic processes, the presence of Raman scattering by water molecules, and the presence of both Raman scattering and fluorescence of chlorophyll a pigment. Fluorescence of colored dissolved organic matter was omitted from all simulations. For each of these three simulation scenarios, the simulations were made for three sun zenith angles of 0, 30, and 60^∘ assuming clear skies, standard atmosphere, and a wind speed of 5 m s⁻¹. Thus, overall 29 880 RT simulations were performed. The output results of these simulations include radiance distributions, plane and scalar irradiances, and a whole set of AOPs, including remote-sensing reflectance, vertical diffuse attenuation coefficients, and mean cosines, where all optical variables are reported in the spectral range of 350 to 750 nm at 5 nm intervals for different depths between the sea surface and 50 m. The consistency of this new synthetic database has been assessed through comparisons with in situ data and previously developed empirical relationships involving IOPs and AOPs. The database is available at the Dryad open-access repository of research data (https://doi.org/10.6076/D1630T, Loisel et al., 2023).

1 Introduction

Investigating the propagation of natural light in the ocean can be addressed experimentally through in situ measurements and theoretically through numerical radiative transfer (RT) simulations. Understanding the relationships between the radiometric quantities (i.e., radiance and irradiances) that characterize the light fields within and leaving the ocean and the inherent optical properties (IOPs) of the water column, as well as boundary conditions at the sea surface (i.e., surface illumination conditions and sea state) and at the ocean bottom (i.e., bottom depth and albedo), requires comprehensive datasets of multiple variables acquired over a broad range of environmental conditions. For example, of particular interest are the relationships between the spectral remote-sensing reflectance of the ocean (in sr⁻¹), R_rs(λ), which is an apparent optical property (AOP) derivable from radiometric quantities, and the seawater IOPs that are directly linked to various seawater constituents because these relationships form the cornerstone of various applications of optical (ocean color) remote sensing. Recent technological developments and broader accessibility of optical in situ instrumentation have led to a significant increase in optical datasets collected across diverse oceanic environments, and efforts have been undertaken to merge data from various sources within publicly available databases (e.g., Werdell and Bailey, 2005; Valente et al., 2019; Casey et al., 2020). Although the importance of field data collection across diverse environments cannot be overstated, existing database compilations are subject to certain limitations. In addition to typical measurement errors, it is difficult to ensure consistent data quality and characterization of uncertainties across all merged data because individual datasets are often obtained with different instruments as well as measurement and data processing methods. Also, even large databases such as NASA's SeaWiFS Bio-optical Archive and Storage System (SeaBASS, https://seabass.gsfc.nasa.gov/, last access: 1 March 2023) cannot ensure the balanced representativeness of collected field data in terms of a broad range of optical conditions across diverse ocean environments. In this context, radiative transfer (RT) simulations, which are free of measurement errors, provide a useful tool to generate comprehensive synthetic databases and complement the existing datasets of field measurements in support of studies in ocean optics and optical remote sensing.

Over the past decades, various radiative transfer models that employ different numerical solution techniques have been developed and used to address a wide range of problems related to the optics of natural water bodies (e.g., Mobley et al., 1993; Mobley, 1994; Stamnes et al., 2017). Since the 1990s, the HydroLight code based on the invariant imbedding technique (Mobley, 1989; Mobley et al., 1993; Mobley, 1994) has been among the most commonly used radiative transfer models in oceanographic optics. The HydroLight code solves the scalar (i.e., polarization of light is not included) time-independent radiative transfer equation for a horizontally homogeneous water body in which the IOPs can vary with depth and under given boundary conditions at the surface and at the bottom of the water body. Inelastic radiative processes within the water column that include Raman scattering by water molecules, fluorescence of chlorophyll a pigment, and fluorescence of colored dissolved organic matter (CDOM) can be included in HydroLight simulations.

Radiative transfer simulations with HydroLight code have proven useful for generating synthetic databases of light field characteristics (i.e., radiance and irradiances) within and leaving the ocean and the AOPs derived from the simulated radiometric quantities for various scenarios of seawater IOPs that provide input to the simulations. In particular, as a result of efforts dedicated to inverse bio-optical algorithms and coordinated under the auspices of the International Ocean Colour Coordinating Group (IOCCG Report, 2006), a widely used publicly available synthetic database was generated within the spectral range 400 to 800 nm with a 10 nm resolution for clear sky conditions with three different sun zenith angles (0, 30, and 60^∘), a sea-surface state corresponding to a wind speed of 5 m s⁻¹, and 500 different IOP combinations driven by chlorophyll a concentration, Chla (in units of mg m⁻³), within the surface ocean layer. The input IOP data included the spectral absorption coefficients of phytoplankton, a_ph(λ), non-algal particles (also referred to as depigmented or detrital particles, which can include various types of particles such as organic detritus, mineral particles, heterotrophic bacteria, and depigmented phytoplankton cells), a_d(λ), colored dissolved organic matter (CDOM), a_g(λ), and the spectral backscattering coefficients of phytoplankton, b_b-ph(λ), and non-algal particles, b_b-d(λ) (λ represents the wavelength of light in a vacuum expressed in units of nm and the IOP coefficients are typically expressed in units of m⁻¹). The output parameters provided by simulations available in the public database included the following AOPs: spectral remote-sensing reflectance, R_rs(λ), remote-sensing reflectance just below the sea surface, r_rs(λ), irradiance reflectance just below the sea surface, $R(z={\mathrm{0}}^{-},\mathit{\lambda })$, and the diffuse attenuation coefficient for downwelling plane irradiance, K_d(λ,z), at the depths $z={\mathrm{0}}^{-}$, 5, and 10 m (where 0⁻ indicates a depth of just beneath the sea surface).

Another synthetic database that is publicly available was developed as part of the CoastColour Round Robin project (Nechad et al., 2015). This project was focused on coastal waters and IOPs were described by 5000 combinations of Chla, a_g(λ), and the mass concentration of mineral particles. The HydroLight simulations were run from 350 to 900 nm at 5 nm intervals for a cloudless sky, three sun zenith angles (0, 40, and 60^∘), and a wind speed of 5 m s⁻¹. The output parameters included in the publicly available database are the water-leaving reflectance, RL_w(λ)=πR_rs(λ), K_d(λ), photosynthetically available radiation, PAR, and the euphotic depth, z_eu. Most recently, a synthetic database was also developed by the first NASA PACE (Plankton, Aerosol, Cloud, ocean Ecosystem) Science Team, where the ocean contribution to top-of-atmosphere radiances were simulated by HydroLight (Craig et al., 2020). These simulations were performed from 350 to 800 nm with a 5 nm step for a cloudless sky, three sun zenith angles (10, 30, and 60^∘), a wind speed of 5 m s⁻¹, and a set of 720 IOP combinations driven by a_ph(λ). The publicly available output of these HydroLight simulations is R_rs(λ).

While these existing synthetic databases have provided valuable information to the ocean color radiometry (OCR) community, especially for the purpose of algorithm development where the ocean AOPs are linked to IOPs, there are several reasons that have motivated the present study, which is aimed at generating a new optical synthetic database. First, inelastic Raman scattering and fluorescence processes were ignored in the previous RT simulations. These inelastic radiative processes are known to be important for simulating realistic characteristics of light fields within and leaving the ocean, including R_rs(λ), which is a primary optical quantity used in ocean color remote sensing. For example, Raman scattering by water molecules may have an important influence on light within and leaving the ocean and AOPs, especially in the green and red parts of the spectrum (e.g., Marshall and Smith, 1990; Stavn, 1993; Sugihara et al., 1984; Westberry et al., 2013). Second, the three synthetic databases described above are based on the use of the spectral absorption, a_w(λ), and scattering, b_w(λ), coefficients of pure seawater in the visible part of the spectrum as defined by Pope and Fry (1997) and Morel (1974), respectively. However, more recent measurements and theoretical considerations provide new recommendations for spectral values of a_w(λ) and b_w(λ) (IOCCG Protocol Series, 2018; Zhang and Hu, 2009; Zhang et al., 2009). Third, the probability distributions of different IOPs that were used as input to previous RT simulations do not appear to match well with the IOP distributions observed in extensive field datasets or satellite-derived datasets representing the global ocean. This issue may have a biasing effect when the synthetic databases are used to develop optical algorithms based on AOP vs. IOP relationships, especially when the underlying goal is to represent a broad range of IOPs encountered within the global ocean, even if the primary interest is in open-ocean pelagic environments. Finally, the previous synthetic databases were developed specifically for OCR-oriented studies and publicly accessible data generally include only the surface reflectances, R_rs(λ), R(λ), r_rs(λ), and K_d(λ), at selected depths. These databases do not include many of the various output variables obtained from RT simulations, such as various underwater AOPs, which can be useful in supporting a broader range of studies in ocean optics beyond ocean color remote sensing.

In this article, we present a new synthetic optical database generated using RT simulations that addresses some of the limitations of similar databases developed in the past. First, we describe the development of the synthetic IOP dataset that is required to run RT simulations. The key roles in this development are played by the measured data of the phytoplankton absorption coefficient and the desired consistency between the probability distributions of synthetic IOPs and the global distributions based on satellite observations. Following this, we describe different configurations of RT simulations that were performed with the HydroLight code. The next section is dedicated to consistency between the new optical synthetic database and in situ data, including some previously reported empirical relationships. We provide example illustrations of consistency for both the IOP and AOP data. The closing section summarizes the structure of synthetic database files and provides an example illustration of one output radiometric variable, the spectral downwelling plane irradiance, calculated with RT simulations.

2 Development of a synthetic dataset of seawater inherent optical properties

2.1 General overview of methodology

The scope of the synthetic database generated with RT simulations and the degree of its representativeness of diverse marine optical environments within the global ocean depend most critically on the dataset of seawater IOPs used as input to RT simulations. In the present study, our approach to generate the IOP dataset was driven largely by the underlying goal to obtain the probability distributions of IOPs that are generally consistent with the distributions observed in the global ocean dominated by open-ocean pelagic environments. The key IOPs involved in the creation of our IOP dataset include the spectral absorption and backscattering coefficients associated with the main categories of seawater constituents representing suspended particulate matter and CDOM. Specifically, the absorption coefficients of the different constituents are the spectral absorption coefficients of phytoplankton, a_ph(λ), non-algal particles, a_d(λ), and CDOM, a_g(λ). Note that the sum $a_{\mathrm{ph}}\left(\mathit{\lambda }\right)+a_{\mathrm{d}}\left(\mathit{\lambda }\right)=a_{\mathrm{p}}\left(\mathit{\lambda }\right)$ represents the particulate absorption coefficient with combined contributions of phytoplankton and non-algal particles, and the sum $a_{\mathrm{d}}\left(\mathit{\lambda }\right)+a_{\mathrm{g}}\left(\mathit{\lambda }\right)=a_{\mathrm{dg}}\left(\mathit{\lambda }\right)$ represents the non-phytoplankton absorption coefficient with combined contributions of non-algal particles and CDOM. The backscattering coefficient of the different constituents are the spectral backscattering coefficients of phytoplankton, b_b-ph(λ), and non-algal particles, b_b-d(λ), such that the sum $b_{\text{b-ph}}\left(\mathit{\lambda }\right)+b_{\text{b-d}}\left(\mathit{\lambda }\right)=b_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ is the particulate backscattering coefficient.

Among these constituent IOPs, the phytoplankton absorption coefficient, a_ph(λ), plays the most fundamental role in the creation of the synthetic dataset of IOPs in this study. The a_ph(λ) spectra in this dataset were derived from actual measurements of phytoplankton absorption made on near-surface samples collected across diverse oceanic environments. Thus, the a_ph(λ) data are not “synthetic” in the sense that these data were not obtained from a modeling approach, although some spectral interpolation or extrapolation was applied to measured data as described in more detail below. In contrast, the remaining four constituent IOPs in the IOP dataset, i.e., a_d(λ), a_g(λ), b_b-ph(λ), and b_b-d(λ), are “synthetic” in the sense that they are entirely based on calculations using a modeling approach with some assumptions about the magnitude and spectral behavior of the modeled IOPs. Importantly, the measured values of a_ph(λ) were used in the calculations of these IOPs. These calculations are also described in detail below. Thus, each combination of the five constituent IOPs in the synthetic IOP dataset consists of the measured a_ph(λ) and the calculated a_d(λ), a_g(λ), b_b-ph(λ), and b_b-d(λ), where the results of these calculations depend on the measured a_ph(λ). As a result of this approach, it would seem justifiable to refer to the created IOP dataset as a quasi-synthetic dataset. For simplicity, however, we refer to it as a synthetic IOP dataset while bearing in mind that a_ph(λ) spectra were derived from measurements.

2.2 Description of the in situ dataset

Figure 1a depicts the location of oceanographic stations where the near-surface measurements of a_ph(λ) were made. As shown, these measurements were collected across diverse open-ocean and coastal environments. Their total number is 4382, constituting the initial field dataset of a_ph(λ) considered in this study. Figure 1 also shows the location of stations where coincident measurements are available for the pairs of IOP coefficients, namely a_ph(λ) and a_g(λ) (Fig. 1b), a_ph(λ) and a_dg(λ) (Fig. 1c), and a_ph(λ) and b_bp(λ) (Fig. 1d). We recall that while the in situ data of a_g(λ), a_dg(λ), and b_bp(λ) were not used in the development of the synthetic IOP dataset, they were assembled for the purpose of comparison with corresponding coefficients that were calculated and included in the synthetic IOP dataset. Many in situ data of the IOP coefficients used in the present study were collected in previous studies (e.g., Reynolds et al., 2001; Babin et al., 2003; Loisel et al., 2007; Claustre et al., 2008; Huot et al., 2008; Stramski et al., 2008; Lubac et al., 2008; Loisel et al., 2009; Bricaud et al., 2010; Loisel et al., 2011; Antoine et al., 2011; Neukermans et al., 2012; Uitz et al., 2015; Neukermans et al., 2016; Reynolds et al., 2016; Aurin et al., 2018; Reynolds and Stramski, 2019; Stramski et al., 2019). Some data are described in publications devoted to compilation of various datasets (Valente et al., 2019; Casey et al., 2020) and are included in several databases (e.g., SeaBASS, CoastlOOC, BOUSSOLE, and GOCAD). As the IOP coefficients in the in situ dataset were measured over a broad range of trophic and environmental conditions, their spectral values span more than 3 or 4 orders of magnitude. This large dynamic range is illustrated in terms of probability distributions at selected light wavelengths, i.e., at 440 nm for the constituent absorption coefficients and 550 nm for the particulate backscattering coefficient (Fig. 2).

Figure 1 Locations of oceanographic stations where in situ measurements were collected for (a) a_ph(λ), number of measurements N=4382; (b) a_ph(λ) and a_g(λ), number of matchup measurements N=2206; (c) a_ph(λ) and a_dg(λ), number of matchup measurements N=813; and (d) a_ph(λ) and b_bp(λ), number of matchup measurements N=775.

Figure 2 Histograms and relevant statistical parameters of field measurements of (a) a_ph(440), (b) a_g(440), (c) a_dg(440), and (d) b_bp(550). N is the number of measurements, and $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{̃}}{x}$ are the mean and median values of each IOP, respectively.

2.3 Generation of the dataset of hyperspectral a_ph(λ)

The first task necessary in developing the synthetic IOP dataset was to assemble data on the hyperspectral absorption coefficient of phytoplankton, a_ph(λ), from field measurements collected across diverse open-ocean and coastal environments (Figs. 1a, 2a). These a_ph(λ) data were obtained with the filter-pad spectrophotometric method as a difference between the measurements of a_p(λ) and a_d(λ) (Kishino et al., 1985; IOCCG Protocol Series, 2018). Historically, most of these measurements were acquired with the transmittance configuration of the filter-pad method; such measurements are included in our dataset. However, some data in our dataset were obtained with the inside integrating-sphere configuration of the filter-pad method, which is superior to the transmittance configuration of measurement (Stramski et al., 2015; IOCCG Protocol Series, 2018).

A significant portion (23.7 %) of the initial dataset of a_ph(λ), consisting of 4382 measurements, covers a spectral range of 400 to 750 nm with a high spectral resolution of data reported at 1 nm intervals. In some cases, the original measurements extended to the near-UV spectral region and/or longer wavelengths in the near-infrared spectral region (800 or 850 nm). The data beyond 750 nm are not used in this study because our RT simulations target the spectral range of 350 to 750 nm. It is notable that the absorption measurements of marine particles and phytoplankton are generally unavailable or are not reported in the UV because of increased methodological challenges and uncertainties in this spectral region (Stramski et al., 2015; IOCCG Protocol Series, 2018; Kostakis et al., 2021). As a result, only a relatively small fraction of a_ph(λ) measurements in the initial dataset were reported in the near-UV region. In addition, the initial dataset included a relatively large fraction of a_ph(λ) measurements that were reported at wavelength intervals larger than 1 nm. These lower-resolution data (hereafter referred to as multispectral) ranged from a small wavelength interval of 2 nm to data reported at more limited numbers of wavelengths (as few as <10) within the visible spectral range. It is likely that the multispectral data available from some data sources that we used in this study were originally measured at higher spectral resolution but were eventually reported only for some selected wavelengths, such as those corresponding to the spectral bands available on satellite ocean color sensors.

The first objective of the analysis of a_ph(λ) was to consider the initial a_ph(λ) dataset within the 400–750 nm range and convert the measurements that were reported at lower spectral resolution to uniformly hyperspectral data at 1 nm intervals. In this analysis, all measurements originally available at 1 nm intervals were considered to provide reference spectral shape functions of a_ph(λ). The originally multispectral data of a_ph(λ) were converted to hyperspectral data using several different approaches depending on the spectral features of lower-resolution data. One approach utilized the reference spectral shape functions of a_ph(λ) and was applied to multispectral a_ph(λ) data if they were reported for fewer than 100 wavelengths. In this case, a given multispectral spectrum of a_ph(λ) was converted to a hyperspectral spectrum using a specific hyperspectral measurement that exhibited the highest correlation with the multispectral measurement under consideration. The correlation coefficient was calculated using the spectral data available at common wavelengths of the considered pair of spectra. A necessary condition to proceed with conversion of a given multispectral spectrum to a hyperspectral spectrum was a correlation coefficient of 0.95 or higher. If this condition was satisfied, the multispectral data were converted to hyperspectral data so that the created hyperspectral spectrum maintained a multispectral measurement magnitude in the range of 440–450 nm and had the spectral shape of the reference hyperspectral measurement. An alternative approach to convert multispectral data to hyperspectral data involves a linear interpolation of multispectral data. This approach was used when the multispectral data were reported at relatively small wavelength intervals (at least 100 spectral data points available between 400 and 750 nm) or when the correlational analysis described above did not yield a correlation coefficient of 0.95 or higher (5.2 % of the multispectral data). The original multispectral spectra which did not include data below 450 nm or fell into the category of data subject to linear interpolation, but had no data above 700 nm, were rejected from further analysis. For all hyperspectral spectra that passed the above-described analysis and criteria (i.e., 2204 spectra: the 593 original hyperspectral measurements and 1611 hyperspectral spectra created from the multispectral data), the null-point correction was applied by subtracting the average value of a_ph(λ) in the 745–750 nm range from all spectral values in the 400–750 nm range.

The next step of the analysis was to extend all null-point corrected spectra of a_ph(λ) that cover the 400–750 nm range into the UV spectral region. The primary focus was on the 350–400 nm range because our RT simulations were designed to provide output results in the 350–750 nm range. For this purpose, we used a separate subset of reference hyperspectral measurements of a_ph(λ) that includes the near-UV spectral region. This reference subset of data consisted of 233 measurements collected across bio-optically diverse marine environments in the Pacific and Atlantic oceans and western Arctic seas. The majority of these 233 spectra (170) were collected with the inside integrating-sphere configuration of filter-pad method, while the remaining 63 measurements were done using either the transmittance or transmittance-reflectance filter-pad configuration (Zheng et al., 2014). A correlational analysis was applied to pairs of spectra, each consisting of a spectrum covering the 400–750 nm range and a reference spectrum covering the 350–750 nm range. The correlation coefficient was calculated using data at common wavelengths from the 400–750 nm range. The reference spectrum that yielded the highest correlation with the investigated 400–750 nm spectrum was selected as a basis for extrapolation of the investigated spectrum into the 350–400 nm range. This extrapolation ensured that a given investigated spectrum maintained its magnitude at 400 nm and that its extrapolated near-UV portion had the spectral shape of the selected reference spectrum. The final aspect of extrapolation in the UV is related to the spectral range 300–350 nm. The rationale for IOP data extending to 300 nm is to ensure that the results of RT simulations that start at 350 nm account for possible effects of Raman scattering by water molecules in the UV spectral region. Therefore, for the 300–350 nm range, we simply assumed that a_ph(λ) in this range is equal to a_ph(350). The limitation associated with this assumption is not considered to be serious given the limited role of the 300–350 nm range in the RT simulations and the weak Raman scattering effects in UV spectral region. Example spectra of a_ph(λ) in the 350–750 nm range from contrasting marine environments are presented in Fig. 3. These examples show significant variation in both the magnitude and spectral shape of a_ph(λ).

Figure 3 (a) Two example spectra of a_ph(λ) from contrasting oceanic environments. For each example a_ph(λ), two spectra are displayed, namely the measurement from the initial a_ph(λ) dataset shown at the original wavelength intervals (red points) and the spectrum after interpolation to 1 nm intervals (if required) and null-point correction (continuous lines). The UV portion of the latter was obtained by extrapolation based on reference data in the UV (see text for details). (b) Example of normalized a_ph(λ) spectra illustrating the variability of the spectral shape of the a_ph(λ) database. These spectra have been normalized to their integral.

2.4 Generation of the complete IOP dataset

In the next step of the analysis, the subset of 2204 a_ph(λ) spectra that was created from the initial a_ph(λ) dataset as described above was subject to additional modifications to ensure that the final a_ph(λ) dataset was characterized by the probability distribution that resembles the distribution representative of the global ocean. This process and background information on the motivation for such adjustments in the probability distribution are described below.

When the end goal is to achieve a high degree of representativeness of the global ocean as in this study, the process of assembling in situ datasets of IOPs is unavoidably subject to limitations, even if relatively large amounts of data from many field experiments and cruises are considered. This is mainly because the global ocean is dominated by vast areas of open-ocean pelagic environments and the amount of IOP data collected in these environments is disproportionally limited compared to the amount of data collected in coastal regions, which represent a relatively small portion of the global ocean. Thus, the probability distributions based on in situ datasets, such as those presented in Fig. 2, are expected to deviate from the probability distributions representative of the global ocean. In particular, the maxima of probability distributions and the measures of central tendency, such as the median and mean values, obtained from compilations of a relatively large amount of in situ IOP data (such as in Fig. 2) are expected to be shifted to larger values compared to actual global distributions because the IOPs exhibit a general tendency of higher values in coastal regions compared to open-ocean environments. While this issue has been recognized, it has not been addressed or resolved in various studies that focus on global ocean color applications. For example, current global ocean color algorithms for estimating chlorophyll a concentration (Chla) are based on relatively large amounts of in situ data whose probability distribution is shifted significantly to higher Chla compared with the global Chla distribution (O'Reilly and Werdell, 2019). Similarly, in the development of previous synthetic optical databases with RT simulations (e.g., IOCCG Report, 2006), no special attempt was made to ensure consistency between the probability distributions of input IOP data and the distributions expected for the global ocean. In the recent development of refined global ocean color algorithms for estimating the concentration of particulate organic carbon (POC), the in situ dataset was assembled with a goal to achieve reasonable consistency with a global POC distribution (Stramski et al., 2022). This goal was, however, achieved at the expense of a significant reduction in the amount of accepted in situ data compared to the size of the overall pool of available in situ data.

In this study, our goal was to create a relatively large synthetic IOP dataset based on the initial dataset of several thousand measurements of spectral a_ph(λ) so that the probability distributions of IOPs in the final synthetic dataset are reasonably consistent with the expected distributions representative of the global ocean. As described above, the initial field dataset in support of this process consisted of 4382 spectra of a_ph(λ); this number was further reduced to 2204 spectra that were accepted as a result of analysis and some criteria applied to the initial dataset. This reduced dataset of accepted a_ph(λ) spectra was then further modified to ensure that the final probability distribution of a_ph(440) resembles the global distribution of a_ph(440). The global probability distribution of a_ph(440) was estimated using retrievals of a_ph(440) from satellite ocean color data. Specifically, we used global satellite observations made with the ocean color sensor OLCI (Ocean and Land Colour Instrument) deployed on the Sentinel-3 mission (Donlon et al., 2012) from the period 1 December 2020 through 30 November 2021. The weekly data product of remote-sensing reflectance R_rs(λ) at 4 km² spatial resolution was used as input to the three-step semi-analytical algorithm (3SAA) to derive a_ph(443) as described in Jorge et al. (2021). The a_dg(443) and b_bp(λ) coefficients were also derived from this algorithm. In general, the 3SAA first derives the diffuse attenuation coefficient for downwelling plane irradiance averaged within the surface layer down to the first attenuation depth, 〈K_d(λ)〉₁, from R_rs(λ), and then utilizes the inverse model LS2 (Loisel et al., 2018) to derive the total absorption, a(λ), and backscattering, b_b(λ), coefficients from R_rs(λ) and 〈K_d(λ)〉₁. After subtracting the pure seawater contributions, the non-water absorption, a_nw(λ), and the particulate backscattering, b_bp(λ), coefficients are obtained. Finally, a_ph(λ) and a_dg(λ) are derived from a_nw(λ) using an optimization algorithm of Zhang et al. (2015) with modifications that account for differences in optical water types defined in terms of different spectral shapes of R_rs(λ) (Mélin and Vantrepotte, 2015). While the original classification of Mélin and Vantrepotte (2015) includes 16 optical water classes (OWCs), the derivation of a_ph(λ) and a_dg(λ) from the 3SAA additionally included a 17th OWC to improve the representation of ultra-oligotrophic waters such as those found in the South Pacific Gyre (Morel et al., 2007; Claustre et al., 2008; Stramski et al., 2008) and in some areas of the Mediterranean Sea in summer (Loisel et al., 2011). This 17th OWC is described in Jorge et al. (2021).

The 3SAA does not yield the separate contributions of CDOM, a_g(λ), and non-algal particles, a_d(λ), to the overall non-phytoplankton absorption coefficient, a_dg(λ). Therefore, we also used another semi-analytical model (CDOM-KD2) described in Bonelli et al. (2021) to estimate a_g(443) from OLCI-derived R_rs(λ). Having a_dg(λ) from the 3SAA and a_g(λ) from CDOM-KD2, non-algal particulate absorption, a_d(λ), was obtained as the difference a_dg(λ)−a_g(λ). As a result of this analysis, we obtained a dataset of satellite-derived constituent absorption coefficients, a_ph(443), a_g(443), a_d(443), and a_dg(443), as well as the particulate backscattering coefficient, b_bp(550), where we focused on the spectral band near 440 nm for absorption and 550 nm for backscattering.

For illustrative purposes, Fig. 4a depicts the spatial distribution of 17 optical water classes (OWCs) over the global ocean obtained from satellite OLCI data following the methodology of Mélin and Vantrepotte (2015). For further illustrative purposes, these 17 OWCs were grouped into three optical water groups (OWGs). Group 1 consists of OWC1 and OWC2, which are characterized by high water turbidity such as in coastal areas affected by discharge from large rivers. Although the focus of this study is to create synthetic datasets representative primarily of open ocean and moderately turbid coastal waters, an explicit identification of Group 1 data that represent very turbid waters is of interest for comparisons with the database developed specifically for coastal waters by Nechad et al. (2015). The second OWG, Group 2, includes six OWCs: from OWC3 through OWC8. This group represents mainly productive waters in both coastal and open-ocean environments, such as those encountered in the North Atlantic during phytoplankton bloom periods (Lévy et al., 2005). Finally, Group 3 included the remaining nine OWCs: from OWC9 through OWC17. These water types are observed mainly in mesotrophic and oligotrophic regions of the global ocean. Based on this classification, 79.6 % of OLCI water pixels in Fig. 4a belong to Group 3, 10.8 % to Group 2, and 9.6 % to Group 1. The histograms of OLCI-derived R_rs(443) associated with these three groups of data are shown in Fig. 4b. For comparative purposes we also assembled a dataset of in situ measurements of R_rs(λ), which were collected at various locations within the global ocean (Fig. 4c). The histograms of in situ R_rs(443) associated with Groups 1, 2, and 3 are depicted in Fig. 4d, which show a similar pattern to that in Fig. 4b. For the in situ dataset of R_rs(λ), 69.2 % of data belong to Group 3, 15.7 % to Group 2, and 15.1 % to Group 1.

Figure 4 (a) Global map illustrating the distribution of 17 optical water classes estimated from monthly R_rs(λ) values derived from satellite observations with the ocean color sensor OLCI from December 2020 through November 2021 (weekly products at 4 km²). The color bar scale refers to optical water classes. (b) Histogram of OLCI-derived R_rs(443) for the three optical water groups (see text for details). (c) Location of oceanographic stations where in situ measurements of R_rs(λ) were collected and used to analyze the consistency of the synthetic dataset with field measurements. (d) Histograms of in situ measurements of R_rs(443) for the three optical water groups.

The probability density functions (PDFs) of global satellite-derived a_ph(440), a_g(440), a_dg(440), and b_bp(550) are depicted in Fig. 5. We refer here to satellite-derived absorption coefficients at 440 nm although they were derived from OLCI reflectances at 443 nm, which is a minor difference that is inconsequential for the purpose of this study. Comparison of Figs. 2a and 5a indicates that the distribution of measured a_ph(440) from our initial field dataset (Fig. 2a) is shifted towards higher values compared to the global distribution of satellite-derived a_ph(440) (Fig. 5a). The probability distribution of reduced dataset of measured a_ph(440) (N=2204) that was created from the initial field dataset of a_ph(λ) shows similar deviations from the global distribution (not shown). Thus, to create the final dataset of a_ph(λ) that has the probability distribution of a_ph(440) consistent with the global satellite-derived distribution, we adjusted the number of a_ph(440) measurements in each bin of the histogram of the reduced dataset either by removing the measurements from any given bin or by adding the measurements to this bin. The removal or addition of a_ph(440) measurements associated with any given bin was done by subjecting all a_ph(440) measurements originally contained within a given bin to random selection. Specifically, in the case of addition the randomly selected a_ph(440) was added as a replicate of a_ph(440) to a given bin. In the case of removal, the randomly selected a_ph(440) was simply removed from a given bin. As a result of this process, we obtained a modified distribution of measured a_ph(440) that is fairly consistent with the satellite-derived distribution of a_ph(440). Both the modified histogram and the corresponding modified PDF of measured a_ph(440) are depicted in Fig. 5a for comparison with the global satellite-derived distribution. In total, this modified distribution consists of 3320 measurements of a_ph(440); obviously, each of these measurements at 440 nm has an associated full spectrum of a_ph(λ) values between 300 and 750 nm. These 3320 spectra of a_ph(λ) represent one IOP component of the final synthetic IOP dataset.

Figure 5 Histograms showing the distribution of the synthetic IOP data used in the present study. The synthetic and satellite-derived probability density functions (PDFs) for each IOP are represented by the solid and dashed curves, respectively.

The full synthetic IOP dataset created in this study consists of 3320 combinations of measured a_ph(λ) and synthetically generated a_d(λ), a_g(λ), b_b-ph(λ), and b_b-d(λ). Below is a description of calculations of a_g(λ), a_d(λ), b_b-ph(λ), and b_b-d(λ). We note that all IOP coefficients are expressed in units of m⁻¹ and the light wavelength is given in nanometers.

The four IOP coefficients were calculated using a similar methodology to that applied in previous studies aiming at generation of synthetic ocean optical databases (IOCCG Report, 2006; Craig et al., 2020). Specifically, we used the measured values of a_ph(440) as the main driver of calculations of a_g(λ), a_d(λ), b_b-ph(λ), and b_b-d(λ). Thus, the variability in the measured a_ph(440), as depicted by the probability distribution of measured a_ph(440) in Fig. 5a, is the main source of variability in these four co-existing IOP coefficients. It is notable that the replicate values of a_ph(440) present within any given bin of the a_ph(440) distribution result in the generation of different values of the four IOP coefficients because the formulas involved in these calculations contain random parameters. The coupling between a_ph(440) and CDOM absorption coefficient was defined as

$$\begin{array}{}\text{(1)}& a_{\mathrm{g}}\left(\mathrm{440}\right)={\mathrm{10}}^{(P_{\mathrm{1}}+\mathit{\gamma })},\end{array}$$

where P₁ is a parameter related to a_ph(440) and γ is randomly selected from a predetermined range of values (Table 1). The spectral values of a_g(λ) are subsequently determined from

$$\begin{array}{}\text{(2)}& a_{\mathrm{g}}\left(\mathit{\lambda }\right)=a_{\mathrm{g}}\left(\mathrm{440}\right)e^{-S_{\mathrm{g}}\left(\mathit{\lambda }-\mathrm{440}\right)},\end{array}$$

where the spectral slope parameter, S_g in units of nm⁻¹, is randomly selected from a predetermined range of values (Table 1). The absorption coefficient of non-algal particles was modeled in a similar fashion:

$$\begin{array}{}\text{(3)}& {\displaystyle }{a}_{\mathrm{d}}\left(\mathrm{440}\right)=P_{\mathrm{2}}{a}_{\mathrm{ph}}\left(\mathrm{440}\right),\text{(4)}& {\displaystyle }{a}_{\mathrm{d}}\left(\mathit{\lambda }\right)=a_{\mathrm{d}}\left(\mathrm{440}\right)e^{-S_{\mathrm{d}}\left(\mathit{\lambda }-\mathrm{440}\right)},\end{array}$$

where P₂ is a parameter related to a_ph(440) and the spectral slope parameter S_d (nm⁻¹) is randomly selected from a predetermined range of values (Table 1). The parameterizations of P₁ and P₂ were chosen to match relationships observed in the in situ dataset assembled in this study.

Table 1 Symbols of variables, mathematical expressions, and corresponding equations in the text of the paper. rng(0,1) is a random number between 0 and 1.

Particulate backscattering is not modeled in terms of the single coefficient b_bp(λ) but instead as separate contributions by phytoplankton, b_b-ph(λ), and non-algal particles, b_b-d(λ), so that their sum yields b_bp(λ). In order to calculate b_b-ph(λ), first the formula that couples a_ph(440) with the beam attenuation coefficient of phytoplankton at 550 nm, c_ph(550), is used:

$$\begin{array}{}\text{(5)}& c_{\mathrm{ph}}\left(\mathrm{550}\right)=P_{\mathrm{3}}{\mathrm{Chla}}^{\mathrm{0.57}}=P_{\mathrm{3}}{\left[{\displaystyle \frac{a_{\mathrm{ph}}\left(\mathrm{440}\right)}{\mathrm{0.05582}}}\right]}^{\mathrm{0.57}},\end{array}$$

where Chla is the concentration of chlorophyll a in units of mg m⁻³, 0.05582 (m² mg⁻¹) is the value of chlorophyll-specific absorption coefficient of phytoplankton at 440 nm, $a_{\mathrm{ph}}^{\ast }\left(\mathrm{440}\right)$ (Maritorena et al., 2002), and P₃ is a parameter with a randomly selected value from a predetermined range (Table 1). The exponent value of 0.57 is based on the study of Voss (1992). Subsequently, the spectral values of phytoplankton beam attenuation coefficient are calculated from

$$\begin{array}{}\text{(6)}& c_{\mathrm{ph}}\left(\mathit{\lambda }\right)=c_{\mathrm{ph}}\left(\mathrm{550}\right){\left({\displaystyle \frac{\mathrm{550}}{\mathit{\lambda }}}\right)}^{S_{\text{c-ph}}},\end{array}$$

where the spectral slope parameter, S_c-ph (dimensionless), is calculated using both a_ph(440) and a random number generator (Table 1). Next, the spectral scattering coefficient of phytoplankton is determined:

$$\begin{array}{}\text{(7)}& b_{\mathrm{ph}}\left(\mathit{\lambda }\right)=c_{\mathrm{ph}}\left(\mathit{\lambda }\right)-a_{\mathrm{ph}}\left(\mathit{\lambda }\right),\end{array}$$

where the spectral values of a_ph(λ) are from the same measured spectrum as the value of a_ph(440) in Eq. (5). Finally, the spectral backscattering coefficient of phytoplankton is calculated from

$$\begin{array}{}\text{(8)}& b_{\text{b-ph}}\left(\mathit{\lambda }\right)=\mathrm{0.01}{b}_{\mathrm{ph}}\left(\mathit{\lambda }\right),\end{array}$$

where 0.01 is the value of the backscattering ratio of phytoplankton, ${\stackrel{\mathrm{̃}}{b}}_{\text{b-ph}}$, assumed to be constant and independent of the light wavelength in the present study (IOCCG, 2006; Loisel et al., 2007). Laboratory measurements performed on various phytoplankton cultures have shown, however, that ${\stackrel{\mathrm{̃}}{b}}_{\text{b-ph}}$ can exhibit a slight spectral variation, with the value at 442 nm ranging from 0.0035 to 0.029 (Whitmire et al., 2010). We note that b_b-ph(λ) is not required as input to our radiative transfer simulations but b_ph(λ) is needed.

To calculate the backscattering coefficient of non-algal particles, b_b-d(λ), phytoplankton absorption at 440 nm is first coupled with the scattering coefficient of non-algal particles at 550 nm, b_d(550), using the following relationship:

$$\begin{array}{}\text{(9)}& b_{\mathrm{d}}\left(\mathrm{550}\right)=P_{\mathrm{4}}{\mathrm{Chla}}^{\mathrm{0.766}}=P_{\mathrm{4}}{\left[{\displaystyle \frac{a_{\mathrm{ph}}\left(\mathrm{440}\right)}{\mathrm{0.05582}}}\right]}^{\mathrm{0.766}},\end{array}$$

where the parameter P₄ is randomly selected from a predetermined range (Table 1) and the value of 0.05582 is $a_{\mathrm{ph}}^{\ast }\left(\mathrm{440}\right)$ as explained in relation to Eq. (5). The exponent value of 0.766 is based on the study of Loisel and Morel (1998). Then, the spectral values of non-algal scattering coefficient are calculated from

$$\begin{array}{}\text{(10)}& b_{\mathrm{d}}\left(\mathit{\lambda }\right)=b_{\mathrm{d}}\left(\mathrm{550}\right){\left({\displaystyle \frac{\mathrm{550}}{\mathit{\lambda }}}\right)}^{S_{\text{b-d}}},\end{array}$$

where the spectral slope parameter, S_b-d (dimensionless), is calculated using both a_ph(440) and a random number generator (Table 1). In the final step, the spectral backscattering coefficient of non-algal particles is calculated as

$$\begin{array}{}\text{(11)}& b_{\text{b-d}}\left(\mathit{\lambda }\right)=\mathrm{0.018}{b}_{\mathrm{d}}\left(\mathit{\lambda }\right),\end{array}$$

where the constant 0.018 is the backscattering ratio of non-algal particles, ${\stackrel{\mathrm{̃}}{b}}_{\text{b-d}}$. This value was proposed by Mobley (1994) and was derived by averaging three particle phase functions measured in oceanic waters by Petzold (1972). Again, we note that b_b-d(λ) is not required as input to radiative transfer simulations but b_d(λ) is needed. The spectral slope of b_bp(λ), γ, where b_bp(λ) is obtained as the sum of b_b-ph(λ) and b_b-d(λ), has a mean and standard deviation of −1.10 ± 0.34 and exhibits a trend from oligotrophic (where γ is around −2) to eutrophic waters (where the b_bp(λ) spectrum is nearly flat). These results are in good agreement with previous studies (Morel and Maritorena, 2001; Loisel et al., 2006; Antoine et al., 2011).

The variability of measured a_ph(440) illustrated in Fig. 5a, along with the dynamic range of parameters P₁, P₂, P₃, P₄, the spectral slopes S_g, S_d, S_c-ph, and S_b-d, and the degree of randomness in the selection of these parameters for any given value of a_ph(440) that initiates the process of calculating a_g(λ), a_d(λ), b_b-ph(λ), and b_b-d(λ), resulted in the generation of the synthetic dataset of IOP coefficients covering a wide dynamic range consistent with in situ and satellite observations over the global ocean. Figure 5b, c, and d compare the probability distributions of satellite-derived a_g(440), a_dg(440), and b_bp(550) with the distribution of these coefficients from the final synthetic IOP dataset. This comparison supports the general consistency of the distributions of these IOP coefficients, which is in line with the desired consistency achieved for a_ph(440) (Fig. 5a) as discussed earlier in this section. It is also noteworthy that in contrast to this newly created synthetic IOP dataset, the previous synthetic datasets exhibit significant differences between the probability distributions of synthetic IOPs and global distributions based on satellite observations (Fig. 6).

Figure 6 Histograms showing the distribution of IOPs from the synthetic datasets of the IOCCG Report (2006) and Craig et al. (2020) in green and pink, respectively. The IOP distributions estimated from satellite ocean color observations by the OLCI sensor over the global ocean are represented by the black line.

Overall, the above-described synthetic IOP dataset includes 3320 scenarios of non-water IOPs, i.e., IOPs associated with the variable contributions of phytoplankton, non-algal particles, and CDOM to the optical properties of seawater. In addition to the non-water absorption coefficients, a_ph(λ), a_d(λ), and a_g(λ), as well as the non-water scattering coefficients, b_ph(λ) and b_d(λ), the radiative transfer simulations required input of scattering phase functions of particles, specifically for phytoplankton and non-algal particles. We assumed the particulate phase functions proposed by Fournier-Forand (1994) with the backscattering ratio ${\stackrel{\mathrm{̃}}{b}}_{\text{b-ph}}=\mathrm{0.01}$ for phytoplankton and ${\stackrel{\mathrm{̃}}{b}}_{\text{b-d}}=\mathrm{0.018}$ for non-algal particles. Note that while the backscattering ratios are assumed to be spectrally constant, the phase functions vary with light wavelength because of spectral variations of b_ph(λ) and b_d(λ). All IOP data in the final synthetic IOP dataset cover the spectral range of 300 to 750 nm with a 5 nm interval. This wavelength interval is consistent with the intended output of our radiative transfer simulations.

The radiative transfer simulations also required input of the absorption and scattering properties of pure seawater. For the spectral absorption coefficient of pure seawater, a_w(λ), we used the values recommended in the IOCCG Protocol Series (2018). This recommendation includes the values from Jonasz and Fournier (2007) for the spectral range 300–330 nm, Morel et al. (2007) for 340–415 nm, Pope and Fry (1997) for 420–725 nm, and Kou et al. (1993) for 730–750 nm. The spectral volume scattering function of pure seawater (from which the spectral scattering coefficient and scattering phase function can be obtained) was calculated following Zhang et al. (2009), assuming a water temperature of 18 ^∘C and a salinity of 35 ‰. The temperature of 18 ^∘C is consistent with the mean sea-surface temperature (SST) calculated from the monthly global NOAAv2 SST database at 1^∘ spatial resolution from December 1991 through November 2021 (Jérôme Vialard, personal communication, 2023). The salinity of 35 ‰ is also consistent with the global surface average (Durack et al., 2013).

3 Radiative transfer simulations

The IOP dataset described in Sect. 2, which includes 3320 combinations of non-water IOPs, provided the key input to radiative transfer (RT) simulations that were performed with the HydroLight v5.0 radiative transfer code (Mobley and Sundman, 2008). All RT simulations were run assuming vertically homogeneous IOPs within the water column and an infinitely deep ocean, i.e., no seafloor effect on the light field within the water column. For all simulations, the computed radiometric and AOP variables were saved into the output data files at 10 cm depth intervals between the ocean surface and 1 m depth, and at 1 m intervals between 1 and 50 m depth. Thus, the primary focus of our RT simulations is on the ocean surface layer that can potentially contribute to light leaving the ocean and which has significance to remote sensing by spaceborne or airborne optical instruments. All simulations were carried out in the spectral range of 300 to 750 nm using 5 nm spectral bands, and the results were produced for the nominal wavelengths of each of the 81 bands, that is at 350, 355, 360, …, 745, 750 nm. The results in the 300–350 nm range were not retained in the output files (that include seawater IOPs, radiometric quantities, and AOPs) because this spectral region was included primarily to account for the potential effects of inelastic processes at wavelengths longer than 350 nm and, additionally, because it is known that the uncertainties in the characterization of seawater IOPs can increase significantly at wavelengths shorter than 350 nm.

For 3320 scenarios of input IOPs, we performed several separate sets of RT simulations that differed in terms of assumed sea-surface boundary conditions and the inclusion or exclusion of inelastic radiative processes within the water column. The assumptions regarding the sea-surface boundary conditions were the same as in the previous RT simulations described in Loisel et al. (2018). Specifically, all simulations were made under the same assumption of a wind speed of 5 m s⁻¹, which determines the sea-surface roughness involved in the calculations of transmission and reflection of light at the air–water interface. In all simulations the sky conditions were also assumed to be the same, i.e., clear skies and standard atmosphere. However, three distinct sets of simulations were made for the three sun zenith angle values of 0, 30, and 60^∘. With regard to consideration of inelastic processes, we also performed three distinct sets of simulations. The first of these sets assumed the absence of inelastic processes in water, that is no Raman scattering by water molecules, no fluorescence by chlorophyll a, and no fluorescence by CDOM. The second set of these simulations included Raman scattering by water molecules. Finally, the third set included both Raman scattering and chlorophyll a fluorescence, and this scenario of inelastic processes is expected to generally provide the most realistic simulations of radiative transfer in the ocean surface layer. We note, however, that fluorescence by CDOM was not included in any simulations. The Raman scattering coefficient, phase function, and wavelength distribution function were set to the default values described in HydroLight technical documentation (Mobley, 2012). The quantum efficiency of chlorophyll a fluorescence, which may exhibit significant variability (nearly five-fold, between about 0.01 and 0.05) in ocean waters (Maritorena et al., 2000; Morrison et al., 2003), was also set to its default value of 0.02 in the HydroLight code. For each scenario of sun zenith angle and inelastic processes, we performed 3320 RT simulations, each for a different combination of seawater IOPs. Thus, given the three sun zenith angles, the three scenarios of inelastic processes, and 3320 combinations of IOPs, overall we performed 29 880 simulations. The combination of the synthetic IOP dataset used as input to RT simulations (Sect. 2) and the results of the radiance, other radiometric quantities, and AOPs obtained from these 29 880 simulations (described in this section) constitute the synthetic ocean optical database developed in this study.

4 Comparisons of the synthetic database with in situ data

In this section, we compare the selected spectral IOP coefficients from the synthetic IOP dataset with the in situ data of the IOPs and the selected spectral AOPs from the synthetic database generated with the RT simulations with the in situ data of the AOPs. In these comparisons, we also include some empirical relationships between the IOPs and between the AOPs that were established in previous studies based on the analysis of in situ data.

Table 2 Symbols, variables, and units for the various quantities included in the final synthetic optical database.

^* All optical variables in the database are spectral and provided at different light wavelengths between 350 and 750 nm at 5 nm intervals and different depths within the water column between the sea surface and the 50 m depth, except for R_rs, and L_w, which are defined at the sea surface.

Figure 7 depicts scatter plots of IOP coefficients, specifically a_g(440) vs. a_ph(440) (Fig. 7a), a_dg(440) vs. a_ph(440) (Fig. 7b), a_d(440) vs. a_ph(440) (Fig. 7c), and b_bp(550) vs. a_ph(440) (Fig. 7d). The scatter plots include two datasets, the in situ dataset and the synthetic dataset, as described in Sect. 2. We recall that in both types of datasets, a_ph(440) plotted on the x axis is the same because the phytoplankton absorption data used in this study were obtained from field measurements with no modeling involved. The scatter plots show a significant degree of overlap which indicates general consistency between the synthetic and in situ datasets. Similar patterns are observed when the $a_{\mathrm{g}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, $a_{\mathrm{dg}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, $a_{\mathrm{d}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, and $b_{\mathrm{bp}}\left(\mathrm{550}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$ ratios are plotted as a function of a_ph(440) (not shown). For illustrative purposes, the data from the synthetic IOP dataset are color coded to indicate the partitioning of data into the three OWGs, i.e., Groups 1, 2, and 3, which were defined using the synthetic spectra of R_rs(λ) generated through RT simulations with input of the synthetic IOP data. As expected, the data with the generally lowest values of IOPs belong to Group 3, the data with intermediate values of IOPs to Group 2, and the data with the highest IOPs (most turbid waters) to Group 1. We also note that the in situ dataset exhibits a somewhat wider dynamic range of variability than the synthetic dataset, especially when the IOP ratios, $a_{\mathrm{g}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, $a_{\mathrm{dg}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, $a_{\mathrm{d}}\left(\mathrm{440}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, and $b_{\mathrm{bp}}\left(\mathrm{550}\right)/a_{\mathrm{ph}}\left(\mathrm{440}\right)$, are relatively high. While this result can reflect some degree of intrinsic difference in the dynamic range covered by the two datasets, it must also be recognized that some variability in the in situ dataset may be associated with the fact that these data were collected on numerous cruises by different groups of investigators using methodologies (instrumentation, data processing, data quality control, etc.) that are unavoidably different across the different data sources.

Figure 7 (a) a_g(443), (b) a_dg(443), (c) a_d(443), and (d) b_bp(550) as a function of a_ph(443) for the in situ dataset (black data points) and the synthetic dataset (colored data points). The black polygon lines in each panel delimit approximately the scatter of the in situ data points (black dots). Each color refers to the optical water group as indicated (139, 262, and 2919 data points for Group 1, 2, and 3, respectively). Empirical relationships previously developed for (a) a_g(443) vs. a_ph(443) and (d) b_bp(550) vs. a_ph(443) are also displayed for comparison. The original relationships were formulated as a function of Chla and the presented relationships were obtained by converting Chla to a_ph(443) using the chlorophyll-specific phytoplankton absorption at 443 nm from Bricaud et al. (1998).

For additional comparative purposes, Fig. 7a and d include a few empirical relationships between the IOPs in question, which were established in previous studies based on considerable amount of field measurements collected mostly in open-ocean environments. As seen, the relationships between a_g(440) and a_ph(440) based on the studies of Morel (2009) and Bricaud et al. (2010) agree quite well with the central tendency of variation within our synthetic dataset. We note that Morel (2009) and Bricaud et al. (2010) reported relationships between a_g(440) and Chla that were very similar in these studies. For the purpose of illustration in our Fig. 7a, we replaced Chla with a_ph(440) using the formula Chla =a_ph(440)$/\mathrm{0.05582}$. Similarly, the studies of Huot et al. (2008) and Antoine et al. (2011) reported on empirical relationships between b_bp(λ) and Chla. After converting Chla to a_ph(440) as mentioned above, these two relationship are plotted in Fig. 7d. Although these two relationships have different slopes, they are both generally consistent with the average trend of variation in the synthetic dataset.

Radiative transfer is driven mainly by two ratios of IOPs, which are the scattering to absorption ratio, $b\left(\mathit{\lambda }\right)/a\left(\mathit{\lambda }\right)$, which controls the number of scattering events (Morel and Gentili, 1991), and the molecular to total scattering ratio, $b_{\mathrm{w}}\left(\mathit{\lambda }\right)/b\left(\mathit{\lambda }\right)$, which is the parameter controlling the weighted sum of the particle scattering and molecular scattering phase functions (Morel and Loisel, 1998; Loisel and Stramski, 2000). Figure 8 shows the distribution of these two ratios at 440 nm for the synthetic dataset. The $b_{\mathrm{w}}\left(\mathrm{440}\right)/b\left(\mathrm{440}\right)$ and $b\left(\mathrm{440}\right)/a\left(\mathrm{440}\right)$ ratios range between about 0 and 0.2 and between 0.5 and 10, respectively, which is consistent with previous models developed for Case 1 waters (Figs. 2 and 3 in Morel and Gentili, 1991; Fig. 2 in Morel and Loisel, 1998).

Figure 8 (a) Histograms of (a) $b_{\mathrm{w}}\left(\mathrm{440}\right)/b\left(\mathrm{440}\right)$ and (b) $b\left(\mathrm{440}\right)/a\left(\mathrm{440}\right)$ for the synthetic dataset.

For comparing the AOPs from the synthetic database with in situ data, we have chosen two AOPs, the spectral remote-sensing reflectance, R_rs(λ), and the spectral diffuse attenuation coefficient of downwelling plane irradiance averaged within the water column from the sea surface to the first attenuation depth, 〈K_d(λ)〉₁, as well as the maximum band ratio (MBR) of reflectance. The scatter plot of our synthetic data of R_rs(555) vs. R_rs(443) is depicted in Fig. 9a. For comparison, the range of in situ data is illustrated by the dashed contour lines. The maximum value of R_rs(443) reached 0.0165 sr⁻¹, which is in good agreement with in situ measurements performed in ultra-oligotrophic waters in the South Pacific Gyre during the BIOSOPE cruise (see Fig. 3 in Stramski et al., 2008). These results are once again illustrated using color coding to represent different optical water types, specifically Groups 1, 2, and 3. As seen, there is relatively good agreement between the synthetic data and the range of variability of the in situ data for Groups 2 and 3 (Fig. 9a). For Group 1 (very turbid waters), however, the synthetic data exhibit a smaller range of variability compared with in situ data. This result is not unexpected because our primary goal was to generate the synthetic database that is most representative of open-ocean pelagic environments as well as coastal areas where water turbidity is low to moderate rather than very high. As described in Sect. 2, turbid waters of Group 1 correspond to OWCs 1 and 2 as defined in Mélin and Vantrepotte (2015). It is interesting to note that the synthetic optical database that was developed by Nechad et al. (2015) for coastal waters shows relatively good consistency between the synthetic and in situ data for Group 1 (Fig. 8b). However, in contrast to our synthetic database, the synthetic data of Nechad et al. (2015) exhibit a limited range of variability compared with in situ data for Groups 2 and 3. Thus, the synthetic data of Nechad et al. (2015) for turbid waters in Group 1 can provide useful complementarity to our synthetic database whose main focus is on water types from Groups 2 and 3.

Figure 9 (a) R_rs(555) as a function of R_rs(443) for the synthetic dataset (colored data points) and in situ dataset (colored contours). (b) Same as panel (a) but for the synthetic dataset of Nechad et al. (2015) that was developed for coastal waters. The color coding refers to the optical water groups as indicated.

The scatter plot of the synthetic data of 〈K_d(490)〉₁ as a function of the blue-to-green band ratio of reflectance, $R_{\mathrm{rs}}\left(\mathrm{490}\right)/R_{\mathrm{rs}}\left(\mathrm{555}\right)$, is shown in Fig. 10. These synthetic data are again color coded according to optical water classes defined in terms of Groups, 1, 2, and 3. For comparison, a few empirical relationships between these AOP variables established in previous analyses of field measurements are also displayed in Fig. 10 (Mueller, 2000; Werdell, 2005; Werdell, 2009). The relationship of Mueller (2000) was formulated during the early phase of the SeaWiFS satellite mission to serve as an operational global algorithm for estimating K_d(490) from ocean color observations. Werdell (2005) provided an updated relationship with a primary goal to improve the estimation of K_d(490) at low values of K_d(490) that correspond to high reflectance band ratio values. Figure 10 shows that these two relationships are generally consistent with our synthetic data across the entire range of variability encompassing data from Groups, 1, 2, and 3. This is reassuring given that the main purpose of our synthetic database and these two empirical relationships is similar in the sense of targeting the optical variability within the global ocean dominated by open-ocean environments. Figure 10 also includes the relationship of Werdell (2009) that represents the most recent update of global empirical algorithms for estimating K_d(490) from different ocean color satellite sensors. Specifically, the relationship of Werdell (2009) presented in Fig. 9 is referred to as KD2S and is based on SeaWiFS spectral bands. In contrast to the relationships of Mueller (2000) and Werdell (2005), the relationship of Werdell (2009) deviates significantly from our synthetic data within the range of relatively high values of 〈K_d(490)〉₁, which correspond to relatively low values of $R_{\mathrm{rs}}\left(\mathrm{490}\right)/R_{\mathrm{rs}}\left(\mathrm{555}\right)$. It is remarkable that this deviation occurs within the range where our synthetic data are classified as Group 1, so these optical water types are associated with high water turbidity. Another remarkable result illustrated in Fig. 10 is that the relationship of Werdell (2009) in this range is quite consistent with the main trend observed within the synthetic database of Nechad et al. (2015), which was developed for coastal environments. This result further supports the potential complementarity between our synthetic database and that of Nechad et al. (2015).

Figure 10 Scatter plot of K_d(490) vs. the blue-to-green reflectance ratio, $R_{\mathrm{rs}}\left(\mathrm{490}\right)/R_{\mathrm{rs}}\left(\mathrm{555}\right)$, for the synthetic database. The red, green, and blue data points represent OWGs 1, 2, and 3, respectively. The black cross data points are from the Nechad et al. (2015) synthetic dataset. The curves represent the relationships developed by Mueller (2000), Werdell (2005), and Werdell (2009). The K_d(490) data points represent 〈K_d(490)〉₁ for the present synthetic database (colored data points) and the near-surface K_d(490) calculated within the top 1 cm layer for the Nechad dataset (black data points).

The scatter plot of Chla vs. the MBR of reflectance for the synthetic database is shown in Fig. 11. The monotonically decreasing trend of Chla with increasing MBR is consistent with the SeaWiFS-specific OC5 algorithm for estimating Chla from MBR (O'Reilly and Werdell, 2019). For this illustration, we estimated chl a using the relationship between a_ph(660) and Chla from Bricaud et al. (1998), which is unavoidably affected to some extent by natural variability in this relationship.

Figure 11 Scatter plot of Chla vs. the blue-to-green maximum band ratio (MBR) of remote-sensing reflectance (i.e, $R_{\mathrm{rs}}\left(\mathrm{412}\mathit{>}\mathrm{443}\mathit{>}\mathrm{490}\mathit{>}\mathrm{510}\right)/R_{\mathrm{rs}}\left(\mathrm{555}\right)$) for the synthetic database. The red, green, and blue data points represent OWG 1, 2, and 3, respectively. The solid black line represents the OC5 algorithm developed by O'Reilly and Werdell (2019) for SeaWiFS spectral bands. For this illustration, Chla was calculated from a_ph(660) using the chlorophyll-specific phytoplankton absorption at 660 nm from Bricaud et al. (1998).

Figure 12 Examples of depth profiles of E_d(z,λ) for a given IOP scenario from (a) OWG 1, (b) OWG 2, and (c) OWG 3. Radiative transfer simulations were performed for a sun zenith angle of 30^∘ and included Raman scattering by water molecules and chlorophyll a fluorescence. The black line depicts the first optical attenuation depth, z₁

5 Data availability

The synthetic optical database described in this study is publicly available at the Dryad open-access repository of research data (Loisel et al., 2023; https://doi.org/10.6076/D1630T).

6 Summary

We have generated a new synthetic database that consists of seawater IOPs as well as corresponding radiometric quantities and AOPs within the ocean surface layer down to a depth of 50 m and at the sea surface. The radiometric quantities and AOPs were obtained from radiative transfer (RT) simulations performed with HydroLight code using the IOPs as input to the calculations. The list of variables included in the database is provided in Table 2. Because of the use of the absorption and scattering properties of pure seawater (assuming the salinity of 35 ‰) in the simulations, the present database cannot be applied to freshwater environments and special caution should be exercised for applications where water salinity is significantly less than 35 ‰ because of the decrease in pure seawater scattering. This database is organized following an easy-to-read NetCDF structure and divided into two subsets of data for which the file name identifies the sun zenith angle and the RT simulation scenario related to the presence or absence of inelastic radiative processes within the water column. The first subset of data includes the seawater spectral absorption and backscattering coefficients as well as sea-surface radiometric quantities relevant to ocean color radiometry, R_rs(λ), L_w(λ), $E_{\mathrm{d}}(z={\mathrm{0}}^{+}$, λ), and $L_{\mathrm{u}}(z={\mathrm{0}}^{+}$, λ), where $z={\mathrm{0}}^{+}$ is just above the surface. The surface and depth profile values of several spectral radiometric quantities and AOPs, as well as PAR, are included in the second subset of data. The spectra of z_eu and z₁ are also provided in the second file. More details on the organization and content of the database are included in the README file provided in the database.

In closing, we present an example illustration of one of the radiometric variables included in the output data files generated by RT simulations. We recall that the primary result of HydroLight simulations is the spectral radiance that provides comprehensive information about the angular distribution of the light field, from which different irradiances and AOPs are calculated. However, it is the spectral downwelling plane irradiance, E_d(z,λ), that has been the most commonly measured radiometric quantity in ocean optics, so in Fig. 12 we have chosen to illustrate HydroLight-simulated E_d(z,λ) within the ocean surface layer down to a depth of 50 m. These results are presented for three different scenarios of IOPs which are representative of three different optical water types defined in terms of Groups 1, 2, and 3 (see Sect. 2). These RT simulations were performed for the sun zenith angle of 30^∘ in the presence of Raman scattering by water molecules and chlorophyll a fluorescence in the water column. In addition to significant differences in the variation of the spectral E_d(λ) as a function of depth z between Groups 1, 2, and 3, Fig. 12 also illustrates distinct differences in the magnitude and spectral behavior of the first optical attenuation depth, z₁. This quantity is equivalent to the inverse of the diffuse attenuation coefficient, 〈K_d(λ)〉₁. As expected, the first attenuation depth z₁ is located much closer to the ocean surface for data from Group 1 (Fig. 12a) compared with Group 2 (Fig. 12b) and Group 3 (Fig. 12c), especially across the blue-to-green region of the spectrum. In the red part of the spectrum, where pure water absorption dominates the attenuation of E_d(λ), the differences between the three groups are small. It is also notable that the spectral behavior of z₁ for Group 3 (Fig. 12c) that represents relatively clear ocean waters is remarkably similar to the spectral shape of the pure water absorption coefficient.