Simultaneous atmospheric correction and quantification of suspended particulate matter in the Guadalquivir estuary from Landsat images

Earth observations (EOs) following empirical and/or analytical approaches are a feasible alternative to obtain spatial and temporal distribution of water quality variables. The limitations observed in the use of empirical approaches to estimate high concentrations of suspended particulate matter (SPM) in the estuarine water of Guadalquivir have led the authors to use a semi-analytical model, which relates the water constituents’ concentration to the water leaving reflectance. In this work, the atmospheric correction has been carried out simultaneously and the aerosol reflectance and backscattering coefficients of SPM obtained. The results are validated using in situ SPM data series provided by a monitoring network in the study area. The results show that the model allows us to successfully estimate backscattering coefficients of SPM in the estuary, differentiating clear and turbid water and using two ε( 4,5) .These considerations improve the value of R2 from 0.68 (single ε( 4,5) ) to 0.86 (two ε(4,5)) on 18 May 2009. This method could be used as a preliminary approach to obtain SPM concentration in the Guadalquivir estuary with the limitations that the model shows for turbid waters.


INTRODUCTION
River and sea dynamics interact in estuaries, which are complex aquatic systems and are among the most important areas for sediment and nutrient exchange between oceans and continents (Contreras et al., 2012).Suspended matter is an environmental indicator of the state of water quality and affects its transparency, turbidity and colour.The study of water quality variables, such as suspended sediments, chlorophyll and dissolved matter, is a key tool for characterizing the quality and dynamics of these environments, especially when regulations upstream are a constraint for the equilibrium between fresh and saltwater.
Traditionally, water quality variables have been measured in situ using conventional techniques implying direct measurements of water quality and representing a high cost both in time and money, along with their restricted spatial and temporal domain.However, the employment of remote sensing sources following empirical and/or analytical approaches has become a feasible, fast and low-cost alternative for studying the spatial distribution of water quality variables.The use of empirical algorithms presents some limitations regarding the specific from near-infrared (NIR) reflectances in the Guadalquivir estuary using Landsat satellite images, thus extending their use in Mediterranean areas.

Study area
The Guadalquivir River estuary (Southwestern Spain) stretches for 105 km between the Alcalá del Río dam, upstream, and the river mouth at Sanlúcar de Barrameda (Fig. 1).The contribution of freshwater to the estuary is subject to intense regulation by the Alcalá del Río Dam, the lowest one of a dense reservoir network throughout the 47900 km 2 contributing area, which blocks the tidal wave upstream.The sediments in the estuary are of a very fine texture due to the great length of the river and, mainly, to the extreme trapping efficiency of the dense reservoir network upstream, with an average value of 0.5-4.5 g L -1 for the suspended sediment range along the estuary with extreme values of up to 16 g L -1 .
Fig. 1 Location of Guadalquivir River estuary and of the CTD buoys.

Available turbidity data
Continuous turbidity measurements were taken in situ every 30 minutes with eight turbidimeters (Turner Designs, model Cyclops-7) using CTD Seabird Electronics SBE16plus equipment with other external sensors along the Guadalquivir axis.These eight stations were installed on navigation buoys between the river mouth and Seville harbour, strategically positioned along the estuary, providing data from February 2008.The CTD number 2 was not used because there are no date for the study period.To obtain SPM concentration with this remote monitoring network, a relationship between SPM (measured by gravimetric method) and turbidity (NTU) was established (Navarro et al., 2011).

