1. Introduction

The Southern Ocean (SO, south of 50°S), occupying ∼10% of the world’s ocean surface area, is a typical high-nutrient low-chlorophyll (HNLC) region [1], yet an important carbon sink, accounting for ∼25% of the global ocean carbon dioxide uptake [2–4]. Currently, SO is undergoing unprecedented physical and chemical changes owing to global warming [5, 6]. Thus, there is a high demand for accurate assessments of primary production (PP) and carbon export in the region, as well as their long-term trends, so that a better understanding of the global carbon cycle and the response of SO to climate change can be possibly achieved. In this respect, satellite remote sensing is the only feasible means of evaluating PP and carbon export in the SO at large spatio-temporal scales.

Inherent optical properties (IOPs), including the spectral absorption and backscattering coefficients of different water components, determine light propagation within the water column and are closely related to the biogeochemical properties of the ocean. For instance, the absorption coefficient of phytoplankton (a_ph(λ), in m^-1) is particularly important for accurate estimation of PP [7–9], while the particulate backscattering coefficient (b_bp(λ), in m^-1) is a useful proxy for the concentration of particulate organic carbon (POC) in global oceans [10,11]. Colored dissolved organic matter (CDOM) contributes significantly to the dissolved organic carbon (DOC) pool, and absorption by CDOM (a_g(λ), in m^-1) is conventionally employed to remotely estimate DOC concentration, especially in estuarine and coastal waters [12,13]. Previous studies have investigated the spatial-temporal variabilities of PP and POC using remotely sensed IOPs in the SO [9,10]; however, the robustness of satellite IOPs products has not been well investigated. There have been extensive efforts towards a better understanding of the remotely sensed chlorophyll-a concentration (Chla) in SO [14–17], while evaluations of the inversion algorithms for IOPs in SO have scarcely been reported.

Over the last several decades, several semi-analytical algorithms (SAAs) have been developed to derive IOPs from remote sensing reflectance (R_rs(λ), in sr^-1) measured in situ or by satellite [18–21]. These SAAs can be categorized into two types according to their bottom-up and top-down strategies [22]. Top-down algorithms, such as the widely-used quasi-analytical algorithm (QAA) [23], first estimate the total absorption coefficient a(λ) and backscattering coefficient b_b(λ), and then partition a(λ) and b_b(λ) into IOPs by different water components. Bottom-up algorithms, on the other hand, derive component IOPs simultaneously with assumed IOPs spectral shapes via spectral optimization [19,21,24,25], where the Generalized IOP model (GIOP) is mostly employed by the international community [19]. However, these SAAs were mainly developed based on measurements collected from mid-low-latitude oceans, and their applicability in polar oceans has not been systematically evaluated.

In this study, we compiled a relatively large bio-optical dataset for SO from both field campaigns and published data for the evaluations of SAAs in SO. The main objectives of this effort include (a) evaluating the performance of QAA and GIOP using both field- and satellite-measured R_rs, and (b) quantifying and discussing the uncertainties associated with the IOPs inversions in the SO, particularly for the inversion of absorption coefficients by different water constituents.

2. Data and methods

2.1 in situ bio-optical dataset

We compiled a dataset with a wide dynamic range of optical properties in SO from two sources. The first one includes field measurements in the coastal waters of the Antarctic continent obtained during the 38th Chinese National Antarctic Research Expedition aboard the R/V XUELONG icebreaker, which is hereafter termed SO-38. The sampling locations of SO-38 are illustrated in Fig. 1 with red circles. The second one was screened from version 3 of the global bio-optical dataset (GBD-v3) by Valente et al. [26] for measurements collected in the Southern Ocean (south of 60°S). For simplicity, the subset of GBD-v3 in the Southern Ocean is termed SO-GBD hereafter, while the rest measurements in GBD-v3 are termed GBD-Global. Note that the bio-optical measurements in SO-GBD were mainly collected from the nearby waters of the Antarctic Peninsula (see the blue squares in Fig. 1).

 

Fig. 1. Sampling locations of the compiled bio-optical dataset in the Southern Ocean with red circles and blue squares standing for measurements from the SO-38 and SO-GBD datasets, respectively. Panel (b) shows the sampling locations in the nearby waters of the Antarctic Peninsula (AP, the black box in panel (a)). The green triangles in panel (a) represent stations with matched measurements from both in situ and satellite.

Download Full Size | PDF

2.1.1 SO-GBD dataset

