1. Introduction
The US Environmental Protection Agency [
1] ranked siltation and its attendant turbidity as one of the most widespread pollutants in the United States. It affects rivers, lakes, reservoirs, and coastal waters, adversely affecting aquatic habitat, drinking water treatment facilities, recreational use, and commercial navigation and fisheries. High levels of sediment in the water column reduce spawning grounds and reduce the availability of food supplies for aquatic life. High concentrations of sediment may lead to rapid deposition, causing changes in coastal or fluvial morphology, which may lead to navigational hazards or increase flood risks for the region. Regions of low turbidity due to low suspended sediment concentrations) may experience rapid erosion leading to entrenchment and bank failure in a fluvial system or to beach erosion and barrier island collapse in the coastal environment.
Remote sensing of turbidity has a long history in the scientific literature. Descriptions of qualitative analyses relating Landsat Multispectral Scanner System (MSS) surface radiance to turbidity and suspended sediment concentration appear in the literature as early as 1973 (e.g., [
2,
3,
4]). Quantitative algorithms began to appear in the 1970s, and ranged in complexity from linear relationships between MSS Band 5 radiance and suspended sediment concentration to polynomial fits to the ratio of MSS radiances in different bands (see [
5] and references therein). More rigorous, model-driven algorithms were in development by the late 1970s (e.g., [
5]) with the incorporation of diffuse reflection models developed earlier in the decade (e.g., [
6,
7,
8]). Curran and Novo [
9] summarized early efforts, though several significant studies post-date this work (e.g., [
10,
11,
12,
13,
14]). More recently, good results have been obtained in retrieving turbidity using off-nadir Moderate Resolution Imaging Spectroradiometer (MODIS) imagery (e.g., [
15]). They indicate a near linear relation between
in situ turbidity and red reflectance when applied to individual regions. However, formulation of a global algorithm applicable to all regions has not yet been demonstrated.
Spaceborne lidars have been used in earth observation since the mid-1990s when BALKAN flew on the space station MIR [
16]. A large fraction of these lidars have fallen on water. Unlike their airborne counterparts, spaceborne lidars have not possessed the vertical resolution necessary to profile the water column. This paper offers to add value to existing lidar datasets by using the depth-integrated lidar backscatter to extract information about turbidity or scattering within the water column.
The use of a 532 nm laser to study water turbidity was pioneered in the 1980s with theoretical work by Gordon [
17], and by Phillips
et al. [
18]. Experiments into the use of a ship-mounted laser (532 nm) for measuring the optical properties of seawater were carried out by Ivanov
et al. [
19] in the Soviet Union. Empirical studies using the Australian Weapons Research Establishment Laser Airborne Depth Sounder (WRELADS) sensor by Phillips
et al. [
20], and by Billard [
21] showed that estimates of both backscattering and effective attenuation were possible from depth sounding lidars. The 532 nm channel is of particular interest in both turbidity monitoring and in other hydrologic applications, because it lies close to the wavelength achieving maximum penetration into the water column [
22].
Spaceborne lidars have been used very little for oceanography. Lancaster
et al. [
23] used the GLAS (Geoscience Laser Altimeter System) lidar on the ICESat (Ice, Cloud, and land Elevation Satellite) platform to detect water surface reflectance and compared that with the predicted theoretical reflectance from QuikSCAT (Quick Scatterometer) wind speeds. Hu
et al. [
24] used the CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) lidar on the CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations) satellite to measure wind speed, and compared the results to AMSR-E (Advanced Microwave Scanning Radiometer-EOS) wind speed measurements. They used the depolarization ratio available from CALIOP to remove the sea foam effect from the surface measurement.
The study detailed in this paper also uses the CALIOP lidar, which has a 30 m vertical resolution near the surface. CALIOP has a 90 m horizontal resolution, with centers spaced 333 m along track. CALIPSO flies in NASA’s A-Train constellation, following Aqua by 73 ± 43 seconds. [
25] As an active sensor, it is capable of making night-time measurements, and the ability to separate atmospheric returns from surface returns makes local atmospheric correction possible.
The overall approach is to estimate subaqueous backscatter as the residual term within an expression of total depth-integrated attenuated backscatter. Specifically, total, depth-integrated attenuated backscatter is formulated as the sum of all surface and subsurface effects. With total backscatter obtained from CALIOP, and other surface components modeled, the subsurface backscatter can be retrieved. Sensitivity analysis on the accuracy of the CALIOP data and on individual model components is used to delineate the limits and potential accuracy of the subsurface backscatter retrieval. Theoretical sensitivity analysis is followed by a test case over the Tampa Bay region of Florida.
2. Model Basis
The depth-integrated attenuated backscatter (at a wavelength of
λ, in nm) received by a satellite from the surface bins can be represented as an attenuated sum of scattering returning from components of the target:
where
Tλ is the atmospheric transmittance (including aerosol effects, and along the lidar look direction), the variables
,
,
,
, and
represent, respectively, the integrated backscatter due to specular reflection of the laser from the water surface, laser light scattered by foam on the water surface, the specular reflection of sunlight from the water surface, sunlight scattered by foam on the water surface, and particles in the underwater environment. The magnitudes of
,
,
, and
are all primarily controlled by the wind and wave field.
is controlled by the number and properties of particulates present in the water column. Other possible components beyond the scope of this study and not listed in (1) include bottom scattering effects, scattering from surface slicks, and scattering by spray in regions of high winds. Details of each component are described below.
The goal of the modeling presented here is to allow the estimation of from satellite measurements of , and , and the modeling of the other components. Surface wind speed and visibility estimates are used as inputs into the models.
2.1. Atmospheric Correction
A principal factor affecting the radiometric surface measurement is the atmospheric transmittance, which can be further subdivided into molecular, aerosol, and cloud transmittance. One advantage of using lidar technology compared with a passive sensor is the ability to separate the scattering due to the atmospheric effects from the scattering at the surface based on the time of the photons’ arrival at the satellite. The CALIPSO release products include an aerosol layer optical depth and a cloud layer optical depth product. These data, combined with temperature and pressure data from NASA’s Global Modeling and Assimilation Office, and included in the CALIPSO release product, were used in the current work to develop a LOWTRAN 7 atmospheric model.
In the current study, optical depth of the cloud layers and aerosol layers were accounted for using the CALIOP data. However, optical depths of the clear air interstices, not available from the CALIOP archive, were evaluated using the LOWTRAN 7 model, resulting in a continuous profile of optical depth throughout the atmospheric column. The optical depths for the layers were summed (taking into account that the aerosol and cloud layers sometimes overlap), and a resulting total atmospheric optical depth and transmittance were obtained. To calculate the clear-air optical depths, default settings were used for the LOWTRAN 7 model, except for the input of relative humidity, temperature, pressure and ozone concentration data extracted from the Level 1 and Level 2 CALIOP products, and surface visibility was interpolated from historical airport records from around Tampa Bay. A sensitivity analysis showed that assuming a constant surface visibility of 23 km produced errors in the subaqueous integrated backscatter of less than 5%.
2.2. Specular Surface Reflection
Specular reflection from the water surface is by far the largest contribution to the total integrated attenuated backscatter. Specular reflection decreases as wave steepness increases and is extremely large in very calm waters. Hu
et al. [
24], citing Platt [
26], Menzies
et al. [
27], and Tratt
et al. [
28], modeled the specular surface reflection as:
where
is the Fresnel specular reflection coefficient (
,
),
is the wave slope variance, and
is the zenith angle of the sensor. In the present study, the small angle approximation is not made, as it was made by Hu
et al., because it can produce significant errors at low wind speed even for small angles such as the 0.3° pointing angle for CALIOP. Hu
et al.’s composite model for the wave slope variance as a function of wind speed is employed:
Although this model provides a reasonable guess at the backscattering when given a known wind speed, it is very sensitive to the wind speed. To overcome this sensitivity, one can exploit the fact that the only wavelength dependent term in (2) is
ρλ, and also the fact that infrared light does not penetrate the ocean surface;
i.e.,
is zero. The specular surface reflection in the green can then be formulated in terms of the total integrated attenuated backscatter in the infrared:
where
,
, and
are described in the subsequent sections. The resulting model is only weakly dependent on the accuracy of the wind speed estimate, especially for low wind speeds and for night measurements, when
and
are zero.
An added benefit of this infrared-based approach is that deviations from the fully-developed sea models of Cox and Munk [
29] or Bréon and Henriot [
30] are accounted for implicitly, without the need for estimating fetch, wind direction, or bottom topography. This is particularly important in coastal zones, where sheltered waters can have very short fetches, and shallows can cause rapid and localized steepening of waves.
2.3. Lambertian Scattering from Foam
At relatively high wind speeds, the scattering of the lidar beam from whitecaps and foam streaks on the water surface can be a significant factor, particularly as the wave surface return becomes darker as the waves steepen. At wind speeds higher than about 10 m/s, the magnitude of foam scattering approaches that expected for subsurface scattering in low turbidity natural waters. As noted above, Hu
et al. [
24] used the depolarization ratio from the CALIOP lidar to model the foam contribution to the backscattering. However, depolarization may occur during scattering from foam or from hydrosols, so this type of model may eliminate the potential to retrieve information from the subsurface. The foam is therefore modeled as a Lambertian scattering process, from foam covering a fractional area,
W, described by Callaghan and White [
31]:
Moore
et al. [
32] modeled the reflectance of foam as a function of wind speed. In this model the reflectance of the foam is expressed as an “additional” reflectance, representing the increased reflectance of the ocean surface due to the foam. In the visible wavelengths, including 532 nm and at the upper limit of 670 nm, they expressed this additional reflectance as:
For the near infrared at 860 nm they expressed the additional reflectance, RSAR
860, as
The laboratory data of Frouin
et al. [
33] suggests an exponential decrease in the reflectance of foam with wavelength (λ, in nm) as shown in
Figure 1:
At a given wind speed, the coefficient A and the decay constant k can be estimated by fitting to Equations (6) and (7). Repeating this process over a range of wind speeds, and using least squares polynomial regression, one can estimate A[U] and k[U] as follows.
For wind speeds less than about 10 m/s, A[U] and k[U] can be approximated by their constant first terms.
Figure 1.
Data of Frouin
et al. [
33], which suggests an exponential function of wavelength. The correlation is not significant, with only four data points, but it suggests that exponential decay is a starting point upon which we can construct a model.
Figure 1.
Data of Frouin
et al. [
33], which suggests an exponential function of wavelength. The correlation is not significant, with only four data points, but it suggests that exponential decay is a starting point upon which we can construct a model.
Using these results, and multiplying by
to convert between reflectance and Lambertian scattering, one can write the expression for the scattering due to the foam. Because the model presented above calculates the additional reflectance due to the foam, an estimate of the sea surface scattering without foam is also included or Equation (2). It is possible to use Equation 10(b) for the 532 nm lidar by replacing −1,064 k with −532 k in the exponential, and the same results will be obtained. Equation 10(a) can be used when efficiency is a concern. A graphical representation of these models is presented in
Figure 2.
Figure 2.
Model of backscatter due to Lambertian scattering from foam. Three wavelengths are shown. The black curve shows the model of Moore
et al. [
32] for visible wavelengths up to 670 nm (Equation (6)). This model is used for the 532 nm laser light. The blue curve is the model of Moore
et al. [
32] for 860 nm shortwave infrared (Equation (7)). The red curve shows an extrapolated model for 1,064 nm, created as described in the text as an exponential decay with wavelength perfectly fitting the two shorter-wavelength models. (Equation (8), where
λ = 1,064 nm.)
Figure 2.
Model of backscatter due to Lambertian scattering from foam. Three wavelengths are shown. The black curve shows the model of Moore
et al. [
32] for visible wavelengths up to 670 nm (Equation (6)). This model is used for the 532 nm laser light. The blue curve is the model of Moore
et al. [
32] for 860 nm shortwave infrared (Equation (7)). The red curve shows an extrapolated model for 1,064 nm, created as described in the text as an exponential decay with wavelength perfectly fitting the two shorter-wavelength models. (Equation (8), where
λ = 1,064 nm.)
2.4. Sun Glint
At low solar zenith angles, and particularly for lidar sensors which may use low-power micropulse sensing, the presence of specular reflection of solar radiation, or sun glint, may affect the measurement. To estimate the sun glint, the following model after Kay
et al. [
34] is employed:
W is the fractional area covered by foam,
Eλ is the extra-atmospheric solar irradiance,
ρ(ω,λ) is the Fresnel reflection coefficient of the water surface for the given solar and sensor geometry,
p(zx,zy) is the probability of a given surface element having the required surface slope to specularly reflect light from the sun into the detector,
θs is the solar zenith angle, and
θv, is the sensor viewing angle. The ω and β are angular measures accounting for the geometry of the sun and sensor, given the azimuths of the sun,
𝜑s, and sensor,
𝜑v:
Kay
et al. [
34] provide the Cox and Munk [
29] expression for
p(zx,zy) as well as a more recent version by Bréon and Henriot [
30] which presents coefficients with reduced uncertainty. In this study, the version of Bréon and Henriot is used.
In the absence of foam, sun glint, while dependent on wind speed through the probability density function for wave slopes, is relatively insensitive to wind speeds over the normal range, changing only by 0.3% between 0 and 30 m/s. The error in the scattering estimate caused by sun glint is also very low compared with the intensity of the lidar backscattering, because only approximately 0.006% of the sun’s energy falls within the bandwidth of the CALIOP detector. In the presence of foam, the sun glint falls off as the wind speed increases because less glint is possible from foam-covered slopes. This model assumes that foam is distributed over the surface in a pattern uncorrelated with water surface slope. Even in the case where this assumption is invalid, the estimated sun glint for 3.7 m/s represents an upper bound on the sun glint for wind speeds greater than this value, as it represents a foamless case. For winds less than 3.7 m/s, this assumption is not necessary as foam is assumed to be absent.
The sun glint model is relatively insensitive to uncertainty in the parameters. Using published [
30] uncertainties in the parameters for the wave slope probability, and 10% uncertainties in the parameters of the foam model as above, the model exhibits at most 3% error in the predicted sun glint.
2.5. Sunlit Foam
Sunlit foam also is modeled as a Lambertian scatterer. Using a method similar to Moore
et al. [
32], scattering is modeled as an additional scattering term dependent upon the foam reflectance:
The first term in Equation (13) represents the surface return in the absence of foam, while the second term represents the additional scattering due to foam.
W and
have been previously defined in Equations (5) and (11), respectively.
RSARλ is the additional foam reflectance, either Equation (6) or (8) depending on wavelength,
Eλ is the extra-atmospheric solar irradiance,
Ω is the solid angle of the lidar detector, A is the area of the ground spot,
I is the laser energy,
tp is the pulse length, and
Tλ is the atmospheric transmittance along the lidar look direction. See
Table 1 for the lidar-specific parameters applicable to CALIOP. The only wind-speed dependence for this term is the areal coverage of the foam. For all wind speeds less than 30 m/s, the foam contribution is an order of magnitude less than the sun glint.
Table 1.
CALIOP lidar-specific parameters for Equation (13).
Table 1.
CALIOP lidar-specific parameters for Equation (13).
Parameter | Value | Note |
---|
Ω | 1.26 × 10−13 sr | |
A | 3.85 × 103 m2 | |
I | 1.1 × 10−1 J | |
tp | 2.0 × 10−8 s | |
θv | 5.236 × 10−3 (0.3°) | (June 2006–November 2007) |
5.236 × 10−2 (3.0°) | (November 2007–) |
The sensitivity to model parameters in the sunlit-foam model is dependent on wind speed. The parameters of the sunlit foam model are the same as those of the lidar foam model, and relate to the coverage and reflectance of the foam. However, the sunlit-foam model is less affected by the model parameters than the lidar foam model. At low-to-moderate wind speeds, the model is fairly insensitive to the model parameters, with uncertainties less than a factor of two for wind speeds less than about 22 m/s.
2.6. Subaqueous Scattering
The full model for subaqueous scattering (
) is obtained by rewriting Equation (1) and combining with Equation (4), or:
Further combining Equation (14) with Equations (6, 8, 10, 11 and 13) yields the expanded form:
For nighttime overpasses, sun glint and sunlit foam in both wavelengths are zero simplifying to:
2.7. Comparison of Backscatter Terms
Figure 3 shows model results for the magnitude of the components of the total integrated attenuated backscatter without subsurface scattering, to give a sense of the relative importance of each term. A solar zenith angle of 30° is assumed. The specular wave scattering shown is from Equation (2), in contrast to the method outlined in
Section 2.2 to allow a theoretical dependence on wind speed for comparison to the other terms. In the case study presented in
Section 3, subsurface scattering ranged from 0.0 to 0.13 sr
−1.
Figure 3.
Theoretical CALIOP scattering off water surface as a function of wind speed for 532 and 1064 nm wavelengths. Wave scattering is modeled using Equation (2). Foam scattering and sun-on-foam are modeled by Equation (10). Sun glint is modeled as Equation (11). Graphs illustrate the relative importance of the different components to the total integrated scattering. The solar zenith angle is assumed to be 30°.
Figure 3.
Theoretical CALIOP scattering off water surface as a function of wind speed for 532 and 1064 nm wavelengths. Wave scattering is modeled using Equation (2). Foam scattering and sun-on-foam are modeled by Equation (10). Sun glint is modeled as Equation (11). Graphs illustrate the relative importance of the different components to the total integrated scattering. The solar zenith angle is assumed to be 30°.
2.8. Uncertainties
Uncertainty propagation analysis allows the evaluation of the sensitivity of the model to input error. There are five major inputs into the model: First, an estimate of wind speed is required to estimate the foam component. Second, transmittances are needed for both the green and the infrared channels. Third, the total integrated attenuated backscatter from the lidar measurements is needed in both green and infrared. Uncertainty is introduced both during the measurement process and also during any numerical integration with depth.
2.8.1. Wind Speed Estimate
The use of Equation (4) instead of Equation (2) for the specular surface return results in a modeled subaqueous integrated backscatter that is relatively insensitive to the wind speed estimate.
Figure 4 shows this lack of sensitivity, which even at very high wind speeds never exceeds 0.002 s m
−1 sr
−1.
Figure 4.
Error in subaqueous integrated backscatter due to 1 m/s overestimate in wind speed estimate. Because the wind speed estimate is only used to estimate foam cover and reflectance, the error at low wind speeds is negligible.
Figure 4.
Error in subaqueous integrated backscatter due to 1 m/s overestimate in wind speed estimate. Because the wind speed estimate is only used to estimate foam cover and reflectance, the error at low wind speeds is negligible.
2.8.2. Transmittance Estimate
Two transmittance estimates are used in the atmospheric correction of the lidar data; at 532 nm and at 1,064 nm. The accuracies of these estimates are important considerations in determining the accuracy of the subaqueous integrated backscatter estimate. From Equation (15), it can be shown that:
Equation (17) shows that the sensitivity to the transmittance increases linearly with the total backscatter.
Figure 5 shows the dependence upon both the total integrated backscatter and the transmittance. It is clear from both Equation (17) and from
Figure 5 that the accuracy of the transmittance estimate plays an important role in determining the uncertainty of the final result, and that low transmittance data may provide largely meaningless results.
Because of the opposing signs of (17a) and (17b), positively correlated errors in T532 and T1064 will tend to cancel. This will tend to limit the error introduced by the surface visibility parameter to the model.
Figure 5.
Error in subaqueous integrated backscatter due to 0.01 overestimate in the 532 nm transmittance estimate. The error due to an overestimate in the 1064 nm transmittance is proportional, with a constant of proportionality of . The error is dependent upon the total integrated backscatter received; the two curves represent low (0.01 sr−1, blue) and high (0.1 sr−1, red) integrated backscatters.
Figure 5.
Error in subaqueous integrated backscatter due to 0.01 overestimate in the 532 nm transmittance estimate. The error due to an overestimate in the 1064 nm transmittance is proportional, with a constant of proportionality of . The error is dependent upon the total integrated backscatter received; the two curves represent low (0.01 sr−1, blue) and high (0.1 sr−1, red) integrated backscatters.
2.8.3. Integrated Backscatter
The attenuated backscatter measured by the lidar sensor is numerically integrated to obtain the total integrated attenuated backscatter. There is uncertainty associated both with the original attenuated backscatter measurements and with the numerical integration scheme used. The sensitivity of the model to these uncertainties can be summarized by the following:
Like the sensitivity to the transmittance, the sensitivity to the total integrated attenuated backscatter is dependent on the transmittance (
Figure 6). The error in the subaqueous integrated backscatter is always greater than the error in the total integrated attenuated backscatter, and is greater by more than a factor of four for transmittances less than 0.5. This sensitivity makes the accuracy of the total integrated attenuated backscatter input at least as important as the accuracy of the transmittance estimate in determining the uncertainty in the final result.
Figure 6.
Error in the subaqueous integrated backscatter for a 0.001 sr−1 error in the total integrated attenuated backscatter. This level of uncertainty was observed in the CALIOP data when comparing different numerical integration schemes. The errors associated with uncertainties in both bands are shown.
Figure 6.
Error in the subaqueous integrated backscatter for a 0.001 sr−1 error in the total integrated attenuated backscatter. This level of uncertainty was observed in the CALIOP data when comparing different numerical integration schemes. The errors associated with uncertainties in both bands are shown.
As for the transmittances, the opposing signs (and roughly equal magnitude; ) of (18a) and (18b) will cause the effects of positively correlated errors in γ532 and γ1064 to cancel. This should help to ameliorate the effect of the choice of numerical integration scheme, which might introduce systematic errors.
4. Conclusions
A theoretical method has been presented for estimating subaqueous integrated backscatter from the near-nadir viewing satellite lidars such as CALIOP. External inputs to the algorithm are principally wind speed and surface visibility. The algorithm takes into account specular reflection of laser light at the water surface, laser scattering by wind-generated foam as well as sun glint and solar scattering from the foam. An important feature of the algorithm is the use of the infrared backscattering as a basis for predicting the specular green backscattering. This allows the implicit correction for water surface geometry that may not follow simplistic model assumptions.
Figure 10.
Scatter plot of CALIOP 532 nm integrated subsurface backscatter and MODIS 645 nm remote sensing reflectance. The Pearson’s correlation coefficient is 0.11. The exceptionally low points (γu < −0.05 sr−1) are all from 07 May 2007, from points of low transmittance. Fourteen outliers identified with Peirce’s criterion with γu > 1.0 sr−1 are not shown.
Figure 10.
Scatter plot of CALIOP 532 nm integrated subsurface backscatter and MODIS 645 nm remote sensing reflectance. The Pearson’s correlation coefficient is 0.11. The exceptionally low points (γu < −0.05 sr−1) are all from 07 May 2007, from points of low transmittance. Fourteen outliers identified with Peirce’s criterion with γu > 1.0 sr−1 are not shown.
Theoretical analysis indicates that the model is insensitive to wind speed, but relatively sensitive to errors in atmospheric transmittance and total integrated attenuated backscatter. However, errors in the transmittance or total attenuated integrated backscatter that are positively correlated between the two bands will tend to cancel out because of the reversed sign of the sensitivity.
In a case study, the method was applied to nighttime CALIOP data over Tampa Bay, using interpolated wind speed and visibility data from airports and buoys, and comparison was made to MODIS 645 nm remote sensing reflectance. Good correlation was found on six of seven dates with relatively high atmospheric transmittances, although the correlation decreases when data were composited over all dates. Results were found to be very insensitive to the wind speed and visibility inputs as predicted. The implication is that the 532 channel responds differently to the seasonally varying biogeochemical scatterers in the water column. Thus, while the single CALIOP 532 total integrated attenuated backscatter contains information about the subsurface, more channels may be needed to capture the full range of scattering affecting turbidity, at least within Tampa Bay. CALIOP was designed to measure principally optical properties of the atmosphere, but the current analysis indicates a potential for extracting subsurface scattering and turbidity as well.