Methodology
Atmospheric correction.The total reflectance measured by the sensor at the top of the atmosphere (TOA) was obtained from the calibrated Digital Numbers (DNs).This total reflectance at a given wavelength λ,    , can be written as the sum of several components (Gordon, 1997): where Tg λ and Tv λ are, respectively, the gaseous transmittance and the viewing diffuse transmittance from water to sensor; ρsfc (λ) is the sea-surface reflectance; ρa (λ) is the aerosol reflectance resulting from multiple scattering by aerosol in the absence of air; ρr (λ) is Rayleigh reflectance resulting from multiple scattering by air molecules in the absence of aerosol; ρra (λ) is the reflectance from the interaction between air molecules and aerosol; ρw (λ) is the water-leaving reflectance resulting from the interaction between the light and the water column.To derive the water-leaving reflectance, ρw (λ) , all other terms of equation ( 1) must be quantified.
The standard approach by Gordon and Wang (1994) assumes zero water-leaving reflectance in the NIR (ρw (5) =0).The calculation of Rayleigh reflectance ρr (λ) is well described in terms of geometry and atmospheric pressure (Gordon et al., 1988a).The coupled term ρra (λ) can be neglected at the NIR part of the spectrum (Gordon and Castano, 1987).Sea-surface reflectance, ρsfc (λ) , is estimated from the Fresnel reflectance equation.Gaseous transmittance Tg λ is calculated from ancillary data on ozone and water vapour concentrations by using the transmittance models of Goody (1964) and Malkmus (1967).The diffuse transmittance Tv λ is approximated following the model of Wang (1999).A detailed description of every step can be found in Salama and Shen (2010) and Salama et al. (2012).
Basically, equation (1) has two unknowns; the water-leaving reflectance, ρw (λ) and the aerosol reflectance, ρa (λ) .The Rayleigh corrected reflectance is computed from Aerosol scattering reflectance, ρa (λ) , is directly estimated from the aerosol ratio at two NIR bands, s and l (s short, l long wavelengths), which can be related to aerosol optical thickness and type, and is considered to be constant for the whole image.Salama et al. (2004) suggested an automated approach to determine the aerosol ratio based on eigenvector decomposition of the NIR bands: (,) = cot �0.5 tan −1 �     −  ��, where C is the correlation between Rayleigh corrected reflectance of short band (s) and long band (l) of Landsat-ETM.In this work, this approach was applied to data from Landsat-ETM bands 4 (837 nm) and band 5 (1625 nm) (Danbara, 2014).Once the value of ε(4,5) is estimated, ρw (4) can be derived from equation (1), assuming ρw 5 = 0 in accordance with the standard approach: It should be noted that this approach is not valid for high turbidity values (Kuchinkeet al.,2009) and a second parameter should be included as the ratio of water-leaving reflectances (α) at band 4 and 5, (Carder et al., 1999a,b).However, the great spatial variability of SPM in the study area, which is accentuated on certain dates coinciding with extreme turbidity events, makes it difficult to employ a spatially homogeneous ε(4,5) for the whole image.This variability justifies the fact that, at certain temporal moments, the problem is solved with different ε(4,5).In this work, different regions were identified in the corrected reflectance at band 4 versus corrected reflectance at band 5 plots to derive specific values for each group.The calculation of the ρw (4) for these areas is made not only using a specific ε(4,5) for each of them but also the most suitable algorithm.
Deriving water inherent optical properties ρw (λ) is proportionally related to physical and biological water properties according to the semi-analytical model of Gordon et al. (1988b): where ρw λ is the water-leaving reflectance normalized to the solar transmittance from sun-to-target To (λ) (Moral, 1991); l1=0.0949 is a constant coefficient; Q is the ratio between upwelling radiance and irradiance (Mobley, 1994).Assuming an isotropic light field, the value of Q is taken as being equal to π sr.The constant number 0.54 describes the fraction of transmitted light from below the water surface.a(λ) and bb(λ) are bulk absorption and backscattering, respectively.At the NIR part of the spectrum, the ρw (λ) column is assumed to be optically governed by the SPM backscattering bb(SPM)(λ) and the water absorption aw(λ) coefficients: where b*b(SPM)(λ) and Cspm are the specific backscattering coefficient and SPM concentration respectively.