The GBD-v3 dataset can be accessed directly from https://doi.pangaea.de/10.1594/PANGAEA.941318, which is a compilation of several publicly available datasets, such as the Marine Optical Buoy (MOBY), SeaWiFS (Sea-Viewing Wide Field-of-View Sensor) Bio-optical Archive and Storage System (SeaBASS), and NASA bio-Optical Marine Algorithm Data set (NOMAD), for measurements from global oceans between 1997 and 2021 [26]. We first extracted matched R_rs and IOPs data from their respective spreadsheets aggregated within +/- 6 nm of the VIIRS bands. Measurements in the Southern Ocean (south of 60°S) were then screened. As our main focus is on the descriptions of the methods used to determine the bio-optical properties in SO-38, readers are referred to the corresponding references for the protocols and practices of the measurements of different bio-optical properties in SO-GBD [26]. In this effort, we further screened the R_rs measurements in SO-GBD using a quality assurance (QA) system [27]. Specifically, the QA system assigns a score between 0 and 1 to each R_rs spectrum, with 1 indicating perfect quality and 0 for questionable spectrum. R_rs measurements with a QA score less than 0.5 were considered low quality and were removed from further analysis. As a result, the SO-GBD dataset consists of 88 valid measurements of R_rs(λ) and a_ph(λ) and 84 valid measurements of a_dg(λ), with these measurements mainly collected between 2000 and 2007.

2.1.2 SO-38 dataset

The 38th Chinese National Antarctic Research Expedition was conducted between December 2021 and March 2022. During this expedition, radiometric measurements were carried out and water samples were collected at 142 stations, with their locations shown in Fig. 1(a). Note that at some of these stations, only water samples were collected, and radiometric measurements were absent due to rough sea conditions or bad weather.

Spectral R_rs(λ), defined as the ratio of water-leaving radiance (L_w(λ), in W·m⁻²·nm⁻¹·sr⁻¹) to the downwelling plane irradiance just above the sea surface (E_d(λ, z = 0⁺), in W·m⁻²·nm⁻¹), were measured using the above-water approach [28]. At each station, radiance measurements of the upwelling radiance above the surface (L_u(λ)), standard Spectralon plaque (L_plaque(λ)), and downwelling sky radiance (L_sky(λ)) were acquired sequentially using a GER1500 spectroradiometer (Spectra Vista Corporation, USA). Note that GER1500 covers a spectral range of 350–1050 nm with a spectral resolution of 3 nm. The viewing geometry of each radiance measurement conformed to NASA protocols [29]. The field-measured R_rs(λ) were then calculated as

The absorption coefficients of suspended particles (a_p(λ)), non-algae particles (a_d(λ)), and phytoplankton (a_ph(λ)) were determined using the filter-pad technique (QFT) [32]. Specifically, surface water samples were collected and immediately filtered through a 25 mm diameter Whatman GF/F glass fiber filter (pore size of 0.7 µm) under low vacuum. The volume of the water sample was typically between 0.1 and 2 L to ensure that the absorptance of the particle-loaded filter did not exceed 0.4, as recommended by the IOCCG (International Ocean Color Coordination Group) Protocol [33]. The sample filters were stored at -80°C refrigerator until laboratory analysis. Briefly, the optical density (OD(λ)) of the sample filters between 250 and 800 nm was determined by a dual-beam PE Lambda 950 spectrophotometer equipped with an integrating sphere (150 mm in diameter) using the Transmittance-Reflectance (T-R) method [34,35], with the measured OD(λ) later converted to the spectral a_p(λ) by accounting for the path-length amplification factor (β). Stepwise procedures for obtaining a_p(λ), including the calculation of β, can be found in the latest IOCCG protocol [33]. a_d(λ) were calculated from absorptance measurements on the same sample filters after pigment extraction with methanol, whereas a_ph(λ) were simply calculated by subtracting a_d(λ) from a_p(λ).

 

Fig. 2. Spectra of field-measured (a) R_rs(λ), (b) a_dg(λ), and (c) a_ph(λ) from both SO-38 and SO-GBD datasets.

Download Full Size | PDF

Measurements of a_g(λ) were performed according to the IOCCG protocol [36]. Briefly, water samples were filtered through Millipore membranes (pore size of 0.2 µm) under a low vacuum immediately following collection. Sample filtrates were then collected in amber glass bottles and stored in the dark in a -20°C refrigerator during the cruise. The sample filtrates were unfrozen in a dark environment when returning to the laboratory, and the absorbance of the sample filtrates was measured using the same PE Lambda 950 spectrophotometer with a 10-cm cuvette, where Milli-Q water was used as the standard reference water blank.

Note that both a_p(λ) and a_g(λ) samples were prepared in duplicate during the cruise, and the final a_p(λ) and a_g(λ) were computed as an average of the duplicated samples after removing questionable spectra via visual inspection. In addition, a baseline correction was used to remove the residual scattering errors in the measured component absorption coefficients by subtracting the values of a_p(750), a_d(750), and a_g(750) from the respective a_p(λ), a_d(λ) and a_g(λ) spectra. Finally, the total absorption coefficient of seawater, a(λ), was calculated as the sum of a_p(λ), a_g(λ), and the absorption coefficient of pure seawater, a_w(λ). In this effort, a_w(λ) values were taken from Pope and Fry [37] for 550 - 750 nm, but employed values of Lee et al. [38] for 350–550 nm, as the latter is closer to the ‘true’ absorption coefficient of pure seawater in the 350–550 nm domain [39]. In this effort, a_w(λ) values are kept constant, although they vary slightly with temperature [40]. a_dg(λ) were calculated as the sum of a_d(λ) and a_g(λ). In total, 72 concurrent measurements of R_rs(λ) and the component absorption coefficients were obtained for SO-38.

