Study of the Profile Distribution of the Diffuse Attenuation Coefficient and Secchi Disk Depth in the Northwestern South China Sea

Abstract

The diffuse attenuation coefficient for the downwelling irradiance $K_d\left(z,\lambda \right)$ is a critical parameter in terms of the optical properties of the ocean. In the northwestern South China Sea, there are complex physical processes, and the accurate estimation of $K_d\left(z,\lambda \right)$ in the northwestern South China Sea is critical for the study and application of the underwater light field and water constituents. In this study, using Hydrolight 6.0 (HL60) software, $K_d\left(z,\lambda \right)$ was simulated based on the inherent optical properties (IOPs) and chlorophyll a concentration dataset in the northwestern South China Sea. The simulations were in good agreement with the results calculated by the model of Lee (2005), and the spectral characteristics of $K_d\left(z,\lambda \right)$ were consistent with several oceanic types according to Jerlov’s classification. The horizontal and vertical distribution characteristics of $K_d(z,490)$ were studied in the two typical upwelling areas of eastern Hainan Island and eastern Vietnam. $K_d(z,490)$ in eastern Hainan Island exhibited an overall decreasing trend from west to east at the same depth, while the vertical depth of the maximum value of $K_d(z,490)$ in eastern Hainan Island was found to increase from west to east, which was significantly associated with the distribution trend of the temperature and salinity. $K_d(z,490)$ in eastern Vietnam exhibited unique horizontal and vertical distribution characteristics due to upwelling, with a low temperature and high salinity. A satisfactory linear relationship between $K_d(z,490)$ and $K_d\left(z,\lambda \right)$ was found from $420 \mathrm{t}\mathrm{o} 580 \mathrm{n}\mathrm{m}$ with $R^2$ > 0.76, root mean square (RMSE) $\le$ 0.010 ${\mathrm{m}}^{-1}$, and mean absolute percentage error (MAPE) < 9%, and this result indicated that $K_d\left(z,\lambda \right)$ from $420~580 \mathrm{n}\mathrm{m}$ could be estimated with $K_d(z,490)$. The regression accuracy sharply decreased after 580 nm, indicating that $K_d\left(z,\lambda \right)$ estimation based on $K_d(z,490)$ can be more suitably achieved from 420~580 nm and becomes inaccurate after 580 nm. Based on the simulations, an empirical relationship for $Z_{SDD}$ estimation involving $K_{d,10\%}\left(490\right)$ was developed, and $Z_{SDD}$ in the northwestern South China Sea was calculated, with a range of 5–23 m and a suitable agreement with $Z_{SDD}$ obtained via the method of Lee (2018).

1. Introduction

The diffuse attenuation coefficient for the downwelling irradiance, $K_d\left(z,\lambda \right),$ as a function of the depth, $z$, and wavelength, $\lambda ,$ can be used to characterize the decrease in the ambient downwelling irradiance, $E_d\left(z,\lambda \right)$ [1]. This parameter can be defined as [2]:

$$K_d\left(z_1\leftrightarrow z_2,\lambda \right)=\frac{1}{z_2-z_1}\ln \left(\frac{E_d\left(z_1,\lambda \right)}{E_d\left(z_2,\lambda \right)}\right)$$

(1)

where $z_1$ and $z_2$ denote two profile depths that are very near to each other.

For $E_d\left(z,\lambda \right)/E_d(0^{+},\lambda )<1\%$, $K_d\left(z,\lambda \right)$ cannot be calculated with Equation (1) [3]. $E_d(0^{+},\lambda )$ denotes the downwelling irradiance just above the surface.

This quantity has been considered one of the key indicators to describe the attenuation of the downwelling solar radiation within the water column, as well as a critical factor for the Secchi disk depth (SDD) and euphotic zone depth [4]. All these parameters are important in studying the underwater photovoltaic (PV) power generation [5], heat budget of the oceanic mixed layer [6], photosynthesis of phytoplankton, and warning of harmful algal blooms [7,8].

The South China Sea (SCS) is the largest trophic marginal sea in the western Pacific Ocean, and it plays an important role in regulating the regional climate and carbon cycle with its extensive area and volume [9,10]. In the northwestern South China Sea, there are physical processes such as Qiongdong upwelling, eastern Vietnam upwelling, western Guangdong coastal currents, freshwater from the Pearl River estuary, and some subtropical bay influenced by human activity [11]. These processes lead to complex physical and biochemical changes in the northwestern South China Sea. The accurate estimation of $K_d\left(z,\lambda \right)$ and its distribution in the northwestern South China Sea is critical for the study and application of the underwater light field and water constituents [7].

At present, $K_d\left(z,\lambda \right)$ can be obtained via in situ measurements, empirical algorithms, semi-analytical algorithms, and analytical algorithms [3,12,13].

Regarding in situ measurements, a profile of the downwelling irradiance, $E_d\left(z,\lambda \right)$, is first obtained from instrumentation, and $K_d\left(z,\lambda \right)$ can then be calculated via Equation (1). In situ measurements are generally classified into two types: (1) surface radiation measurements based on optical buoy platforms and (2) profile radiation measurements involving radiometer profilers. Because irradiance measurement is very sensitive to radiometer self-shading, platform stability, and environmental factors such as clouds, aerosols, and solar radiation, it is difficult to control the data quality of in situ $E_d\left(z,\lambda \right)$ measurements [14,15,16].

