1 INTRODUCTION

Coastal ecosystems are among the most productive ecosystems in the world, playing a major role in water, carbon, nitrogen, and phosphorous cycles between land and sea. Furthermore, coastal regions are home to about two thirds of the world's population (Cracknell, 1999). The suspended sediment concentration (SSC) and chlorophyll a (Chl) concentration are among the most important water quality components. They directly affect light distribution and, therefore, the ocean's primary production and growth of ocean agriculture (Gordon and McCluney, 1975). SSC data are essential for understanding water quality because suspended sediments can control the transport of nutrients, reduction of light penetration and fluxes of metals, radio-nuclides and organic micropollutants (Ouillon et al., 2004 ; Wang et al., 2009). Besides, degradation of water quality due to elevated SSC also affects aquatic organisms and has become a major global environmental problem (Schiebe et al., 1992).

The environment of East China Sea (ECS) has been impacted by huge stresses from anthropogenic activities and population growth in the Changjiang (Yangtze) River drainage basin and adjacent areas along the coasts. Improper use of natural resources and short-term economic objectives have resulted in severe environmental degradation in a fairly short time, and the degradation has now reached a level where the health and well being of the coastal populations are threatened (Li and Dag, 2004).

The concentrations of SSC and Chl in the water strongly influence the turbidity and color of water bodies (Harrington Jr et al., 1992 ; Pozdnyakov et al., 1998). In recent years, there is an increasing interest in satellite retrieval of water substances concentrations for optically complex Case 2 coastal waters. For this aim, numerous different algorithms based on surface reflectance were suggested. However, there still exists an insufficient understanding among the aquatic optics community as to the advantages and disadvantages of such algorithms, especially in the case of extremely turbid waters. Thus, it seems important to develop new SSC retrieval algorithms, which are adapted for turbid waters with a wide range SSC and Chl, to cover real situations in the Chinese inland, coastal, and estuarine waters.

Most of the present remote-sensing algorithms for SSC retrieval are based on using so-called spectral remote-sensing reflectance R_rs (λ) or some spectral combinations of this parameter (Ma and Dai, 2005 ; Zhou et al., 2006 ; Ouillon et al., 2008). This parameter is determined for the plane just above the air-water surface as the ratio of the spectral water-leaving radiance L_w (λ) to the spectral incident irradiance E_d (λ). Thus, an accurate determination of this parameter from the marine radiometric measurements is an important issue. However, there is still no consensus in the aquatic optics community about the standard method for R_rs (λ) determination. We will analyze different approaches to estimation of R_rs (λ) and a related parameter, the sky radiance spectral reflectance ρsky (λ).

In this paper three important issues necessary for satellite water quality remote sensing are considered: 1) an atmospheric correction method (Section 2.2); 2) radiometric models for R_rs (λ) and ρ_sky (λ)(Section 2.3); and 3) SSC versus R_rs (λ) models (section 2.4). We have solved these problems in a connection to the East China natural waters [mainly Changjiang River estuary (YRE) and its adjacent coastal area (ACA)], and generated several SSC maps for the region in 2008-2013.

2 MODELING APPROACH 2.1 Study area, ground and remote-sensing data

The East China inland, coastal and estuarine waters are characterized by suspended sediments with concentrations varied from ~10-100 g/m³ in the ECS waters, up to ~2 500 g/m³ in Changjiang River mouth and up to ~5 000 g/m³ in the Subei Bank of the Yellow Sea and Hangzhou Bay (Shen et al., 2010a, 2010b, 2013, 2014 ; He et al., 2012, 2013 ; Sokoletsky et al., 2014a).

Three databases were used in the present work for model development and validation purposes: 1) laboratory tank measurements performed in July 2006(118 samples), when the sediment samples were collected in Changjiang and Huanghe (Yellow) River estuaries (Shen et al., 2010a); and 2) two series of in situ measurements performed in May 2011(19 samples) and in February-March 2014(94 samples) in the estuary of the Changjiang River and coastal zone of the East China Sea. However, the simultaneous field measurements of Chl, SSC, and the spectral radiometric fluxes were made only during February 17 to March 10, 2014 scientific cruise.

Figure 1 shows the locations, including the Changjiang River and its estuary, Hangzhou Bay, Taihu Lake, Yellow Sea, and the ACA of the ECS. The wind speed W was measured by the Handheld anemometer and ranged from 0.5 to 8.0 m/s at average value (±standard deviation) of 2.8(±1.5) m/s. The solar zenith angle (SZA)θ₀ was estimated from the geographical coordinates of the sites and the times of measurements by the standard methods with the small modifications as described in Sokoletsky et al.(2011). The range of variability for θ₀ was 35 to 85 deg at average value (±standard deviation) of 54(±13) deg.