Note that in situ measurements of b_bp(λ) were absent in both SO-GBD and SO-38; therefore, we focused mainly on the evaluation of derived absorption coefficients by different water components. In general, the in situ data used for algorithm evaluation cover a variety of coastal and oceanic water regimes in the Southern Ocean with a wide range of different absorption properties (see Fig. 2(b) and 2(c)). For all measurements in the SO-GBD and SO-38, termed SO-All herein, a_nw(443) (total absorption coefficient without pure seawater, in m^-1) ranged from 0.0084 to 0.31 m^-1, and the ratios of a_ph(443)/a_nw(443) and a_dg(443)/a_nw(443) varied in the range of 13.3% – 87.2% (51.8% ± 17.0%) and 12.8% – 86.7% (48.2% ± 17.0%), respectively, with the numbers in the brackets referring to the mean value and the standard deviation.

2.2 Satellite data

The standard Level-2 R_rs(λ) products of Visible Infrared Imaging Radiometer Suite (VIIRS), downloaded from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov), were employed in this study to perform the matchup analysis. VIIRS, onboard the Suomi National Polar-orbiting Partnership (SNPP) satellite, provides radiometric measurements, at a spatial resolution of 750 m, at six spectral bands in the visible domain centered at 410, 443, 486, 551, 671, and 745 nm, respectively. Note that the standard NASA Level-2 R_rs(λ) was atmospherically corrected following the concept of the ‘black-pixel’ assumption in the near-infrared (NIR) bands [41], but using an iterative bio-optical modeling approach to accurately estimate the water-leaving radiance in the NIR bands [42]. In addition, R_rs(λ) products of R2022.0 reprocessing were used in this effort.

Here, we computed the median value of all valid satellite-measured R_rs(λ) from a 5-by-5-pixel box covering the in situ site as the matched VIIRS R_rs(λ) measurements. Note that a strict comparison between satellite and in situ data usually requires a time difference of less than 3 h between the two measurements [43]. However, due to persistent cloud coverage and occasional sea ice appearance in our sampling locations, only a small number of matchups can be acquired under this strict criterion. Thus, the time difference was set to 24 hours in this effort to acquire matchups. Additionally, we used the QA score to assess the spectral quality of the VIIRS R_rs(λ), and only R_rs(λ) spectra with QA score ≥ 0.5 were retained. Note that radiometric measurements were not always concurrently collected with water sampling in SO-38, especially when some water samples were collected at night. In total, we obtained 24 matched measurements between VIIRS R_rs(λ) and field-measured R_rs(λ), and 35 matched measurements between VIIRS R_rs(λ) and field-measured component absorption coefficients.

2.3 Semi-analytical algorithms

Two widely used SAAs, GIOP and QAA, were employed to evaluate their performance in the Southern Ocean. The GIOP model employs a quadratic expression to describe the relationship between the remote sensing reflectance just below the sea surface (r_rs(λ)) and IOPs (i.e., a(λ) and total backscattering coefficient b_b(λ)) [44]:

In Eq. (2), a(λ) can be calculated as the sum of all component absorption coefficients [45,46], and the absorption by each component can be expressed as a function of the eigenvector (mass-specific absorption spectrum, a*(λ)) and eigenvalue (magnitude or concentration, M):

Similar to a(λ), b_b(λ) can be expressed as,

In this study, the default configuration of GIOP, referred to as GIOP for simplicity, was employed. Specifically, the eigenvector a*_dg(λ) is denoted as exp(-S_dg × (λ – 443)), where S_dg describes the spectral slope and is set to 0.018 nm^-1. The eigenvector a*_ph(λ) was estimated using the approach described by Bricaud et al. [48] with remotely sensed Chla, and further with a*_ph(443) forced as 0.055 m² mg^-1 [19]. The band-ratio algorithm of OC3V was employed to estimate Chla from R_rs(λ) [49]. The eigenvector b*_bp(λ) is expressed as (551/λ)^Sbp, where S_bp determines the steepness of the power-law variation in b*_bp(λ) and is determined by the blue-to-green ratio of R_rs(λ) [23]. The eigenvalues, i.e., M_dg, M_ph, and M_bp, were estimated from R_rs(λ) using spectral optimization with the Levenberg-Marquardt optimization. The component absorption and backscattering coefficients were subsequently derived from the eigenvalues and eigenvectors.