In empirical algorithms, the diffuse attenuation coefficient is directly described as a function of the ratio of the water-leaving radiance or remote-sensing reflectance between several different wavelengths. In 1981, based on the water-leaving radiance retrieved from Coastal Zone Color Scanner (CZCS) data, Austin et al. [17] built an empirical relationship between $K_d\left(490\right)$ and the ratio between the blue and green water-leaving radiances. Subsequently, to provide a high accuracy and regional adaptability, the empirical relationship determined by Austin et al. was improved and optimized [18,19,20]. In 2020, based on time series optical buoy data, according to the favorable correlation of the $K_d\left(490\right)$ and $R_{rs}$ band ratios with $R_{rs}\left(650\right)$/$R_{rs}\left(510\right)$ and $R_{rs}\left(555\right)$/$R_{rs}\left(510\right)$, a dual-band ratio empirical algorithm for $K_d\left(490\right)$ was established by Zhang et al. [21], and this algorithm is superior to other empirical algorithms because the $R_{rs}$ selected by this algorithm can fully reflect water body information and can be adapted to changes in the water composition in the study sea area. All these empirical algorithms are easy in operation and beneficial to calculate the $K_d\left(\lambda \right)$ at a large spatial scale in the surface, but they are limited in accuracy and generalization and cannot explain the variation of $K_d\left(z,\lambda \right)$ and the underwater light field.

To overcome the disadvantages of empirical algorithms, based on the radiative transfer equation (RTE), several analytical and semi-analytical algorithms for $K_d\left(z,\lambda \right)$ estimation have been proposed. In 1976, Preisendorfer [22] established a comprehensive theoretical framework for ocean color and radiative transfer. Analytical algorithms were built based on the RTE, while the developed semi-analytical algorithms were based on the simplified RTE [23,24,25]. Based on the RTE and the increasing computing ability, a numerical simulation model, namely Hydrolight, was established by Mobley [26] in 1994. Hydrolight provides exact solutions to estimate the apparent optical properties (AOPs: diffuse attenuation coefficient $K_d$, remote sensing reflectance $R_{rs}$, etc.) from the inherent optical properties (IOPs: absorption coefficient $a$, attenuation coefficient $c$, volume scattering function $\beta$, etc.) [8], and this model is a powerful tool to study underwater light fields, facilitating progress in the study of the ocean primary productivity, ocean color remote sensing, environmental and ecological changes, and global energy cycles [2,27,28]. With the use of Hydrolight numerical simulations, Huang et al. [29] constructed a remote sensing reflectance dataset for Taihu Lake and compared it to measured remote sensing reflectance spectra, and based on the RTE and Hydrolight numerical simulations, a semi-analytical model for $K_d\left(z,\lambda \right)$ was developed by Lee et al. [2], which is called the model of Lee (2005) in this study, and this model can estimate $K_d\left(z,\lambda \right)$ based on the absorption and backscattering coefficients, which could be adapted well to both oceanic and coastal waters. Based on the simulations, Lee et al. [2] constructed a look-up table (LUT) for the values of the model parameters at specific solar altitude and profile depth ranges.

Previous studies of the diffuse attenuation coefficient, especially $K_d\left(490\right)$, have been conducted in the South China Sea more notably involving empirical algorithms and in situ downward irradiance data than involving analytical algorithms, and the original input parameters include the remote sensing reflectance, Chla, and in situ radiometric measurement data [30,31,32]. Uncertainties exist in both in situ and satellite-derived measurements of $K_d(z,\lambda )$. To address these issues, based on in situ measurements of the IOPs and boundary conditions such as the light field and meteorological parameters, the spectral and vertical distributions of $K_d\left(z,\lambda \right)$ in the northwestern South China Sea, especially the profile distribution characteristics in two typical upwelling regions, were studied via Hydrolight 6.0 (HL60) numerical simulations. Then, linear relationships between $K_d\left(z,\lambda \right)$ and $K_d\left(z,490\right)$ were derived by regressing the simulation data, and the Secchi disk depth in the northwestern South China Sea was estimated via $K_{d,10\%}\left(490\right)$ and $K_{PAR,10\%}$. This study provides strong support for the evaluation of the underwater light field from IOPs.

2. Methods and Data

2.1. Hydrolight

Hydrolight solves the RTE using a numerical technique to compute the radiance distribution in a natural hydrosol based on the IOPs of the hydrosol itself and appropriate boundary conditions at the surface and bottom of the water body, along with the radiance incident on the water surface [33].