Methods of Chl and SSC measurements were described in Shen et al.(2010c) and Sokoletsky and Shen (2014), respectively. The Chl observations were performed only during February 17 to March 10, 2014(93 samples) with the average and standard deviation values for subsurface samples of 2.6(±1.2) mg/m³. The similar statistical measures for all 191 subsurface SSC samples are 243(±539) g/m³. Subsurface values of Chl and SSC during February 17 to March 10, 2014 varied in the ranges from 0.2 to 6.7 mg/m³ and 3.9 to 475 g/m³, respectively (Fig. 2). Values of SSC varied from 3.9-24(with the mean of 10 g/m³), 8.1-70(with the mean of 25 g/m³), and 61- 475(with average of 201 g/m³) g/m³ in moderate, turbid, and ultra-turbid waters, correspondingly; however, Chl values do not reveal any correlation with the water types (Fig. 2). A statistical analysis shows that the relation between Chl and SSC is very weak with the determination coefficient R² =0.017 and two-tailed significance level P =0.21.

It is important to note that these water types were determined based on the diffuse attenuation coefficient estimated at 550 nm, just below the sea surface, and for the sun at zenith, K_d (0°, 550 nm, 0-)(Sokoletsky and Shen, 2013, 2014). Exact definition of water turbidity classes accepted here is given in the Fig. 2 caption.

Figure 3 demonstrates color images for SSC and K_d (0°, 550 nm, 0-) derived from the in situ SSC measurements in the area. Both images reveal an ultra-turbid zone with SSC>150 g/m³ and K_d (0°, 550 nm, 0-)>4/m.

Three optical sensors of the Satlantic Inc. ^® Hyperspectral Surface Acquisition System (HyperSAS) were deployed for field measurements: 1) the vertical sky-viewing sensor designed for detection of the spectral downwelling irradiance E_d (λ); 2) the slant sky-viewing sensor, directed to 130° from the nadir, for the spectral incident radiance L_s (λ); and 3) the slant sea-viewing sensor, directed to 40° from the nadir, for the total spectral upwelling radiance, L_t (λ), which includes the spectral reflected radiance just above the air-water surface L_r (λ) and the spectral water-leaving radiance L_w (λ). These observational conditions are close to those advised by Mobley (1999). Auxiliary information, including local time, latitude, longitude, and cloud conditions, is recorded at each station. We used a Whatman G/F 0.45 m m filter to determine the SSC concentration.

The geocoded MODIS/Aqua L1B true color images collected during nine different days within the 2008-2013, also were used for the current study. The images were generated from the Level 1 Atmosphere Archive and Distribution System (LAADS) data. We used the MODIS data with the reduced-resolution (1 000 m) derived by the 4×4 pixel averaging of initial pixels. For converting these images to the SSC (false) images, atmospheric correction (AC) and SSC retrieval algorithms were developed and applied (Sections 2.2 and 2.4, respectively).

An AC procedure demands using several additional parameters, in addition to the spectral radiance received by the satellite's sensor. More specifically, the aerosol optical thickness at 550 nm, τ_a (865) and the Ångström parameter were obtained from the Goddard_Earth Sciences Data and Information Services Center (GES DISC); the atmospheric temperature and pressure (near to the air-sea surface), along with the total column ozone, were downloaded from the European Centre for Medium-Range Weather Forecasts (ECMWF); the relative humidity and precipitable water vapor have been derived from the Earth System Research Laboratory (ESRL).

2.2 Atmospheric correction

AC was implemented based on the analytic radiative transfer algorithm described in details by Sokoletsky et al.(2014a). This algorithm considers atmospheric transmittance separately for the direct (solar) and diffuse (sky) downwellling irradiances in a similar manner to that described by Gueymard (1995, 2001, 2003, 2008). The following atmospheric constituents and their optical properties were accounted: atmospheric molecules (Rayleigh scattering), aerosols (absorption and scattering), water vapor, mixed gases, ozone, and nitrogen dioxide absorption. The classical (a simple exponential) Bouguer-Lambert-Beer (BLB) law for absorption has been applied for the ozone and nitrogen dioxide transmittances, and the extended BLB law was used for the precipitable water vapor and mixed gases transmittances. However, for transmittances by atmospheric molecules and aerosols, the two stream radiative transfer (van de Hulst, 1980 ; Ben-David, 1995, 1997) model has been applied.

The spectral upwelling transmittances T_u (λ) were derived from the spectral downwelling plane transmittances T_d, _p(λ)(see definition of plane transmittance in Sokoletsky et al., 2014b) by applying the reciprocity principle (Sekera, 1970 ; Yang and Gordon, 1997). The upwelling radiance at the top of the atmosphere (TOA) level caused by the molecules (Rayleigh scattering) L_R, _TOA (λ) and aerosol particles L_a, _TOA (λ) was modeled separately, in accordance with the Gordon and Castaño (1987) and Rahman and Dedieu (1994) assumptions. The two-term Henyey- Greenstein (TTHG) scattering phase function (Gordon and Castaño, 1987):

where θ is a scattering angle and g is asymmetry parameter, has been used to present aerosol particles. The parameters α(λ), g₁, and g₂(λ) were set as functions of λ and τ_a (550) following the approximations suggested by Sturm and Zibordi (2002).