On the other hand, QAA algebraically solves IOPs from R_rs(λ) for optically deep water [23]. Here, the latest version (version 6) of the QAA, referred to as QAA herein for brevity, was used. Stepwise descriptions of QAA version 6 can be found at http://www.ioccg.org/groups/Software_OCA/QAA_v6_2014209.pdf. Briefly, QAA separates the inversion process into two consecutive parts. The first part involves the retrieval of a(λ) and b_b(λ) from R_rs(λ). QAA initiates IOPs inversion with the empirical estimation of a(λ) at a reference wavelength (i.e., a(λ₀), λ₀ stands for the reference wavelength), and then solves b_b(λ₀) algebraically from R_rs(λ₀) and a(λ₀). Then, b_bp(λ₀) is calculated as the difference between b_b(λ₀) and b_bw(λ₀), and the spectral b_bp(λ) can be obtained from b_bp(λ₀) and the empirically-estimated S_bp following the power-law function. With known b_bp(λ), b_b(λ) can be obtained and, subsequently, the a(λ). Note that the reference wavelength is set to 551 nm for relatively clear waters with R_rs(670) less than 0.0015 sr^-1 and set to 670 nm for scenarios with R_rs(670) greater than 0.0015 sr^-1. In our compiled SO dataset, all R_rs(670) values were less than 0.0015 sr^-1. In the second part of the QAA, a(λ) is partitioned into a_ph(λ) and a_dg(λ) algebraically after estimating the ratio of ζ and ξ, where ζ and ξ are the ratios of a_ph(410)/a_ph(443) and a_dg(410)/a_dg(443), respectively. Unlike GIOP treating S_dg as a constant, S_dg in QAA is a function of r_rs(443)/r_rs(551). A simplified flow chart illustrating the inversion of IOPs by both GIOP and QAA can be found in Fig. 3.

 

Fig. 3. Simplified flow chart illustrating GIOP and QAA for the inversion of IOPs.

Download Full Size | PDF

Note that inelastic scattering processes, such as Raman scattering, could result in elevated R_rs(λ), particularly for the long wavelengths in oligotrophic waters. Therefore, prior to the implementation of both QAA and GIOP, a simple correction of R_rs(λ) for Raman scattering effects was performed following Lee et al. [50].

2.4 Statistical metrics

The performances of both SAAs were quantified using statistical metrics calculated between derived IOPs and known values from in situ measurements, which included the median percentage difference (ε), Root Mean Square Difference (RMSD), Mean Absolute Percentage Difference (MAPD), and linear regression parameters, namely the slope and coefficient of determination (R²). Linear regressions were performed on the logarithmically transformed IOPs. The first three metrics are defined as follows,

3. Results and discussion

3.1 Evaluation of derived IOPs from field-measured R_rs(λ)

The performance of both QAA and GIOP was first assessed using SO-38 and SO-GBD at the four VIIRS bands in the blue-green domain, i.e., 410, 443, 486, and 551 nm. Here we only present the evaluation results of a(λ), a_ph(λ), and a_dg(λ) derived by both SAAs. Retrieved IOPs at 670 nm are not necessary in this effort, as a(670) is dominated by the absorption of pure seawater in our compiled dataset.

3.1.1 Evaluation of derived a(λ) by both SAAs

We first evaluated the derived a(λ) using the two SAAs, with scatterplots between the derived and measured a(λ) presented in Fig. 4. The detailed statistical metrics used to quantify the performance of the two SAAs are listed in Table 1. To ensure consistency, the same a_w(λ) values are employed in the inversion schemes as those used in the calculation of a(λ) from the lab-measured component absorption coefficients.

 

Fig. 4. Comparisons between in situ measured a(λ) and that derived by QAA and GIOP at 410, 443, 486, and 551 nm, respectively. The solid black line indicates the 1:1 relationship, and the dashed lines represent the 1:1 line ± 20%log₁₀($\bar{a}$(λ)), with $\bar{a}$(λ) as the mean a(λ) of SO-All.

Download Full Size | PDF

Table 1. Statistics of derived a(λ) by both QAA and GIOP at 410, 443, 486, and 551 nm for SO-GBD, SO-38, and SO-All datasets, respectively

View Table | View all tables in this article