In this study, Hydrolight version 6.0 (HL60) was used to simulate diffuse attenuation coefficients such as $K_d\left(z,\lambda \right)$ and $K_{PAR}\left(z\right)$, based on the water IOPs and surface and bottom boundary conditions. In HL60, seven options are available for IOP input: Constant IOPs, Classic Case 1 IOPs, New Case 1 IOPs, Measured IOPs, Case 2 IOPs, New Case 2 IOPs, and A user-defined IOP Model. The option of Measured IOPs, in which it is assumed that IOPs comprise pure water values plus “everything else” (particles and dissolved substances) measured by ac-s (the Wetlabs (now Seabird, in Bellevue, WA, USA) underwater hyperspectral absorption and attenuation meter) [34], HS6 (the HOBI Laboratories (in Bellevue, WA, USA) HydroScat-6 spectral backscattering sensor and fluorometer) [35], or similar instruments, was selected in this study, and Pope and Fry’s pure water absorption and seawater scattering coefficients were used [36] for the pure water component. The absorption and attenuation coefficients for “everything else” were measured via ac-s, and the backscattering coefficients for “everything else” were obtained in situ via HS6. Regarding the Measured IOPs model, the total minus water absorption coefficient ($a_{nw}(z,\lambda )$), total minus water attenuation coefficient ($c_{nw}(z,\lambda )$), particulate backscattering coefficient ($b_{bp}(z,\lambda )$), and chlorophyll a concentration (Chla$\left(z\right)$) are needed as inputs, a Fournier–Forand phase function with the given $b_{bp}(z,\lambda )/b_p(z,\lambda )$ ratio is used to estimate the scattering phase function of particulate matter in HL60, and the particulate scattering coefficient ($b_p(z,\lambda )$) can be obtained by subtracting $a_{nw}\left(z,\lambda \right)$ from $c_{nw}\left(z,\lambda \right)$. Vertical profiles of Chla($z$) are needed to compute chlorophyll fluorescence effects. In addition, the effects of inelastic scatter by colored dissolved organic matter (CDOM) fluorescence and Raman scattering are needed and could be calculated via the default function in HL60. The sky radiance distribution could be obtained from a semiempirical model, which was built in HL60. Moreover, we used the Harrison and Coombes normalized sky radiance model and RADTRANX sky irradiance model to specify the distribution of the sky radiance and partition the total irradiance into diffuse and direct irradiance components, respectively. The bottom boundary was set to an infinite depth. Then, with these settings, $K_d\left(z,\lambda \right)$ within 26 wavebands between 420 and 670 nm, separated by 10 nm, and $K_{PAR}\left(z\right)$ could be obtained in HL60. A simulation flowchart for HL60 is shown in Figure 1.

2.2. Cruise and Sampling

Two cruises were carried out in the northwestern South China Sea to collect IOPs: 9 August to 2 September 2013 (cruise 2013) and 30 June to 17 July 2015 (cruise 2015). Figure 2 shows the locations of the sampling stations. These cruise datasets contain match-ups of vertical profiles of $a_{nw}(z,\lambda )$, $c_{nw}(z,\lambda )$, $b_{bp}(z,\lambda )$, and Chla($z$) for 55 stations, along with the solar zenith angle, computed based on the time, latitude, and longitude of each station. These in situ measurements include a maximum depth in 99 m (offshore) and a minimum depth in 22 m (nearshore). Additionally, various meteorological parameters (e.g., wind speed and relative humidity) were collected simultaneously during sampling. The sampling area mainly involves two typical upwelling areas of eastern Hainan Island and eastern Vietnam with significant bio-optical variability and strong hydrodynamic processes.

2.2.1. Absorption and Attenuation Measurements

Vertical profiles of the total minus water absorption coefficient ($a_{nw}(z,\lambda )$) and total minus water attenuation coefficient ($c_{nw}(z,\lambda )$) were measured via ac-s with 82 wavebands between 401.6 and 744.1 nm and a path length of 25 cm. Milli-Q ultrapure water was used to calibrate the drift of the ac-s instrument and as a reference blank before the cruises. To mitigate any temperature and salinity effects on the ac-s measurements, a simultaneous seabird conductivity–temperature–depth (CTD) sensor (SBE 911plus CTD) was used to obtain temperature and salinity profiles [37]. In addition, all measurements were corrected for scattering by setting the absorption values as zero near the infrared wavelengths (716 nm was used in this study) [38]. As such, first, outlier data, including data less than or equal to zero, data deviating from the normal range, and same data values at the same depth, were excluded. Then, using the piecewise cubic interpolation (pchip) method, the data were interpolated to a profile resolution of 1 m, and finally, the interpolation data were smoothed via the Savitzky–Golay filtering method [39].

2.2.2. Backscattering Measurements

Vertical profiles of the particulate backscattering coefficient ($b_{bp}(z,\lambda )$) were obtained from the HS6 instrument. The HS6 instrument could be used to measure vertical profiles of the volume scattering function ($\beta$) at a single angle along the backward direction ($140°$) at six wavelengths (420, 442, 488, 532, 590, and 676 nm). The measurement data were processed based on the standard processing procedure provided by HOBI Laboratories, and the following sequential steps were employed: correction for the dark current, conversion of the raw signal into values of the volume scattering function $\beta (140,z)$, improvement in the accuracy of $\beta (140,z)$ via sigma correction based on the attenuation coefficient c, calculation of the particulate volume scattering function ${\beta }_p(140,z)$ from $\beta (140,z)$ under subtraction of the volume scattering function of seawater ${\beta }_w(140,z)$ [40], and estimation of the particulate backscattering coefficient, $b_{bp}(z,\lambda )$, from ${\beta }_p(140,z)$ multiplied by the relevant factor, ${\chi }_b$ [35]. Similarly, the same outlier removal, data interpolation, and smoothing processes were applied to the corrected data, similarly to the ac-s measurement data.

2.2.3. Chla Data

To correct for the effect of inelastic scatter of chlorophyll fluorescence, Chla($z$) can be adopted as an input parameter. In this study, Chla($z$) was calculated from $a_{nw}(z,\lambda )$ using the absorption line height method (denoted as the aLH method) and a profile resolution of 1 m in accordance with the ac-s and HS-6 data.