The spectral sun glitter radiance at the TOA level L_g, _TOA (λ) has been estimated by the air-water statistical model developed by Cox and Munk (1954) and extended by Gordon and Voss (2004). A condition of the sun glitter avoidance has been implemented as described by Remer et al.(2005).

Finally, the R_rs (λ) values were derived from the TOA spectral upwelling radiance measured by the satellite sensor L_u, _TOA (λ) by using equation:

where E_d (λ, 0+) is the modeled spectral global (from the sun and sky) atmospheric incident irradiance.

2.3 Deriving R_rs (λ) and ρ _sky (λ) from radiometric measurements

The first step in developing the retrieval SSC model is an analysis of in situ measurements, particularly, radiometric measurements. The spectral water-leaving radiance L_w (λ) may be expressed as

where L_t and L_r are the total above-water radiance and the surface-reflected incident (sun+sky) radiance, respectively. Here we omit an impact of white caps on the L_t because this impact is generally insignificant compare to the total radiance; this is especially true for the case of turbid waters with a substantial inorganic content. Thus, a L_w (λ) can be estimated if an accurate estimate of L_r (λ) can be obtained.

Now we rewrite an Eq.4 in the form:

where ρ(λ) is the spectral radiance reflectance that is affected by the geometrical conditions for sun and sensor, wavelength, wind speed, the sky radiance distribution, and field of view (Mobley, 1999 ; Lee et al., 2010). The ρ(λ) appears to be a highly variable parameter, especially as a function of the illumination conditions (Doxaran et al., 2004 ; Aas, 2010) and wavelength (Lee et al., 2010 ; Sokoletsky and Shen, 2014).

Below we consider six different models for deriving R_rs (λ)= L_w (λ)/ E_d (λ) from radiometric quantities.

2.3.1 Mobley model (MM)

Mobley (1999) numerically calculated the ρ values for various geometrical and cloud conditions and wind speeds W. His results may be described by the following equation:

where ρ_sky (λ) is the sky radiance spectral reflectance. The ρ_sky (λ) is assumed to be 0.028 for the heavily overcast sky at all values of λ, W and clear sky conditions at the wind speed W < 5 m/s, viewing zenith angle (VZA)θ_v =40°, azimuthal viewing direction of ϕ_v =135°, and SZA θ₀ >10°(Fig. 9, Table 1 and S ection 7 in Mobley, 1999). In most cases of measurements at the YRE and its ACA these conditions were nearly satisfied; thus, we used a single value of ρ_sky =0.028 for this method. A selected observational geometry yields the sun glitter angle θ_g >47.5° that allows to neglect a sun glitter (ρ _sun) for the MM and all other models (Remer et al., 2005).

2.3.2 Whitlock model (WM)

In the earlier publications, like Whitlock et al.(1981), ρ_sky has been often accepted as 0.02 for any conditions. We used this value for estimation of R_rs (λ) by Eq.6 in the WM model.

2.3.3 Ruddick-Mobley model (RMM)

This model based on Eq.6 with the ρ_sky calculated as a function of W and taking the ratio of L_s (750)/ E_d (750) into account according to Ruddick et al.(2006) :

2.3.4 Toole model (TM)

For this model (Toole et al., 2000) the Fresnel's reflectance ρ _Fr and two empirical (correction) parameters, c₁ and c₂ were introduced into Eq.6 as follows:

The coefficients c₁ and c₂ are given in the Table 7 by Toole et al.(2000) for the wind speeds W ≤5 m/s and W ＞5 m/s, and for the clear, broken (foggy or cloudy), and diffuse (overcast) skies. The refractive index of the water for calculation of ρ _Fr was accepted as 1.336.

2.3.5 Simis-Olsson model (SOM)

This model R_rs (λ) is estimated by Eq.6 with the ρ_sky calculated as in Sokoletsky and Shen (2014) :

This equation has been derived from 3 000 observations carried out by Simis and Olsson (2013) in the Baltic Sea.

2.3.6 Sokoletsky-Aas model (SAM)

Here we assume that an optical instrument pointing to the sky measures not only sky (diffuse) light, but also sun (direct) light. Then, introducing a parameter d_E (“irradiance diffuseness”) as the ratio of sky (diffuse) incident irradiance E_d, sky to the total (sky+sun) incident irradiance E_d (Woźniak, 1997), we get from Eqs.4 and 6

Here we used the expression for ρ_sky (λ) as a ratio of the spectral reflected sky radiance just above the airwater surface L_r, sky (λ) to the spectral incident sky radiance L_s, sky (λ), and assuming that L_s, sky (λ)= d_E (λ) L_s (λ). For this model, we expressed the spectral ratio L_r, sky (λθ₀ W)/E_d, sky (λ) as follows (Aas, 2010 ; Sokoletsky and Shen, 2014):

where coefficients A₀ - A₂ and B₀ - B₂ are given in Aas (2010) as functions of wavelength λ(at 405, 450, 520, 550, and 650 nm) and wind speed W (at 0, 5, and 10 m/s) at solar zenith angles θ₀ ranging from 37° to 76°. We performed calculations for any values of λ and W by linear interpolations and extrapolations, but taking θ₀ =37°, if θ₀ ≤37°, and θ₀ =76°, if θ₀ ≥76°.