As shown in Fig. 4, retrieved a(λ) by both GIOP and QAA show good agreement with the measured a(λ). Overall, the MAPD values of a(λ) from both algorithms ranged from ∼8% at 551 nm to ∼31% at 410 nm, with remotely sensed a(λ) slightly lower than that from the water samples (see Table 1 for details). This is equivalent to those reported in the literature [18,19,51], but better than those reported in coastal regions such as the South China Sea [52]. A larger deviation of a(λ) at the shorter wavelengths was also observed in the high-latitude Arctic region [51], and there are a few reasons behind this pattern. First, the statistically good performance at 551 nm is largely attributed to the fact that a_w(λ) is the dominant component of a(λ) at 551 nm for SO-All (a_w(551)/a(551) ∼ 0.83). For the blue bands, on the one hand, there is an impact associated with the empirically estimated S_bp. As both QAA and GIOP assume the same spectral shape of b_bp(λ), uncertainties in S_bp could be propagated to the estimated b_bp(λ), and such impacts would be amplified for b_bp(λ) at shorter wavelengths [53,54], which will subsequently affect the analytical inversion of a(λ) from R_rs(λ). For instance, S_bp for waters in SO is generally between 1 and 2 [55], if S_bp changes from 1.5 to 1.8 (increased by 20%), derived b_bp(410) will be elevated by 9.3%, in comparison to a 3.8% increasing of derived b_bp(486). This is particularly important as the empirically estimated S_bp from the band ratio of R_rs(λ) could have relatively large uncertainties [54,55]. On the other hand, field-measured a(λ) could have relatively larger measurement errors in blue bands than in green bands. For instance, the residual scattering errors of lab-measured a_p(λ) and a_g(λ) typically decrease towards longer wavelengths, and the baseline correction approach may not be able to fully correct these errors in the blue bands [56].

Based on the metrics shown in Table 1, for the respective dataset of SO-GBD and SO-38, QAA and GIOP have overall similar performance in retrieving a(λ), as indicated by the values of ε and MAPD. However, between the two algorithms, QAA performed slightly better in SO-GBD in terms of RMSD and R², while GIOP predicted overall more accurate a(λ) in SO-38. Across different datasets, it appears that the uncertainties of a(λ) by both SAAs are relatively smaller for SO-GBD than for SO-38. For example, larger uncertainties of a(410) by both SAAs were observed in SO-38, with greater ε values in SO-38 (ε = -30.2% and -37.4% for QAA and GIOP, respectively) than that in SO-GBD (ε = -21.5% and -20.8%). These results suggest that there could be slightly larger uncertainties in SO-38 than in SO-GBD, either in the measured R_rs(λ), a(λ), or both. As discussed in previous studies, the above-water approach does not provide a direct measurement of L_w(λ), and it is quite challenging to adequately remove the surface-reflected light when the sea surface is roughened by waves, especially in the Southern Ocean, where the reflected light consists of glints from sky, sun, and broken-ice [30,57–59]. Nevertheless, the MAPD and ε values shown in this study were within the range observed in the algorithm evaluation reported in the literature [19,52].

3.1.2 Evaluation of derived a_ph(λ) by both SAAs

Following the same evaluation strategy for derived a(λ), we further compared the derived a_ph(λ) by GIOP and QAA with the measured a_ph(λ), with results presented in Fig. 5. Detailed statistical metrics are given in Table 2.

 

Fig. 5. Same as Fig. 4, but for derived a_ph(λ).

Download Full Size | PDF

Table 2. Same as Table 1, but for derived a_ph(λ)

View Table | View all tables in this article

As shown in Fig. 5, derived a_ph(λ) by both SAAs agree well with measured values at the four wavelengths, where nearly all points are well distributed along the 1:1 line, suggesting that both GIOP and QAA provide satisfying a_ph(λ) retrievals for the compiled dataset. In general, GIOP and QAA have quite comparable performances in predicting a_ph(λ) in SO-All. Statistically, GIOP slightly outperforms QAA in terms of MAPD, where the MAPD of the GIOP-derived a_ph(λ) ranges between 35.8% and 40.9% for 410–486 nm, which is ∼10% smaller than that of QAA-derived a_ph(λ). However, QAA could be deemed to perform better in terms of ε and RMSD (see metrics in Table 2). These metrics are quite consistent with the evaluation results of both SAAs using measurements from the South China Sea and the IOCCG dataset, where MAPD of derived a_ph(λ) in the blue bands ranges from 25% to 60% [19,52]. Relatively larger uncertainties were observed at a_ph(551) by both SAAs with MAPD of 123.0% and 56.8% for retrievals from QAA and GIOP, respectively. This large deviation in derived a_ph(551) is mainly because a_ph(551) values are very small for these waters (a_ph(551) only accounts for ∼ 8% of a(551)). Thus, the retrieval of a_ph(551) is highly sensitive to errors in the derived a(551), which is particularly true for QAA-derived a_ph(λ), where a_ph(551) is solely derived from a(551). In particular, a_ph(551) derived by QAA could even be negative for waters with very low a_ph(551). For GIOP-derived a_ph(551), it is an extension of the estimated a_ph(443) based on the pre-assumed a_ph(λ) spectral shape; thus, it is not highly dependent on the value of a(551). For SO-All, QAA yielded a few negative a_ph(551) retrievals (N = 28), as also observed in Arctic and low-latitude waters [51]. Those negative retrievals were removed from further statistical analysis. Consequently, there were fewer scatter points from QAA-derived a_ph(551) in Fig. 5(d) compared to those from GIOP. However, negative a_ph(551) retrievals could be easily recovered by linking a_ph(551) with a_ph(443) if needed [19,21].