Scholars have shown that there is a power function-based correlation between the absorption line height at $676 \mathrm{nm}$ and Chla($z$) [41,42], and an absorption line height method (i.e., the aLH method) was developed for Chla($z$) estimation from the absorption line height, $a_{LH}\left(z,676\right)$, at $676 \mathrm{nm}$ [42]. The calculation process of the aLH method can be described as follows:

First, the absorption of the baseline at a reference wavelength, $a_{BL}\left({z,\lambda }_{ref}\right)$, of $a_{nw}(z,650)$ and $a_{nw}(z,715)$ can be calculated as:

$$a_{BL}\left({z,\lambda }_{ref}\right)=\frac{a_{nw}\left(z,715\right)-a_{nw}\left(z,650\right)}{715-650}\ast \left({\lambda }_{ref}-650\right)+a_{nw}\left(z,650\right)$$

(2)

Second, the absorption line height at 676 nm, $a_{LH}(z,676)$, can be calculated from the difference between $a_{nw}(z,676)$ and $a_{BL}(z,676)$ as:

$$a_{LH}\left(z,676\right)=a_{nw}\left(z,676\right)-a_{BL}\left(z,676\right)$$

(3)

Finally, Chla($z$) can be estimated from $a_{LH}\left(z,676\right)$ as follows:

$${Chla}_{LH}\left(z,676\right)=A\ast {a_{LH}\left(z,676\right)}^B$$

(4)

where A and B are 108.07 and 1.084, respectively, in the South China Sea [43], and ${Chla}_{LH}\left(z,676\right)$ is the chlorophyll a concentration derived via the aLH method.

2.2.4. Surface Boundary Conditions

The solar zenith angle was computed based on the time, latitude, and longitude of each station. Wind speed, relative humidity, sea level pressure, etc., were collected at the same time as the IOP dataset. Visibility data were collected from a Global Atmospheric Reanalysis (CRA-40) dataset released by the China Meteorological Administration (CMA) (https://data.cma.cn/, accessed on 20 January 2023). The total cloud cover and precipitable water were obtained from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) Reanalysis I dataset (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html, accessed on 20 January 2023), and these reanalysis datasets were obtained by using a climate model to initialize a wide variety of weather observations (ship data, plane measurements, station data, satellite observations, etc.).

2.3. Accuracy Evaluation Indicators

To assess the performance of the results, the following statistical indicators were used: the coefficient of determination ($R^2$), root mean square error (RMSE), and mean absolute percentage error (MAPE), which can be defined as follows:

$$R^2=1-\frac{{\sum }_{i=1}^n(y_i-\widehat{y_i}{)}^2}{{\sum }_{i=1}^n(y_i-\overline{y}{)}^2},$$

(5)

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n(y_i-{\widehat{y}}_i{)}^2}{n}},$$

(6)

$$MAPE={\sum }_{i=1}^n\left|\frac{y_i-{\widehat{y}}_i}{y_i}\right|\times \frac{100\%}{n},$$

(7)

where $y_i$ denotes the true value, ${\widehat{y}}_i$ denotes the estimated value, $\overline{y}$ is the average true value, and $n$ is the number of data pairs.

3. Results and Applications

3.1. Spectral Characteristics of $K_d\left(z,\lambda \right)$

Figure 3 shows the spectral distribution of $K_d\left(z,\lambda \right)$ simulated by HL60 in all 55 station profiles. As shown in Figure 3, the simulation results exhibited significant spectral and spatial distribution characteristics. The minimum value of $K_d\left(z,\lambda \right)$ generally appeared at approximately 490 nm, while the maximum value occurred at approximately 670 nm. According to Jerlov’s classification of water types (as shown in Figure 3b), the spectra of $K_d\left(z,\lambda \right)$ at these stations did not significantly differ from several oceanic types, especially type ΙB, type II, and type III [44]. There seem to be two different trends from 420–470 nm, one with an obvious downward trend and the other with a stabilized trend. The spectral characteristic with the obvious downward trend from 420–470 nm seems to occur in the nearshore stations. At all stations, $K_d\left(z,\lambda \right)$ tended to increase with the wavelength from 550 to 670 nm, with a sharp upward trend from 570 to 600 nm, while the upward trend started to decline from 600 to 670 nm.

3.2. Spatial and Vertical Distribution Characteristics of $K_d(z,\lambda )$

Figure 4a–c show the depth–wavelength distribution of $K_d\left(z,\lambda \right)$ on stations B2 in 2013, Q45 in 2013, and S63 in 2015, which can represent three water types: turbid water, relatively clear water, and clear water. Among all these water types, there were very obvious wavelength distribution characteristics of $K_d\left(z,\lambda \right)$ in the vertical profile. Among all of the profiles shown, the magnitude of $K_d\left(z,\lambda \right)$ uniformly increased with increasing depth. The $K_d\left(z,\lambda \right)$ profile at approximately 590 nm was an exception to this trend, as the magnitude decreased with increasing depth, and this feature could be observed in all of the wavelength bands > 590 nm. As shown in Figure 4a, the spectral variation in $K_d\left(z,\lambda \right)$ of the vertical profile was obvious from 420~670 nm; $K_d\left(z,\lambda \right)$ at approximately 490 nm was the minimum at the surface, the minimum $K_d\left(z,\lambda \right)$ value then moved toward higher wavelengths with increasing vertical depth, the minimum of $K_d\left(z,\lambda \right)$ moved to 560 nm, and the water window did not involve 490 nm but 560 nm. As shown in Figure 4b,c, there occurred relatively little change in the minimum of $K_d\left(z,\lambda \right)$ from 420~670 nm in the vertical profile, and $K_d\left(z,\lambda \right)$ at approximately 490 nm indicated the greatest ability for penetration in seawater, which makes $K_d\left(z,490\right)$ particularly unique in ocean color remote sensing, seawater transparency, and euphotic zone depth applications [33].