We have used the following analytical approximation for the spectral parameter d_E (λ) derived from the atmospheric radiative transfer modeling (Sokoletsky et al., 2014a, c):

where θ₀ in degrees, λ in nm, and all coefficients are given in the Table 1.

The spectral plots of ρ_sky (λ)(Fig. 4) and ρ_sky (λ) d_E (λ)(Fig. 5) demonstrate a strong increasing of these parameters with λ for clear skies, the more temporal increasing for broken (cloudy or foggy) skies, and almost neutral behavior in the case of overcast skies. These facts differentiate the SAM model from the others, for which no dependency of ρ_sky versus λ is assumed. Note also an agreement of this conclusion with findings of Lee et al.(2010) and Sokoletsky and Shen (2014).

2.3.7 Comparison of different R_rs (λ) radiometric models

It is remarkable that models MM, WM, RMM, and SOM have a similar design as Eq.6. The difference between these four models is caused only by different expressions for the sky radiance spectral reflectance ρ_sky (λ). However, the TM model adds two empirical coefficients to better take into account the wind and sky conditions. The SAM model is similar on its structure to the TM model, if to assume that the parameter c₁(from the TM model) is similar to the parameter d_E (from the SAM model) and c₂ =0. However, the d_E is much more sensitive to the diffuseness of incident light than the c₁ in the TM model. Mathematical analysis of Eqs.6, 8, and 10 shows that an Eq.10(i.e., the SAM model) could yield bigger values for R_rs (λ) than Eqs.6 and 8.

This mathematical analysis is confirmed by numerical calculations showing generally the biggest spectral values of R_rs (λ) for the SAM model and the smallest values for the TM model. Moreover, the TM model sometimes shows negative values in visible, ultraviolet, and NIR spectral ranges, while the other models show either only positive values (WM) or negative values only in the NIR spectral range (MM, RMM, SOM, and SAM). For the whole spectral range of 350-900 nm, negative R_rs (λ) values were obtained in 0%, 0.7, 1.1%, 2.6%, 3.1%, and 25% cases for WM, RMM, SAM, SOM, MM, and TM models, respectively. We explain the significantly large percentage of negative R_rs (λ) values in the case of TM model by the large difference between environmental conditions in which the TM model was parameterized and conditions in which our measurements were performed. For example, the SSC values varied from 0.0 to 3.4 g/m³ and from 3.9 to 475 g/m³ in Toole et al.(2000) site of interest (the Santa Barbara Channel) and in our geographical area, respectively. Therefore, the SAM model based on the strong radiative transfer considerations and which yields generally positive R_rs (λ) values was selected as the main model for the further modeling.

Below we compare the six models described above for the R_rs (λ) retrieval for only two observation dates under clear (Fig. 6a) and cloudy (Fig. 6b) sky conditions. The differences between the SAM and all other models are generally smaller in clear sky conditions. To estimate numerically an impact of sky conditions on the differences between five models (the TM model has been excluded here from consideration), all samples were divided by three groups, namely, collected during clear (21 samples), cloudy or foggy (20 samples), and overcast (13 samples) skies. The results were expressed via the spectral variation coefficient Var (λ) for R_rs (λ)(Fig. 7) equal to the standard deviation of the R_rs (λ) normalized by the mean R_rs (λ) value calculated for each wavelength. Results confirm an increase in models' deviations with increasing atmospheric diffuseness. Notice also a general increasing of Var (λ) in the NIR comparative to the visible spectral range.

All symbols referred to the samples collected during February 17 to March 10, 2014, in the East China coastal area. Note that this figure is a type of a scatter plot with the omitted x-axis labels. This type has been chosen intentionally to stress a missing relationship between the values of two coordinates.

Even though numerical differences between R_rs (λ) estimated by different models may be large enough, using spectral ratios significantly decreases such differences (Sokoletsky and Gallegos, 2010 ; Wang et al., 2010b). Below (Fig. 8) we test six abovementioned models for the spectral ratio of R_rs (645.5)/ R_rs (553.6) at different sky conditions. The 645.5 and 553.6 nm are the central wavelengths of the first and fourth MODIS bands, respectively (Neteler, 2014), and their ratio was found to be usable for the remote sensing of water quality parameters in the YRE and ACA (Sokoletsky et al., 2014a). The Fig. 8 clearly shows that using the spectral ratio decreases the difference between the models, but this difference depends strongly on the spectral ratio value. More specifically, the models generally underestimate SAM values when this ratio approximately<1, and overestimate SAM values at R_rs (645.5)/ R_rs (553.6)>1, independently of sky conditions. Overall, the maximal discrepancies were found between the TM and all other models. Excluding again the TM model from consideration, we obtained the mean variation coefficient between the five models of 4.6(±7.5)%, 11(±15)%, and 16(±12)% for clear, broken (cloudy or foggy), and overcast sky conditions, respectively. Thus, accounting that remote sensing imagery is performed effectively under clear skies, using the R_rs (645.5)/ R_rs (553.6) seems to be safe enough for exploring in the YRE and its ACA area of China.