Between SO-GBD and SO-38, it appears that both SAAs present better agreement in a_ph(λ) for SO-GBD. Statistically, retrieved a_ph(λ) by both SAAs showed overall higher values of ε, MAPD, RMSD, and slope, and smaller R² in SO-38 compared to that in SO-GBD. Generally, ε is greater than 0 for both SAAs for SO-38, indicating that the retrievals from these algorithms are higher than in situ measurements for SO-38, but the opposite is true for SO-GBD.

3.1.3 Evaluation of derived a_dg(λ) by both SAAs

A comparison between the derived a_dg(λ) by both SAAs and in situ a_dg(λ) is presented in Fig. 6, with the statistical measures tabulated in Table 3. Compared to the generally quite consistent retrievals of a_ph(λ) and a(λ), derived a_dg(λ) from both SAAs showed large deviations versus in situ measurements for both SO-38 and SO-GBD. Specifically, at the lower end (a_dg(443) < ∼0.01 m^-1), a_dg(λ) from SAAs are systematically higher, but a_dg(λ) from SAAs are lower at the higher end (a_dg(443) > ∼0.03 m^-1). As a result, while in situ a_dg(443) shows a range of ∼0.003 to 0.2 m^-1, a_dg(443) from R_rs(λ) is in the range of ∼0.003 to 0.03 m^-1, i.e., a much narrower dynamic range from remote sensing. In particular, there are a few stations where a_dg(443) from R_rs(λ) is ∼10 times lower than in situ-measured a_dg(443). This result is quite puzzling as very good agreements are observed for a(λ) and a_ph(λ) (see Fig. 4 and Fig. 5), especially in situ a(λ) is simply a sum of a_w(λ), a_ph(λ), and a_dg(λ), and the dynamic ranges for both of a_ph(λ) and a_dg(λ) in SO-All are similar (i.e., a_ph(λ) is not dominating a(λ)). In the following, we try to provide a thorough explanation of this deviation.

 

Fig. 6. Same as Fig. 4, but for derived a_dg(λ).

Download Full Size | PDF

Table 3. Same as Table 1, but for derived a_dg(λ)

View Table | View all tables in this article

As illustrated clearly with the QAA scheme, while the retrieval of a(λ) is governed by the spectral shape of R_rs(λ), the separation of a_ph(λ) and a_dg(λ) from a(λ) is determined by a(λ) in the blue bands and the spectral shapes of a_ph(λ) and a_dg(λ), i.e., the ζ (a_ph(410)/a_ph(443)) and ξ (a_dg(410)/a_dg(443)) parameters. Generally, ζ is within a narrow range, so it does not significantly impact the derivation of a_ph(λ) and a_dg(λ). For instance, if the value of ζ increases by 10%, the QAA-derived a_dg(443) will be changed by less than ∼2%.On the other hand, ξ, equivalent to exp(S_dg(443-410)), could change from 1.39 (S_dg = 0.01 nm^-1) to 1.93 (S_dg = 0.02 nm^-1), and higher S_dg value will result in lower a_dg(443). Therefore, if in situ data are assumed to be accurate, to balance the higher retrieval at the lower end and lower retrieval at the higher end, as shown in Fig. 6, it requires a higher S_dg value at the lower end, while a lower S_dg value at the higher end, rather than a nearly fixed S_dg value (∼0.016 nm^-1 for QAA and 0.018 nm^-1 for GIOP) in both SAAs. Such a spatial variation of S_dg is not yet known in the SO and has not been reported in the literature for other regions. Furthermore, a 100% increase of S_dg from 0.01 to 0.02 nm^-1 will only result in a_dg(410) increasing by ∼39% and a_dg(551) declining by ∼66%. In short, errors in estimated ζ and ξ could not be the primary sources to explain the large deviations observed in derived a_dg(λ).

This is further supported by evaluations of both SAAs with GBD-Global, where Fig. 7 compares the derived a_dg(443) with the measured values. For this more inclusive dataset, we do not observe the “twisted” pattern as shown in Fig. 6. Statistically, the slope and R² values between the estimated and measured log(a_dg(443)) are close to 1.0, which are significantly better than those of the SO datasets. Moreover, as shown in Fig. 7, the points of the SO dataset generally overlap with those of the GBD-Global, indicating that the quality of the SO dataset is overall equivalent to that of the global dataset, or at least no significant regional bias. These comparisons further indicate that the nearly constant S_dg values used in both SAAs are in general acceptable for the retrieval of a_dg.

 

Fig. 7. Comparisons between measured a_dg(443) and that derived by (a) GIOP and (b) QAA for measurements from GBD-Global (tan circles), SO-38 (purple squares), and SO-GBD (green circles) datasets. The solid black line indicates the regression line for all data and the dashed line for the 1:1 relationship.

Download Full Size | PDF