The spatial distribution of $K_d\left(z,490\right)$ at the different water depths, including 0, 5, 10, 15, and 20 m, at all stations (all stations were surveyed to depths of more than 20 m but less than 25 m) is shown in Figure 5. Figure 5, Figure 6 and Figure 7 were generated by Ocean Data View (ODV), and the Weighted-average gridding method was used to interpolate the station values. The contours of $K_d\left(z,490\right)$ in the five water layers revealed consistent spatial distribution characteristics. At all stations in the same vertical layer, $K_d\left(z,490\right)$ exhibited a decreasing trend from nearshore to offshore. In the five water layers at the same station,$K_d\left(z,490\right)$ indicated an increasing trend with increasing vertical depth, and the increase rate decreased with increasing offshore distance, i.e., the farther offshore, the gentler the change was. As shown in Figure 5, the contour lines significantly indicated that $K_d\left(z,490\right)$ exhibited a decreasing trend with increasing offshore distance. The smaller the offshore distance to the same water layer, the denser the contour lines were and the more significant the spatial variation rate of $K_d\left(z,490\right)$ was.

Four longitudinal sections (sections a, b, c, and d, as shown in Figure 2) in eastern Hainan Island were chosen to analyze the vertical and horizontal distribution characteristics of $K_d(z,490)$, as shown in Figure 6a–d. The spatial distribution characteristic of $K_d(z,490)$ exhibited an overall decreasing trend from west to east at the same depth in longitudinal sections a, b, c, and d. There was a vertical profile maximum of $K_d(z,490)$ in all sections. Moreover, the vertical depth of the maximum value of $K_d(z,490)$ in sections a and c (as shown in Figure 6a,c, respectively) was found to increase from west to east, and the difference in the vertical depth of the maximum value of $K_d(z,490)$ within the same section could exceed 20 m. Compared to that in sections a and c, there was a slight difference in the vertical maximum depth between the different stations in sections b and d, respectively, and the maximum depth reached approximately 50 m.

The vertical and horizontal distribution characteristics of $K_d(z,490)$ in the upwelling area of eastern Vietnam (e, f, and g sections in Figure 2) are shown in Figure 7a–c. In these sections, $K_d(z,490)$ revealed unique profile distribution characteristics in sections e, f, and g. The value of $K_d(z,490)$ in section e followed the same trend as that in sections a–d, with larger $K_d(z,490)$ values in the west and smaller $K_d(z,490)$ values in the east, and the maximum vertical depths of $K_d(z,490)$ were located at shallower depths in the west than in the east, with a depth difference of 20 m. The vertical depth of the maximum value of $K_d(z,490)$ in section f showed a convex structure from west to east. In contrast to section f, the vertical depth of the maximum value of $K_d(z,490)$ in section g showed a concave structure from west to east.

Overall, the vertical profile changes in $K_d(z,490)$ were more significant in the nearshore area. In all of the above sections from a to g, the maximum value of $K_d(z,490)$ reached approximately 0.14 ${\mathrm{m}}^{-1}$ at the Q18 station in section c, and the minimum value was 0.035 ${\mathrm{m}}^{-1}$ at the S49 station in section f.

3.3. Spectral Model of the Diffuse Attenuation Coefficient

Jerlov [44] proposed that $K_d(z,\lambda )$ could be expressed as a satisfactory linear function of $K_d\left({z,\lambda }_R\right)$, where ${\lambda }_R$ is the reference wavelength (e.g., 475 or 490 nm), which helped researchers derive $K_d(z,\lambda )$ at one or more wavelengths from $K_d\left({z,\lambda }_R\right)$. Subsequently, based on in situ measurements, with the use of 490 nm as ${\lambda }_R$, researchers have constructed several models similar to the relationship proposed by Jerlov for different sea areas [31,45]. In our study, based on the function established by Austin et al. (Equation (8)) [45], relationships between $K_d(z,490)$ and $K_d\left(z,\lambda \right)$ at $420~670 \mathrm{n}\mathrm{m}$ were established via HL60 $K_d\left(z,\lambda \right)$ simulations.

$$K_d\left(z,\lambda \right)=M\left(\lambda \right)K_d\left(z,490\right)+\mathrm{I}\left(\lambda \right),$$

(8)