2.4 Modeling SSC versus spectral remote-sensing reflectance 2.4.1 Physical fundament for the R_rs (λ) versus SSC relationships

The R_rs (λ) versus SSC relationships follow from 1) relationships between the IOPs (spectral absorption a (λ), scattering b (λ), and backscattering b_b (λ) coefficients) and SSC; 2) relationships between the spectral subsurface remote-sensing reflectance r_rs (λ), IOPs and measurement geometry conditions; and 3) relationships between the r_rs (λ) and R_rs (λ). Below we consider these issues in some more details.

2.4.2 Relationships between the IOPs and SSC

A mathematical modeling has shown that the a (λ) dependency characterized by the spectrally non-flat curve with a shift of the spectral absorption minimum from the blue to NIR at increasing SSC (Fig. 9a) while the both b (λ)(Fig. 9b) and b_b (λ)(Fig. 9c) dependencies are spectrally flat with the strong decrease toward the NIR spectral range at any SSC. Mathematically, these dependencies may be described as follows:

where all optical properties are expressed in inverse meters, SSC in g/m³, and λ in nm; the subscripts w, p, and d refer to water, particles, and dissolved matter, respectively. The a_w (λ) values were interpolated from the Table V-1 by Lu (2006) at λ<600 nm, from Table 3 by Pope and Fry (1997) at 600≤λ<736 nm, from Table 1 by Buiteveld et al.(1994) at 736≤λ<801 nm, and from the webpage by Prahl and Jacques (1998) [based on the refractive indices by Kou et al.(1993) ] at λ≥801 nm. An expression for the b_w (λ) has been derived from the optical model by Zhang et al.(2009) at the temperature of 20℃ and salinity of 30 g/kg, λ=300 to 1 000 nm (with the mean absolute percentage error of 1.4% and R² of 0.998 2). The b_bw (λ) has been assumed to be exactly a half of the b_w (λ).

The Eqs.13-20 have been used for estimation of the spectral Gordon's parameter G (λ)= b_b (λ)/ [ a (λ)+ b_b (λ)](Fig. 9d), which may be considered as a spectral proxy for the R_rs (λ) in the first approximation (Haltrin and Gallegos, 2003 ; Sokoletsky et al., 2012 ; Sokoletsky and Fang, 2014).

We calculated this parameter (Table 2) at two MODIS wavelengths: 553.6 and 645.5 nm as a function of measured SSC for the May 2011 and February-March 2014 datasets. The accuracy of the model for G (553.6), G (645.5), and G (645.5)/ G (553.6) has been determined by the following statistical measures: determination coefficient R², normalized root-mean-square error NRMSE, and the linear percentage error LPE. The two last measures were calculated as follows (Lee et al., 2002 ; Sokoletsky et al., 2012):

where y _i and ${\tilde{y}}$_i are observed and predicted IOPs for the i th sample, respectively; n is the sample size; y is the observed mean value; and RMSE_log10 is the rootmean- square error in log ₁₀ scale. Table 2 contains these results along with the standard deviation, SD (y_i) and variability, Var (y_i)=[SD (y_i)/ y]×100% for G (553.6), G (645.5), and G (645.5)/ G (553.6).The table shows that using the spectral ratio G (645.5)/ G (553.6)[and, hence, the R_rs (645.5)/ R_rs (553.6)] is generally more numerically accurate than using the single wavelengths G (λ) values. A similar conclusion can be drawn for the correlation between predicted and measured G (λ).

2.4.3 Relationships between the r_rs (λ), IOPs and measurement geometry conditions

The subsurface remote-sensing reflectance r_rs (λ)= L_u (λ, 0-)/ E_d (λ, 0-) may be estimated from IOPs and measurement geometry conditions, for instance, by any of the five accurate approximations suggested and tested in Gordon et al.(1988), Lee et al.(1999), Albert and Gege (2006), Sokoletsky et al.(2009, 2012), and Sokoletsky and Shen (2014). We have used the combined Sokoletsky et al.(2012) and Albert and Gege (2006) approximation for r_rs (λ)(Albert and Gege, 2006 ; Sokoletsky et al., 2012 ; Sokoletsky and Shen, 2014):

where μ₁ and μ₂ are cosines of SZA and VZA, respectively, just below the air-sea surface.

2.4.4 Relationships between the r_rs (λ) and R_rs (λ)

The model R_rs (λ) values were derived from the r_rs (λ) by the following approximation:

where α(λ)= 0.3638 + 8.776 × 10^-4 λ - 9.193 × 10 ^{- 7} λ ² + 3.174 × 10 ^{- 10} λ ³ and β(λ)=1.357+8.608 × 10^-4 λ-6.347 × 10 - ⁷ λ ²(λ in nm). Equation 23 has been derived from the radiative transfer Aas-Højerslev model (Højerslev, 2001 ; Aas, 2010 ; Sokoletsky and Shen, 2014) at the solar zenith angle θ₀ =40°, W =5 m/s, and wavelengths of 400-800 nm with R² =0.999 5 and R² =0.990 3 for α(λ) and β(λ), respectively.