The “twisted” pattern shown in Fig. 7 is more likely due to the inconsistency between measured R_rs(λ) and measured a_dg(λ) in the blue bands, as revealed between measured R_rs(λ) and modeled R_rs(λ) using measured IOPs (see Fig. 8). For example, for a station where the estimated a_dg(λ) is much lower than the measured a_dg(λ) (the higher end as shown in Fig. 6(b)), Fig. 8(a) shows the comparison between measured R_rs(λ) and modeled R_rs(λ), where the latter one was simulated used in situ a(λ) along with an optimized b_b(λ) for the best match in the green-red wavelengths following Eqs. (2) – (3). In this case, the R_rs(410)/R_rs(443) ratio of the modeled R_rs(λ) is ∼0.9, but this ratio from the measured R_rs(λ) is ∼1.4, i.e., measured R_rs(blue) is significantly higher than modeled R_rs(blue). Since blue bands play the major role in the retrieval of a_dg(λ) and R_rs(λ) is an inverse function of a(λ), this no closure in R_rs(blue) caused a_dg(λ) derived from the two SAAs much lower than measured a_dg(λ). Conversely, for the case at the lower end of a_dg(λ) as shown in Fig. 6(b), the simulated R_rs(λ) with known a(λ) is much higher than the measured R_rs(λ) in the blue bands (see Fig. 8(b)), i.e., retrieved a_dg(λ) is higher than measured a_dg(λ) due to this no closure. Possible sources of errors or uncertainties for the no closure in R_rs(blue) include insufficient corrections of the residual scattering for measured a_dg(λ) [33,36] and incomplete removal of surface-reflected skylight for the measurement of R_rs(λ) [58], where both could significantly impact a_dg(λ) or R_rs(λ) in the blue bands. Nevertheless, Fig. 8 suggests that optical closure between measured IOPs and R_rs(λ) can be reached for the 450–700 nm domain, which supports the overall good agreement observed between measured and retrieved a(λ) and a_ph(λ). Further, due to that a_dg(λ) is just one component of the total absorption coefficient, the relative impact of a_dg(λ) is different between the evaluation of a_dg(λ) and the evaluation of a(λ). For instance, the ratio of a_dg(443) to a(443) is ∼0.42 ± 0.16 (significantly lower for longer wavelengths) for SO-All, i.e., any deviations in a_dg(443) will be about halved for deviations in a(443), which is consistent with the observations shown in Fig. 4 where overall better agreements are found in the longer wavelengths.

 

Fig. 8. Examples of measured R_rs(λ) compared with simulated R_rs(λ) using measured IOPs. (a) for a station with measured a_dg(λ) in the higher end (the black star shown in Fig. 6(b)), while (b) for a station with measured a_dg(λ) in the lower end (also a black star shown in Fig. 6(b)). The solid black line represents measured R_rs(λ), while the red dashed line for simulated R_rs(λ).

Download Full Size | PDF

In addition, the retrieval of a_ph(λ) by SAAs is primarily determined by R_rs(λ) at blue-green bands (i.e., 443 and 551 nm), while a_dg(λ) is largely determined by R_rs(λ) at blue bands (i.e., 410 and 443 nm). Thus, failure or no closure at the blue bands has relatively small impacts on the retrieval of a_ph(λ). In other words, the performance of both SAAs in retrieving a_dg(λ) is largely dependent on the accuracy of measured R_rs(λ) at blue bands, which highlights the importance of accurate atmospheric correction for the global estimation of a_dg(λ). On the other hand, as long as ocean color satellites can provide good R_rs(λ) products at 443 and 551 nm, small errors in R_rs(410) may not significantly affect the retrieval of a_ph(λ) by SAAs.

3.2 Evaluation of VIIRS R_rs(λ) and IOPs products

Since the ultimate objective of the IOPs inverse algorithm is to obtain the spatial-temporal distribution of IOPs in the global ocean, it is of great interest to evaluate the derived IOPs by both SAAs from satellite measurements. In this effort, VIIRS matchups, with the locations of satellite-in situ matchups shown in Fig. 1, were employed for the evaluation.

3.2.1 Evaluation of VIIRS-measured R_rs(λ)

The robustness of satellite-derived IOPs is, of the first order, determined by the accuracy of the atmospherically corrected R_rs(λ) from satellite measurements. Thus, we first evaluated the quality of the standard VIIRS R_rs(λ) product for the Southern Ocean. Specifically, VIIRS R_rs(λ) at 410, 443, 486, and 551 nm were evaluated against matched field-measured R_rs(λ) from the SO-38 dataset (N = 24) (see Fig. 9). Note that only good-quality VIIRS-R_rs(λ) (QA score ≥ 0.5) were employed for the matchup analysis.

 

Fig. 9. Comparisons between VIIRS-measured and matched in situ-measured R_rs(λ) at (a) 410 nm, (b) 443 nm, (c) 486 nm, and (d) 551 nm, respectively.