As shown in Figure 8 and Figure 9, these relationships performed better between 420 and 580 nm, with $R^2>0.76$, RMSE$\le 0.010 {\mathrm{m}}^{-1}$, and MAPE $<9\%$. $R^2$ of these linear relationships, centered at 490 nm, tended to first increase and then decrease with increasing wavelength, while RMSE and MAPE tended to first decrease and then increase with increasing wavelength. The closer $\lambda$ was to 490 nm, the higher $R^2$ and the smaller RMSE and MAPE. The relationships performed the best at approximately 490 nm, with $R^2>0.9$, RMSE $\le 0.004 {\mathrm{m}}^{-1}$, and MAPE $<4\%$ at 450$~$540 nm. However, the linear regression accuracy decreased obviously from 590 nm to higher wavelengths, with the $R^2$ value decreasing from 0.578 at 590 nm to 0.128 at 670 nm, which indicated that the estimation of $K_d\left(z,\lambda \right)$ based on $K_d(z,490)$ is subject to a certain wavelength applicability range (420~580 nm). A poor quality of linear regression at higher wavelengths was previously reported by Austin et al. [45]. This result showed that to ensure the accuracy of $K_d\left(z,\lambda \right)$ estimation based on $K_d(z,490)$ in the northwestern South China Sea, the wavelength should be chosen from 420~580 nm.

We also found that coefficients $M\left(\lambda \right)$ and $I\left(\lambda \right)$ in Equation (8) from 420~580 nm could be described as quadratic functions of the wavelength via regression analysis. As shown in Figure 10, the quadratic regression accuracy of $M\left(\lambda \right)$ versus the wavelength indicated $R^2$ = 0.9995, RMSE = 0.009, and MAPE = 0.9%, and the regression of $I\left(\lambda \right)$ versus the wavelength achieved a higher accuracy with $R^2$ = 0.991, RMSE = 0.003${\mathrm{m}}^{-1}$, and MAPE = 18.9%.

3.4. Algorithm for $Z_{SDD}$ Estimation Based on the Diffuse Attenuation Coefficient

The Secchi disk depth (SDD) is an important parameter to describe the submarine optical field distribution and a key indicator in military–civilian integration fields, such as underwater target detection, communication, fishery production, and marine ecological environment monitoring [46,47,48]. Traditionally, SDD was measured by a white or black-and-white disk with a diameter generally about 30 cm [49]. A depth at which the Secchi disk is no longer visualized by an observer represents a quantitative measure of the transparency of ocean and lake waters. The traditional instrument is low cost and easy to operate [50], but it is limited in accuracy. Therefore, a series of regional algorithms or semi-analytical algorithms have been constructed in terms of the relationship between $Z_{SDD}$ and $K_d\left(\lambda \right)$ [4,49]. In contrast, the euphotic zone depth ($Z_{eu}$) is the depth where 1% of the photosynthetically available radiation (PAR($z$)) just beneath the water surface remains and is often quantified by the diffuse attenuation coefficient of the PAR($z$) ($K_{PAR}\left(z\right)$) [51].

$$Z_{eu}=\frac{4.605}{K_{PAR}\left(Z_{eu}\right)},$$

(9)

where $K_{PAR}\left(Z_{eu}\right)$ is the $K_{PAR}$ at the euphotic zone depth.

In the ocean color community, $Z_{eu}$ can be converted to $Z_{SDD}$ via Equation (10).

$$Z_{eu}\approx k{Z}_{SDD},$$

(10)

where $k$ denotes the proportionality factor.

As reported by Austin et al. [52], $K_{PAR}$ attained a notable linear correlation with $K_d\left(490\right)$, and $K_d\left(490\right)\le 0.25{\mathrm{m}}^{ -1}$.

$$K_{PAR}=M_{PAR}\ast K_d\left(490\right),$$

(11)

where $M_{PAR}$ is the slope of the linear regression equation.

In this study, in order to estimate $Z_{SDD}$, we used the average values of $K_{d,10\%}\left(490\right)$ and $K_{PAR,10\%}$, which were calculated between the surface and the depth where 10% of the $E_d\left(\mathrm{0,490}\right)$ and PAR(0) [4], to represent the $K_d\left(490\right)$ on Equation (11) and $K_{PAR}\left(Z_{eu}\right)$ and $K_{PAR}$ on Equations (9) and (11). Because $K_{d,10\%}\left(490\right)$ simulated by HL60 at all these stations was less than $0.16{\mathrm{m}}^{ -1}$, we performed a regression analysis of the simulation results for $K_{PAR,10\%}$ and $K_{d,10\%}\left(490\right)$ and obtained an $M_{PAR}$ value of 1.603 (Figure 11). $Z_{SDD}$ can be described as a function of $K_{d,10\%}\left(490\right)$ based on Equations (9)–(11) and $M_{PAR}$ in Equation (11), and can be written as:

$$Z_{SDD}=\frac{4.605}{k\ast M_{PAR}\ast K_{d,10\%}\left(490\right)}\approx \frac{0.821}{K_{d,10\%}\left(490\right)},$$

(12)

where $M_{PAR}$ is 1.603 and $k$ is 3.50 [4].

The $Z_{SDD}$ results obtained with Equation (12) ranged from 5 to 23 m, with an average of approximately 14 m. In particular, the values of $Z_{SDD}$ at three stations (B1, B2, and Q18) in 2013, which are nearest to the shore, were less than 10 m, and the value of $Z_{SDD}$ at station S29 in 2015 was greater than 20 m. In order to illustrate the spatial schematic diagram of $Z_{SDD}$ in the investigation area, ODV was also used the same as in Figure 5, Figure 6 and Figure 7, and the Weighted-average gridding method was used to interpolate the station values. The result is shown in Figure 12.

4. Discussion

4.1. Discussion about $K_d\left(z,\lambda \right)$