Atmospheric correction
The lower mean value of ρw (4) estimated using the standard approach by Gordon and Wang (1994) is found in the image of 5 May 2010.On this day, SPM concentration values measured by the network were lower and spatially homogeneous (0.102-0.180 g L -1 ).Whereas the mean maximum value of ρw (4) , 0.126, occurred on 31 March 2009.However, when extending to turbid waters, ρw (5)  was different from 0 and the mean values of ρw both for band 4 and 5 were higher.The mean minimum value of ρw (4) on 5 May, 2010, increased to a value of 0.087 compared to that of 0.075 obtained with the standard approach by Gordon and Wang (1994).A mean maximum value of ρw (4) of 0.182 and of ρw (5) of 0.073 were reached on 31 March 2009.These ρw values are much higher than those observed by other authors using this technique for turbid waters (Ruddick et al., 2000).Making a comparative study of the ρw (4) results obtained by the standard approach by Gordon and Wang (1994) and that extended to turbid waters, it was observed that the values of the latter were around 30% higher than those of the former for all the days, both in minima, maxima and means.Negative values of ρw (4) were observed in three images (0.22% in the most unfavorable case) but they were removed from the analysis.The existence of negative values is associated with an overestimation of ρa (λ) , which underestimates the values of ρw (Mao et al., 2013), or with the existence of deeper waters.
However, the mean values of ε( 4,5) obtained for each of the images do not show the same suitability.On 5 May 2010, the line of fit associated with obtaining a mean value of ε ( 4,5) has an R 2 = 0.85, whereas for 19 August 2008, its R 2 = 0.44.On days on which this trend was observed various areas are differentiated, distinguishing between waters with a low ρw, corresponding to clear waters, and a high ρw corresponding to turbid waters.For each water type we used the corresponding ε( 4,5), thus obtaining for 19 August 2008, a mean ε( 4,5) of 0.533 and mean ρw (4) of 0.214 for turbid waters and a mean ε( 4,5) of 1.042 and a mean ρw (4) of 0.092 for clear waters.

Deriving optical properties
It was demonstrated that the application of the Gordon semi-analytical model (1988b) on all the dates was not possible over the whole estuary area.It was necessary to restrict its use to areas in which ρw (4) < 0.54l1 Q, to prevent negative SPM backscattering coefficients (bb (SPM)).Also, values of ρw (4) which although being lower are very close to 0.54 l1Q give values of bb (SPM) of a high order of magnitude.In the study area, in which the SPM concentration range is high and varied, there are many areas with high ρw (4) which cannot be modelled.
Figure 2 (a) shows the derived bb(SPM) along the estuary on 25 December 2008 for the 90% of the estuary area.Values of bb (SPM) were observed with a minimum of 0.158 m -1 and a maximum of 13.47 m -1 .Large percentage, 83%, of the bb (SPM) values were lower than 2 m -1 .On days when the SPM concentration was low, the bb(SPM) could be obtained for a larger area, reaching an area of 95% on 5 May 2010.
Figure 2 (b) shows the relationship between SPM measured and derived bb(SPM) coefficients at the CTDs for all the days studied.The image taken on 11 February 2009 has been removed from the study due to the extreme conditions found which are at the limit of the optimum working range of CTDs.A positive linear relationship can be observed on most days, which allows us to estimate SPM concentration.Figure 3 shows the measured versus estimated SPM concentration on four days using a mean value of ε(4,5) (Fig. 3(a)) and two values (Fig. 3(b)).It is clear that the employment of two ε(4,5) improves the estimation of the SPM concentration, with better R 2 and RMSE (see Table 1).On 18 May 2009, R 2 varies from 0.68 for the single value solution to 0.86 for the two solution with RMSE of 46 g m -3 and 35 g m -3 respectively.

CONCLUSIONS
The standard approach by Gordon and Wang (1994) which assumes that ρw (5) =0, underestimates the water leaving reflectance at band four ρw (4) , by approximately 30% with respect to those obtained when we do not neglect the NIR reflectance.
The semi-analytical model of Gordon (1988b), employed in band 4 to obtain the water's properties, is valid for a range of ρw (4) , so that ρw (4) < 0.54l1Q.In this case, with high SPM concentrations, it was not possible to apply over the whole study area, and its range was between 80% and 95% of the area for the days 19 August 2008 and 5 May 2010, respectively.This fact justifies the need to define new calibration constants for this model which could be adapted to very turbid waters.
The use of two values of ε( 4,5) associated with turbidity level results in a better accuracy in the retrieved SPM values.The use of the two ε( 4,5) approach improves the results as R 2 = 0.86 and RMSE = 35 g m -3 in comparison to R 2 = 0.68 and RMSE = 46 g m -3 of the single ε( 4,5) value.

Fig. 3
Fig. 3 (a) Measured versus estimated SPM concentration on four days using a single mean value of ε(4,5); (b) measured versus estimated SPM concentration on four days using two mean value of ε(4,5)

Table 1
Values of R 2 and RMSE obtained using a single mean value of ε( 4,5) and two value of ε(4,5).