2.4.5 Relationships between the R_rs (λ) and SSC

Figure 10 shows the model dependencies of R_rs (λ, SSC) at θ₀ =θ v= 40° and wide ranges of λ and SSC computed by three steps in the manner described above, namely: 1) SSC→IOPs (Subsection 2.4.2); 2) IOPs→ r_rs (λ)(Subsection 2.4.3); and 3) r_rs (λ)→ R_rs (λ)(Subsection 2.4.4).

An important feature of the R_rs (λ, SSC) spectral behavior is an almost exact coincidence of a (λ, SSC) spectral minima (Fig. 9a) with the R_rs (λ, SSC) spectral maxima (Fig. 10a). Another important feature is an almost linear character of R_rs (λ) vs. SSC at small SSC (Fig. 10b). However, a relationship between SSC and R_rs (λ) becomes flatter with the SSC increasing. Fig. 10b demonstrates a shift of the SSC boundary between linear and nonlinear relationships with the increasing wavelength. For example, at wavelengths of 412-532, 595-676, 715, and 865 nm the boundary is achieved at SSC about 6, 10, 20, and 600 g/m³, respectively. These findings are in a qualitative accordance with the previous studies (Harrington Jr et al., 1992 ; Han and Rundquist, 1994 ; Wang and Lu, 2010 ; Shen et al., 2010a ; Qu, 2014).

One can use the angles θ₀ and θ_v different from 40°, which may lead to small differences from the results shown in Fig. 10. An analysis of errors in R_rs (λ θ₀, θ_v) at λ=550 nm caused by different θ₀ and θ_v is given in Annex. The relative errors in R_rs (λ=550 nm, θ₀, θ_v) vs. R_rs (λ=550 nm, θ₀ =40°, θ_v =40°) derived from Eqs.A1 and A3 are shown in Fig. 11. Numerically, they do not exceed 5% for the observational conditions during data collection.

2.4.6 Comparison between different R_rs (λ) versus SSC relationships

A large variability of different combinations of spectral reflectance may be used for SSC retrievals (Ma and Dai, 2005 ; Zhou et al., 2006 ; Ouillon et al., 2008), that makes the choice of the spectral SSC model a non-trivial task. Such models are inherently different for specific geographical sites and seasons, and even for different hours of the day.

Below (Table 3) we compare numerous SSC retrieval models developed for the turbid Case 2 waters based on different indices X. Testing these models was realized for the laboratory tank and in situ data (see Section 2.1) with the wavelengths taken from either the middle of the spectral ranges, or from the central wavelengths indicated by the authors. The basic R_rs (λ) values were derived from radiometric measurements, θ₀, τ a (λ) and W by the SAM method (Subsection 2.3.6). Table 3 lists literature models with the corresponding values of R², NRMSE, and LPE, calculated similarly to those measures applied to the IOPs (Eq.21). In case of necessity for modeling, the spectral water-leaving reflectance ρ_w (λ) has been replaced by R_rs (λ) according to the simple equation: ρ_w (λ)=π R_rs (λ)(Jerlov, 1976, p. 79; Gordon and Voss, 2004 ; Ruddick et al., 2006 ; Sokoletsky and Shen, 2014).

The Nechad et al.(2010) model was exploited as a function of the spectral ratios formed for three central MODIS wavelengths: λ₁ =553.6, λ₂ =645.5, and λ₃ =856.5 nm. This model has been derived from the generic spectral equation (Nechad et al., 2010):

where spectral coefficients A (λ), B (λ), and C (λ) were interpolated from Nechad et al.(2010) as follows (Table 4):

Below we present a derivation of the spectral ratios version of this algorithm. It is follows from Eq.24:

Two quadratic equations may be formulated, if to write an Eq.24 for all three wavelengths and to form two ratios, namely, X₁ = R_rs (λ₂)/ R_rs (λ₁) and X₂ = R_rs (λ₃)/ R_rs (λ₁), as functions of $\tilde{S}$. Two different solutions of these equations for $\tilde{S}$($\tilde{S}$ ₁ and $\tilde{S}$ ₂) are as follows:

where

with the coefficients A _i, B _i, and C _i meaning A (λ_i), B (λ_i), and C (λ_i), where i =1, 2, 3 from the Table 4.

A similar approach was applied to the semiempirical model by Shen et al.(2014). This model initially developed for individual wavelengths was adjusted for the R_rs (645.5)/ R_rs (553.6) spectral ratio (Table 3). A pair of coefficients (α and β) necessary for this model at wavelengths of 553.6 and 645.5 nm was derived by linear interpolation of values used in Table 2 by Shen et al.(2014). As noted in Shen et al.(2014), a spectral reflectance ratio approach leads to underestimated or even negative values SSC for relatively clear waters. To avoid this situation, a condition: SSC=1 g/m³, if R_rs (645.5)/ R_rs (553.6)<0.575 has been added to the model.