As demonstrated in Section 3.1 and Section 3.2, $K_d\left(z,\lambda \right)$ presented two significant spectral distributions in the northwestern South China Sea according to Jerlov’s classification; one is more in line with type III, and the other is more like type ΙB and type II. The reason for this difference is that the waters of the northwestern South China Sea are complex and belong to a mixed water where multiple types coexist, and the optical properties are more intricate [31]. As an example, the value of $K_d(z,490)$ can change from 0.03 ${\mathrm{m}}^{-1}$ to 0.14 ${\mathrm{m}}^{-1}$ with the investigated stations. Owing to diverse human and physical activities, the water at nearshore stations is more turbid than offshore stations, thereby causing a greater value and more significant variation of $K_d\left(z,\lambda \right)$ influenced by water constituents.

Regarding the vertical and horizontal distribution characteristics of $K_d(z,490)$ in eastern Hainan Island, these vertical and transverse distribution features of $K_d(z,490)$, which exhibited an overall decreasing trend from west to east, were found to be consistent with the distribution trends of the temperature and salinity in these longitudinal sections, namely, a significantly low temperature and high salinity were found along the western sea transects [53]. Wang et al. [53] had demonstrated the existence of upwelling on the east coast of Hainan Island during the cruise, so the upwelling process may be responsible for this vertical and horizontal distribution of $K_d(z,490)$. Additionally, in section c, the depths of the vertical profile maximum of $K_d(z,490)$ at stations Q23 and Q24 showed a more significant uplift relative to adjacent stations (as shown in Figure 6c), which is consistent with the uplift contours of the temperature and salinity in this section [53]. The reason for this trend may be that these stations are close to the Pearl River Mouth basin, the seabed topography is complicated and undulating, and there is an upwelling occurring in the shelf slope break zone [54]. Additionally, there was a difference at station Q6 in section d, and $K_d(z,490)$ significantly increased at a depth of 70 m in section d because of the small station depth (70 m occurs near the seabed), while significant scattering occurred due to the sediment on the seabed.

To investigate the reason for the vertical and horizontal distribution of $K_d(z,490)$ in eastern Vietnam, the distribution of temperature, salinity, and Chla in sections e, f, and g were compared with $K_d(z,490)$. Taking section f as an example, Figure 7 shows a comparative result similar to the region of eastern Hainan Island, where the distribution of $K_d(z,490)$ was significantly associated with the distribution trend of the temperature, salinity, and Chla. Considering that upwelling is often formed in the region of eastern Vietnam in summer [55] and the vertical distribution of nutrients and the bio-optical properties of water are greatly affected by water mixing [56], the internal reason for the observed convex distribution characteristic in section f might include the contribution of low-temperature and high-salinity upwelling to the promotion of the transport of nutrients from the bottom layer to the upper layer, triggering of phytoplankton growth, and increase in the Chla concentration.

4.2. Comparison of $K_d\left(z,490\right)$

$K_d\left(z,490\right)$, calculated by HL60 and the model of Lee (2005) [2] has been compared. In the model of Lee (2005), several model parameters of the solar zenith angles at $10°$, $30°$, and $60°$ and several profile depth ranges were determined to calculate $K_d\left(z,490\right)$. There are three stations (stations Q20 in 2013, Q47 in 2013, and D22 in 2015) with a solar zenith angle at around $10°$ ($\pm 2.5°$) and five stations (stations Q6 in 2013, Q34 in 2013, Q45 in 2013, Q48 in 2013, and D10 in 2015) with a solar zenith angle at around $30°$ ($\pm 2.5°$) that can be matched for comparison. $K_d\left(z,490\right)$ in 0–20 m profiles obtained from HL60 simulations and the model of Lee (2005) are shown in Figure 13. As shown in Figure 13, there was a suitable agreement between the results obtained from HL60 simulations and the model of Lee (2005). In all $10°$ stations, the $R^2$, RMSE, and MAPE were 0.891, 0.004 ${\mathrm{m}}^{-1}$, and 5.4%, respectively, and there was a better accuracy at all $30°$ stations with $R^2$ = 0.960, RMSE = 0.001 ${\mathrm{m}}^{-1}$, and MAPE = 2.1%. Among all these stations with a solar zenith angle at around $10°$ ($\pm 2.5°$) and $30°$ ($\pm 2.5°$), $K_d\left(z,490\right)$ calculated by the HL60 and the model of Lee (2005) had an excellent accuracy, with $R^2$ = 0.921, RMSE = 0.003 ${\mathrm{m}}^{-1}$, and MAPE = 3.3%.

Although an excellent agreement has been shown in Figure 13, there is also a subtle difference between the simulations and the model of Lee (2005), and the main factors contributing to this difference are: (1) meteorological parameters; (2) particle phase function; (3) error of Lee’s model itself. Firstly, the values of the model parameters were conducted by Lee et al. [2] based on Hydrolight simulations under clear sky conditions and a wind speed of 5 m/s, while the meteorological parameters input into HL60 in our study were not constant. For example, the wind speed input into HL60 for stations S66 and B1 were 11 $\mathrm{m} {\mathrm{s}}^{-1}$ and 2.2 $\mathrm{m} {\mathrm{s}}^{-1}$. The secondary reason for the difference was the particle phase function. The values of the model parameters in Lee (2005) were conducted following the averaged particle phase function, while the Fournier–Forand phase function was used in this study. The major factor for the discrepancy was that the model parameters conducted by Lee et al. [2] are functions of the profile depth, but in fact the constant was used in the driving process, and this meant that $K_d\left(z,490\right)$ has an error due to the constant parameters, just as Lee et al. [2] mentioned (the average error is $~3\%$ between modeled $K_d\left(z,\lambda \right)$ and Hydrolight determined $K_d\left(z,\lambda \right)$).