Download Full Size | PDF

VIIRS-measured R_rs(λ) agree overall well with in-situ R_rs(λ) at the four VIIRS bands with MAPD values ranging between 19.8% and 23.9%, and there are no systematic deviations between VIIRS and in situ R_rs(λ) at these bands (ε in a range of -1.3 – -13.6%). Slightly larger uncertainties were observed at the two blue bands in terms of MAPD and R², which could be probably due to uncertainties associated with the atmospheric correction for blue wavelengths [60]. Note that the error statistics of VIIRS-measured R_rs(λ) in this effort are overall consistent with previous studies on the evaluation of VIIRS R_rs(λ) [61–63], and the cross-sensor validations [64–66]. Thus, it is feasible to use the standard VIIRS R_rs(λ) product screened by the QA score to generate confident IOPs products in the Southern Ocean.

3.2.2 Evaluation of VIIRS-derived IOPs

Both GIOP and QAA were applied to VIIRS R_rs(λ) and the derived IOPs at 410, 443, and 486 nm were then compared with the matched field-measured a(λ), a_dg(λ), and a_ph(λ) (see Fig. 10), with detailed statistical metrics provided in Table 4. In general, as found with field-measured R_rs(λ), both SAAs have quite comparable performance in the retrieval of absorption coefficients by different water components. Moreover, uncertainties in VIIRS-derived a(λ) and a_ph(λ) by both SAAs are overall smaller than that of derived a_dg(λ), which are also consistent with results obtained from filed measured R_rs(λ). Statistically, MAPD of derived a(λ) and a_ph(λ) by both SAAs are ∼30% and ∼40%, respectively, in comparison to MAPD of ∼60% for derived a_dg(λ). In particular, a_ph(λ) from VIIRS by both SAAs explained more than ∼70% of the variability (R² > 0.7), providing confidence to use VIIRS-derived a_ph(λ) for further applications, such as the estimation of chlorophyll concentration [48,67] or primary production [7,8]. For a_dg(λ), however, there is also a “twisted” pattern between the retrieved and measured values, which is reflected in the very small slope and R² values. As indicated in the above section, this pattern is a puzzle, as it was shown in both SO-38 and SO-GBD, which is probably due to uncertainties in a_dg(λ) measurements and the inversion algorithms (i.e., S_dg in GIOP, ζ and ξ in QAA).

 

Fig. 10. Evaluation of VIIRS-derived IOPs by GIOP and QAA using the matchups. The upper, middle, and lower panels show the results for derived a(λ), a_ph(λ), and a_dg(λ), respectively. The orange circles and green triangles represent the results of QAA and GIOP, respectively. The solid black line indicates the 1:1 relationship and the dashed lines represent the 1:1 line ± 20% log₁₀($\overline {\textrm{IOPs}} $), with $\overline {\textrm{IOPs}} $ as the mean value of matched in situ IOPs.

Download Full Size | PDF

Table 4. Statistics of derived a(λ), a_ph(λ), and a_dg(λ) from VIIRS matchups by both QAA and GIOP at 410, 443, and 486 nm, respectively

View Table | View all tables in this article

The statistical measures for derived a(λ) and a_ph(λ) from VIIRS-measured R_rs(λ) are overall worse than those from field-measured R_rs(λ), which could be expected because errors in the atmospherically-corrected R_rs(λ) will inevitably propagate to the derived IOPs, especially for the short blue bands [60]. In addition, the mismatch between satellite and field measurements due to the time and space gap will also introduce deviations between satellite and field measurements [68–70], especially when we employed a relatively relaxed time window (i.e., 24h) to acquire satellite-field matchups (see Section 2.2). Nevertheless, from the evaluations of both satellite and field measurements, we can conclude that QAA and GIOP can derive overall satisfying retrievals of a(λ) and a_ph(λ) in the Southern Ocean, with MAPD of derived a(λ) and a_ph(λ) being ∼20% and ∼40%, respectively. However, large uncertainties are observed in derived a_dg(λ) (MAPD > 60%), and it may require more high-quality field data and regional optimized inversion algorithms to resolve the puzzle associated with the “twisted” pattern about a_dg(λ).

4. Conclusion

In this effort, a relatively large bio-optical dataset for the Southern Ocean was compiled to evaluate the performance of two widely used semi-analytical algorithms (QAA and GIOP) in the inversion of IOPs. To the best of our knowledge, this is the first systematic evaluation of both SAAs for waters in SO. Results of this effort show that both QAA and GIOP can provide overall satisfying retrievals of a(λ) and a_ph(λ), but there is a “twisted” pattern between the retrieved and measured a_dg(λ). While it is always an ongoing effort to validate satellite products and improve remote sensing algorithms, results here advocate strong efforts to acquire high-quality field measurements for the evaluation, and improved inversion of a_dg(λ) in the Southern Ocean.