The currently used model (#28) has been developed as a modification of the earlier derived MODIS model (Sokoletsky et al., 2014a). The new model is derived from the measured SSC and R_rs (λ) values and it numerically close to the model values derived from the IOPs vs. SSC equations (Eqs.13-20), Eq.22 for r_rs (λ) vs. IOPs, and Eq.23 for R_rs (λ) vs. r_rs (λ).

Overall, of the 28 models listed in Table 3, the seven models use the single wavelength (in the red or NIR ranges) for the SSC retrieval while the simple spectral ratios were used in the 17 models. Han et al.(2006), Nechad et al.(2010), Zhang et al.(2010), and Bhatti et al.(2011) used some more complicated functions of the spectral reflectance. The models by Doxaran et al.(2002), Nechad et al.(2010), He et al.(2013), Liu et al.(2013), Meng et al.(2013), Shen et al.(2014), Sokoletsky et al.(2014a), and currently used MODIS model (#28) have been demonstrated to provide the best results (from statistics and visual perception points of view)(Fig. 12).

Interestingly to note that among these eight models, the He et al.(2013), Liu et al.(2013), Meng et al.(2013), Shen et al.(2014), Sokoletsky et al.(2014a), and the current MODIS-AA (#28) models have been derived from the East China natural waters (YRE with its ACA, Yellow Sea, Taihu Lake and Hangzhou Bay), while the Doxaran et al.(2002) and Nechad et al.(2010) models were obtained for the Gironde Estuary and the Southern North Sea, respectively. Overall, the currently developed MODIS-AA model (#28) seems to be the superior model. However, a total accuracy and correlation with the measured SSC values for all models, even for the best ones, is still relatively poor (Table 3); this may be explained mainly by: 1) a mismatch between literature study geographical areas and our areas; 2) a large variability of IOPs vs. SSC relationships (Table 2); and 3) inaccuracies of measurements and models themselves.

2.4.7 The SSC versus spectral R_rs (λ) ratio relationship

Figure 13 confirms the reability of using the red-togreen models by Nechad et al.(2010), Shen et al.(2014), Sokoletsky et al.(2014a), and the current MODIS-AA model (#28) for SSC retrieval in the East China natural waters. The figure also confirms that these four SSC models are robust yielding appropriate results for the datasets collected in different times in the East China coastal area.

Each match-up R_rs (λ) curve has the same color as the associated in situ R_rs (λ) curve, and is numbered following Table 5.

3 RESULT 3.1 Validation of the atmospheric correction and the water quality algorithms

A larger part of AC algorithm has been already confirmed by previous investigations. For example, a REST2 solar irradiance model by Gueymard (2008), many elements of which are used in the present study, has been validated by help of benchmark dataset. The two stream radiative transfer model used in our AC algorithm for transmittances by atmospheric molecules and aerosols has been successfully validated by comparison with numerical results derived from the Herman and Browning (1965) radiative transfer solution. Finally, transmittance of upwelling radiance was modeled following the reciprocity principle theoretically derived and validated by Sekera (1970) and Yang and Gordon (1997).

A comparison of the field-based reflectance spectra with the MODIS-derived R_rs (λ) spectra and SSC is another way to validate our AC and SSC algorithms. Unfortunately, only five in situ-remote sensing couples with the relatively small space (<700 m) and temporal (<1.4 hours) differences could be used for validation (Table 5, Fig. 14). The R_rs (λ) spectra of the field and satellite datasets demonstrate close numerical values, spectral form, and the peak position. The closeness of the spectral forms for in situ and remote sensing data stimulates using the R_rs (λ) ratios (Fig. 15) for remote sensing water quality monitoring in the study area.

3.2 Study area SSC mapping

Figure 16 shows examples of the color (false) images for the SSC in the East China coastal area, 2008 -2013. All images were produced based on the AC and SSC estimation algorithms described above and have shown reasonable values, comparable with the literature data. For example, SSC in the ECS area were found relatively small: ~1 to 70 g/m^3that is close to the values found by Kiyomoto et al.(2001), Shen et al.(2010a, 2013, 2014), Siswanto et al.(2011), and He et al.(2012, 2013). The SSC concentration values in the Taihu Lake lie in the range of ~19-150 g/m^3that corresponds to the values reported by Zhou et al.(2006), Jiang et al.(2009), Hu et al.(2010), He et al.(2012, 2013), Zhao et al.(2013), and Shen et al.(2014).

Figure 16 demonstrates the extremely high SSC values achieved in the YRE (~20-750 g/m³), Hangzhou Bay (~70-3 300 g/m³), and the Subei Bank of the Yellow Sea (~80-2 500 g/m³) that also confirmed by numerous studies, e.g., Li et al.(1994), He et al.(1999), Han et al.(2006), Salama and Shen (2010), Shen et al.(2010a, b, 2013, 2014), Shi and Wang (2010), Wang et al.(2010a), Zhang et al.(2010), He et al.(2012, 2013), Li et al.(2012), and Liu et al.(2013). Note that we rejected here some SSC values appeared too high (totally ~2% of all values), which may be caused by wrong pixel identification (e.g., in case of thin clouds, hardly identified by the AC algorithm).

4 CONCLUSION

The study revealed several important points concerning in situ observations and satellite remotesensing of water quality in the East China natural waters:

1. The area is characterized by the extremely wide range of SSC variations: from ~1 to 3 300 g/m³(Figs. 2, 3, 16). The SSC concentration values are closely related to the quality of water (Figs. 2, 3); however, such correlation was not found for Chl concentration. Also, no significant correlation was found between SSC and Chl. Thus, our study shows that the SSC cconstitutes the main impact on water quality in the area. Therefore, the water quality components estimation methods based on the close relationships between Chl and SSC (like methods described by Sokoletsky et al., 2011, 2012) are probably inadequate in the East China area;

2. Measurements and modeling results demonstrated a strong spectral dependency for the diffuse light reflectance ρ_sky (λ) and ρ_sky (λ) d_E (λ) during clear sky hours; however, such dependency is weakened with the degree of cloudiness (Figs. 4, 5). This fact allows using more elaborated radiative transfer air-water modeling (like has been described in Subsection 2.3.6) for in situ above-water measurements with the aim of better support of remote-sensing observations;

3. A comparison between different radiometric R_rs (λ) models have shown a relatively small variability (below 20% for five best models) between the models in the case of clear sky conditions in the whole 400- 900 nm spectral range (Figs. 6, 7). However, variability increases with increasing atmospheric cloudiness or diffuseness, especially in the NIR range (up to 100% and more). Thus, this modeling shows that the fieldbased estimation of R_rs (λ), like the remote-sensing estimation of this parameter, is more reliable if the radiometric data are collected during clear sky conditions.

A further decrease in discrepancies between different R_rs (λ) models has been found using the spectral R_rs (λ) ratios (Fig. 8 provides an example). In any case, the SAM model appeared as a more elaborated radiative transfer model, yielding only 1.1% negative values in the whole 350-900 nm spectral range;

4. Figures 9 and 10 demonstrate important feature of the R_rs (λ, SSC) relationships, namely, an almost exact coincidence of R_rs (λ, SSC) spectral maxima with the a (λ, SSC) spectral minima. This fact may be useful for understanding the spectral features of R_rs (λ) and for remote-sensing retrieval of IOPs. Another important feature of the R_rs (λ, SSC) relationships is an almost linear increasing with SSC at relatively small SSC; then, with increasing SSC, dependency becomes slower, and finally almost invariable with SSC (Fig. 10). This fact explains the success of different (linear and nonlinear) models for different water situations. For example, at wavelengths of 412-532, 595-676, 715, and 865 nm, a linear part of R_rs vs. SSC relationship lasts up to 6, 10, 20, and 600 g/m³, respectively (Fig. 10b), that allows to develop the linear, but still accurate algorithms for relatively clear waters. However, for a general case of widely varied SSC more complicated methods are needed;

5. The model computations have shown (Table 2) a preference of using the method of the spectral ratios (relative to the single wavelengths approach) in terms of accuracy and correlation between predicted and measured optical properties;

6. An analysis of errors caused by possible variations in the solar and viewing zenith angles has shown their small impact (within 5%-10%) on the spectral reflectance dependency on IOPs (Fig. 11). Moreover, using spectral R_rs (λ) ratios cancels an angular impact almost completely. This confirms the similar conclusions obtained before from the more elaborated numerical analysis (Sokoletsky and Gallegos, 2010);

7. The comparison study of 28 different SSC vs. R_rs (λ)[or ρ_w (λ) or R trs (λ)] models taken from literature have revealed several models (based on investigations in Chinese and non-Chinese natural waters) yielding robust results relatively close to our observational data (Fig. 12);

8. The study confirmed the reliability of using the red-to-green models for SSC retrieval (Fig. 13);

9. A new empirical SSC vs. R_rs (645.5)/ R_rs (553.6) model along with the atmospheric correction based on the analytic radiative transfer algorithm (Sokoletsky et al., 2014a) were applied for SSC mapping of inland, estuarine and coastal waters of the East China (Fig. 16). An analysis of these maps has shown their common compatibility with the published data. A validation of the atmospheric correction algorithm has been provided by comparison with the in situ measured reflectance spectra and spectral R_rs (645.5)/ R_rs (553.6) ratio (Figs. 14, 15).

Thus, it is feasible, given the rigorous nature of our approach, to consider the further application of developed AC-SSC algorithm not only for the area under study, but also for the other areas around the world.

5 ACKNOWLEDGEMENT

We thank Dr. Chris Gueymard from the Solar Consulting Services for the valuable discussion on the optical atmospheric model. We would like also to express our gratitude to the two anonymous reviewers for their extremely constructive and helpful feedback on earlier drafts of this paper.