4.3. Comparison of $Z_{SDD}$

Lee et al. [4] established an empirical equation for seawater transparency calculation based on $K_{PAR}$ (Equation (13)), namely the method of Lee (2018) employed in our study.

$$Z_{SDD}\approx \frac{\alpha }{K_{PAR}},$$

(13)

where $\alpha$ is taken as a constant 1.48 in the study.

The $Z_{SDD}$ calculated by the method of Lee (2018) and the algorithm for $Z_{SDD}$ estimation based on $K_{d,10\%}\left(490\right)$ has been compared in this study. The average values of $K_{PAR,10\%}$ were calculated between the surface and the depth where 10% of the PAR(0) [4] was used to represent the $K_{PAR}$ on Equation (13). As shown in Figure 14, there was a suitable agreement between the results obtained with the two methods with slope = 0.972, $R^2$ = 0.804, RMSE = 1.326 $\mathrm{m}$, and MAPE = 8.9%. The uncertainty in the algorithm for $Z_{SDD}$ estimation based on $K_{d,10\%}\left(490\right)$ was mainly caused by the following: One is $K_{PAR}\left(Z_{eu}\right)$ in Equation (10); $K_{PAR}\left(Z_{eu}\right)$ in Equation (10) represents $K_{PAR}$ at the euphotic zone depth, while $K_{PAR,10\%}$ was calculated between the surface and the depth where 10% of PAR(0) was used in this study, and this approximation may generate some errors. The other uncertainty source is the proportionality factor of $k$ in Equation (11). Equation (11) is an empirical relationship between $Z_{SDD}$ and $Z_{eu}$, developed based on field measurements. The proportionality factor, $k$, is not a constant reported in a range of $~1–10$ [4,57,58]. The value of $k$ taken as a constant 3.5 in this study may introduce some errors to the estimation of $Z_{SDD}$ based on $K_{d,10\%}\left(490\right)$.

5. Conclusions

In this study, we simulated the diffuse attenuation coefficient for the downwelling irradiance, $K_d\left(z,\lambda \right)$, with HL60 based on measured vertical profiles of the IOPs from 9 August to 2 September 2013, and 30 June to 17 July 2015, in the northwestern South China Sea at 55 stations. The results between the simulations and the model of Lee (2005) had an excellent accuracy. The simulation results revealed significant spectral distribution characteristics, with a minimum at approximately 490 nm and a maximum at 670 nm. The spectral characteristics of the simulated $K_d\left(z,\lambda \right)$ values were similar to the spectra of several oceanic types reported by Jerlov [44].

The distribution characteristics of $K_d(z,490)$ in sections a–d in the eastern Hainan Island area showed an overall decreasing trend from west to east at the same depth, which was significantly associated with the distribution trend of the temperature and salinity. $K_d(z,490)$ in sections e–g in eastern Vietnam revealed unique horizontal and vertical distribution characteristics, which was completely consistent with the investigated vertical distribution characteristics of the temperature and salinity at the same time. During these cruises, $K_d(z,490)$ attained a global maximum of approximately 0.14 ${\mathrm{m}}^{-1}$ and a global minimum of 0.035 ${\mathrm{m}}^{-1}$.

The regression results indicated a satisfactory linear relationship between $K_d(z,490)$ and $K_d\left(z,\lambda \right)$ from $420\mathrm{t}\mathrm{o}580\mathrm{n}\mathrm{m}$ with $R^2$ > 0.76, RMSE $\le$ 0.010 ${\mathrm{m}}^{-1}$, and MAPE < 9%, and this result indicated that $K_d\left(z,\lambda \right)$ from $420~580\mathrm{n}\mathrm{m}$ could be estimated by $K_d(z,490)$, and the closer to 490 nm, the higher $R^2$ and the smaller RMSE and MAPE were between $K_d(z,490)$ and $K_d\left(z,\lambda \right)$. This conclusion is consistent with the results reported by Austin et al. [45]. It was also found that the regression slope,$M\left(\lambda \right)$, and intercept,$I\left(\lambda \right)$, of $K_d(z,490)$ and $K_d\left(z,\lambda \right)$ from 420~580 nm can be expressed as a quadratic function of the wavelength,$\lambda$, with $R^2$ values of 0.9995 and 0.991, respectively.

Based on the simulations, a linear relationship between $K_{d,10\%}\left(490\right)$ and $K_{PAR,10\%}$ for the northwestern South China Sea was first constructed, and an empirical relationship for $Z_{SDD}$ estimation based on $K_{d,10\%}\left(490\right)$ was then developed, and $Z_{SDD}$ at all 55 stations was calculated, with a range of 5–23 m and an average of approximately 14 m. $Z_{SDD}$ calculated in this study suitably agrees with the results obtained via the method of Lee (2018), with a slope of 0.972 and $R^2$, RMSE, and MAPE values of 0.804, 1.326 $\mathrm{m}$, and 8.9%, respectively.