1. Introduction

The irradiance illuminating the water surface is composed of two spectrally and geometrically distinct components: the direct irradiance from the Sun, $E_{dd}$, which is an almost parallel beam of light from the direction of the Sun; and the diffuse irradiance of the sky, $E_{ds}$, formed by rays with incidence angles covering the upper hemisphere. A number of measurement techniques and simulation tools exist to determine both components experimentally or by modeling.

In water, the two components cannot be determined reliably, not by measurement or by simulation, because the water surface is almost never perfectly flat. Waves, ripples, and foam incline the water surface with spatially highly variable slopes, and this locally variable geometry of the water surface is changing quickly in time because of wind and currents. Because the surface geometry determines the refraction angle, the angular distribution of downwelling radiation is changing for an observer in water locally and temporally in an unpredictable way. The induced changes depend on wind speed, solar elevation, and depth [1,2]. Observations show that variations are typically of the order of 20% to 40% in the upper few meters concerning intensity [3,4] and 5% concerning spectral shape across the visible [5], but flashes can increase intensity up to a factor of five [6]. Because of this huge variability, in-water measurements of irradiance, and consequently of reflectance, are quite challenging [7,8]. For this reason, some measurement concepts, e.g., the widely used Ocean optics protocols [9], recommend to make all incident irradiance measurements above the surface.

Radiative transfer models usually account for the influence of the water surface on the underwater light field by an empirical relationship between wind speed and the inclination of the wave facets [10–12]. Because of the unpredictable behavior of individual facets, these models can describe the influence on irradiance only as statistical averages but not for individual measurements. The analytic model developed in this paper uses a different approach that is not based on parameters describing wind and waves, but instead it is based on the induced changes of $E_{dd}$ and $E_{ds}$ intensity. $E_{dd}$ and $E_{ds}$ are treated as two spectrally well-defined light sources with unknown intensities. Their spectral shapes at depth $z$ are calculated individually using an analytic model, and their intensities are fit parameters during data analysis. In this way, inverse modeling allows to analyze irradiance measurements even for geometrically undefined radiance distributions. Results of data analysis are a number of model parameters like sensor depth and concentrations of water constituents as well as separation of $E_{dd}$ and $E_{ds}$.

The depth dependency of irradiance is commonly parameterized by the diffuse attenuation coefficient, $K_d$, which is an apparent optical property. The new approach replaces $K_d$ by an inherent optical property, $a+b_b$, where $a$ is the average absorption coefficient and $b_b$ is the average backscattering coefficient of the water column between surface and sensor. The model was validated using the well-established radiative transfer program HydroLight [10] and implemented into the public-domain software WASI [13–15] for forward calculation and inverse modeling.

2. Parameterization of Irradiance

The downwelling irradiance, $E_d$, is the sum of a direct ($E_{dd}$) and a diffuse ($E_{ds}$) component. In this study, these components are defined according to their spectral shapes at the location of the observer in air or in water: $E_{dd}(\lambda ,z)$ represents the spectral irradiance of the Sun disk for an observer at depth $z$, and $E_{ds}(\lambda ,z)$ the average irradiance of the sky excluding the Sun. Upwelling radiation that is reflected at the water surface or scattered in the water in downward direction is neglected. Hence,

The diffuse component at depth $z$ is related to that below the surface as follows:

The direct component is attenuated along a path with length $z/\cos {\theta }_{\text{sun}}^{\prime }$:

The irradiance components beneath the surface ($0-$) are related to those in air ($0+$) through $E_{dd}(\lambda ,0-)=E_{dd}(\lambda ,0+)(1-{\rho }_{dd})$ and $E_{ds}(\lambda ,0-)=E_{ds}(\lambda ,0+)(1-{\rho }_{ds})$, where the reflectance factors ${\rho }_{dd}$ and ${\rho }_{ds}$ describe the losses of the direct and diffuse components of the downwelling irradiance at the air-water interface, respectively. Typical values are ${\rho }_{dd}\approx 0.02$ and ${\rho }_{ds}\approx 0.06$ for small Sun zenith angles. The actual values are calculated as function of the Sun zenith angle as follows:

Equation (4) is the Fresnel equation for unpolarized radiation [16]. Equation (5) was derived from radiative transfer simulations as described below. Both equations are valid for a plane water surface. Waves and foam alter locally the inclination of the water surface, and the resulting effective values of ${\theta }_{\text{sun}}^{\prime }$ lead to changes of ${\rho }_{dd}$ and ${\rho }_{ds}$, which can be calculated only for well-known wind speed and wind direction as statistical averages [11,12]. However, as the wavelength dependencies of ${\rho }_{dd}$ and ${\rho }_{ds}$ are small, the impact of these changes on $E_{dd}$ and $E_{ds}$ can be described by wavelength-independent correction factors ${\epsilon }_{dd}$ and ${\epsilon }_{ds}$, which change $f_{dd}$ and $f_{ds}$ toward $(1+{\epsilon }_{dd})f_{dd}$ and $(1+{\epsilon }_{ds})f_{ds}$. Consequently, the parameterization of $E_d$ through Eq. (1) accounts for changes of surface reflections induced by a wind-roughened water surface.

A suitable model to calculate $E_{dd}(\lambda ,0-)$ and $E_{ds}(\lambda ,0-)$ for undisturbed illumination geometry has been developed by Gregg and Carder [17]. This widely used model is adopted here. The equations used in the software implementation of WASI 4 are recalled; for more details, see Gregg and Carder [17]. The direct component of downwelling irradiance is calculated as follows:

The diffuse component of downwelling irradiance is given by

The aerosol forward scattering probability, $F_a$, is calculated using the empirical equation $F_a=1-0.5\exp \{(B_1+B_2\cos {\theta }_{\text{sun}})\cos {\theta }_{\text{sun}}\}$ with $B_1=B_3[1.459+B_3(0.1595+0.4129B_3)]$, $B_2=B_3[0.0783+B_3(-0.3824-0.5874B_3)]$, $B_3=\ln (1-〈\cos {\theta }_{\text{sun}}〉)$. The asymmetry factor of the aerosol scattering phase function is calculated as $〈\cos {\theta }_{\text{sun}}〉=-0.1417\alpha +0.82$ when $\alpha$ is in the range 0 to 1.2, and set equal 0.82 for $\alpha <0$ and 0.65 for $\alpha >1.2$.

The atmospheric transmittance spectra are calculated as follows:

Aerosol is parameterized in terms of aerosol optical thickness, ${\tau }_a(\lambda )=\beta {(\lambda /{\lambda }_a)}^{-\alpha }$, and aerosol single scattering albedo, ${\omega }_a=(-0.0032\mathrm{AM}+0.972)\exp \{3.06\times {10}^{-4}\mathrm{RH}\}$. The Angström exponent $\alpha$ determines the wavelength dependency, and the turbidity coefficient $\beta$ is a measure of aerosol concentration. The reference wavelength ${\lambda }_a$ is set to 550 nm. $\alpha$ typically ranges from 0.2 to 2, and $\beta$ ranges from 0.16 to 0.50. $\beta$ is related to horizontal visibility $V$ and aerosol scale height $H_a$: $\beta ={\tau }_a(550)=3.91H_a/V$. Typical values are 8 to 24 km for $V$, and 1 km for $H_a$. The parameters of ${\omega }_a$ are air mass type, $\mathrm{AM}$, which ranges from 1 (typical of open-ocean aerosols) to 10 (typical of continental aerosols), and relative humidity, $\mathrm{RH}$, with typical values from 46% to 91%. $WV$ is the total precipitable water vapor content (in units of cm) in a vertical path from the top of the atmosphere to the surface.

The atmospheric path length is $M=1/\cos {\theta }_{\text{sun}}+a{(90°+b-{\theta }_{\text{sun}})}^{-c}]$. The numerical values used by Gregg and Carder ($a=0.15$, $b=3.885°$, $c=1.253$) were replaced by updated values $a=0.50572$, $b=6.07995°$, $c=1.6364$ from Kasten and Young [18]. $M^{\prime }=MP/(1013.25\mathrm{mbar})$ is the path length corrected for nonstandard atmospheric pressure $P$, and $M_{\mathrm{oz}}=1.0035/{({\cos }^2{\theta }_{\text{sun}}+0.007)}^{1/2}$ is the path length for ozone.

Aerosol optical thickness (${\tau }_a$) and the absorption coefficients of ozone ($a_{\mathrm{oz}}$), oxygen ($a_o$), and water vapor ($a_{\mathrm{wv}}$) are wavelength dependent; the other parameters (${\omega }_a$, ${\tau }_a$, $M$, $M$, $H_{\mathrm{oz}}$, $M_{\mathrm{oz}}$, $WV$) are independent of $\lambda$. Gregg and Carder [17] list the values of $a_{\mathrm{oz}}(\lambda )$, $a_o(\lambda )$, $a_{\mathrm{wv}}$($\lambda$) for the range 400–700 nm in 1 nm intervals. I extended these spectra to the range 300–1000 nm, as described in Section 3.

Because the ratio of direct to diffuse irradiance, $r_d$, is the key parameter of $E_d$ variance in water [5], analytic equations are derived that express $r_d$ as function of the parameters of the irradiance model. Just below the water surface, the ratio is given by

Equation (16) shows that the wavelength dependency of $r_d(0-)$ is determined by the scattering components of the atmosphere but not by its absorbing components. Consequently, the distinctive spectral gradients of $E_d$, which are caused by the extraterrestrial solar irradiance and the absorbing components of the atmosphere, are not present in $r_d(0-)$. Rather $r_d(0-)$ has a smooth spectral shape, which is increasing almost linearly from the shortwave to the longwave {see Fig. 3(a) in Gege and Pinnel [5]}. For depth $z$ the following relationship is obtained:

 

Fig. 1. Atmospheric transmission after absorption by ozone, oxygen, and water vapor.

Download Full Size | PDF

 

Fig. 2. Effective absorption coefficients of the atmospheric components ozone, oxygen, and water vapor.

Download Full Size | PDF

 

Fig. 3. Comparison of analytically estimated sensor depths with results from inverse modeling.

Download Full Size | PDF

Hence, the wavelength dependency of $r_d(0-)$ is altered at depth $z>0$ by a factor whose spectral shape is determined by $K_{ds}(\lambda )$ and $K_{dd}(\lambda )$. For the concentrations and depths studied here ($z<5\mathrm{m}$), these changes are relevant only above 700 nm {see Fig. 3(b) in Gege and Pinnel [5]}.

3. Optical Properties of the Atmosphere

The adopted model of downwelling irradiance is a semi-empirical parameterization of the radiative transfer within the atmosphere. In particular, absorption of radiation is simplified by ignoring the temperature and pressure dependencies of the absorption coefficients of the atmospheric gases. These coefficients are replaced by “effective” absorption spectra representing averages over the vertical profile. Gregg and Carder [17] provide their values for the wavelength interval 400–700 nm. Because the measurements of this study cover a wider range, the effective absorption values were derived for an extended spectral range (300–1000 nm) by simulation.

The calculations were performed using version 1.5 of the radiative transfer model MODTRAN-3 [19], which includes highly resolved spectral absorption coefficients from the HITRAN database [20]. Program settings were as follows. Altitude of surface relative to sea level: 0, observer height: 0, Sun zenith angle: 30°, horizontal visibility: 23 km, rain rate: 0, multiple scattering: yes, model: midlatitude summer, aerosol model: rural, clouds: none. The data interval was set to $40{\mathrm{cm}}^{-1}$ for the range 300–500 nm, $20{\mathrm{cm}}^{-1}$ for 500–725 nm, and $10{\mathrm{cm}}^{-1}$ for 725–1000 nm. Figure 1 shows the calculated transmission spectra for ozone, oxygen, and water vapor (labeled WASI4).

The figure shows for comparison transmission spectra calculated using Eqs. (11)–(13) and the absorption spectra $a_{\mathrm{oz}}(\lambda )$, $a_o(\lambda )$, $a_{\mathrm{wv}}(\lambda )$ of Gregg and Carder [17]. Ozone scale height and water vapor concentration were adjusted to match the MODTRAN spectra.

A new set of absorption spectra $a_{\mathrm{oz}}(\lambda )$, $a_o(\lambda )$, $a_{\mathrm{wv}}(\lambda )$ was calculated for the spectral range 300–1000 nm by inverting the MODTRAN transmission spectra using Eqs. (11–13) and parameter values obtained from matching the transmission spectra: $H_{\mathrm{oz}}=0.38\mathrm{cm}$, $WV=2.5\mathrm{cm}$, and $M_{\mathrm{oz}}=1.07479$, $M^{\prime }=M=1.07834$. A comparison of the new absorption spectra (labeled WASI4) with those of Gregg and Carder [17] is shown in Fig. 2.

The comparison indicates that the new spectra not only cover a wider wavelength range but also are spectrally finer resolved. Gregg and Carder provide no information concerning spectral resolution of their data. To illustrate the effect of spectral sampling, a third data set is included in Fig. 2 (labeled HE5), which is used in the software HydroLight [10]. The three data sets are similar but not identical. The differences concern mainly spectral fine structures, in particular those of $a_{\mathrm{wv}}(\lambda )$. The mismatches are caused by the gases’ numerous absorption lines with bandwidths far below 1 nm. For spectral regions with strong gradients, the irradiance measured by a sensor depends very much on the spectral response of each band, i.e., on center wavelength, width, and spectral shape. The HE5 and WASI4 spectra were calculated using slightly different spectral weighting; hence, their noticeable differences indicate that spectral fine structures of $E_d$ measurements using real sensors should be interpreted with care. Neither the WASI4 nor the HE5 spectra $H_0(\lambda )$, $a_{\mathrm{oz}}(\lambda )$, $a_o(\lambda )$, and $a_{\mathrm{wv}}(\lambda )$ are suited to model properly spectral fine structures of $E_d$ measurements of sensors with a resolution of the order of 1 nm or higher.

4. Optical Properties of the Water Body

A beam of light passing a water layer is affected by absorption and scattering processes along its path, resulting in spectrally dependent intensity changes. For irradiances, a diffuse attenuation coefficient $K$ parameterizes the average changes along the various paths. The bulk coefficient for $E_d$, denoted $K_d$, has been studied extensively (see Bukata [21] for an overview), but I am not aware of analytic models for the coefficients $K_{ds}$ and $K_{dd}$ as defined by Eqs. (2) and (3), respectively.

Because an irradiance sensor detects besides the direct also radiation from angles covering a hemisphere, only a part of the photons that are scattered out of the incident direction is lost for detection. For a beam incident perpendicular on an irradiance sensor, these are the backscattered photons. They are parameterized by the backscattering coefficient $b_b(\lambda )$, which measures the resulting decrease of photon flux per length (in units of ${\mathrm{m}}^{-1}$). Hence, the following approximation is made:

The optical properties of water are calculated as follows:

Four types of water constituents are considered: phytoplankton, gelbstoff, detritus, and suspended particles. The first three are parameterized by their spectral absorption coefficients $a_{\mathrm{ph}}(\lambda )$, $a_Y(\lambda )$, and $a_d(\lambda )$, respectively; suspended particles by their spectral backscattering coefficient, $b_{b,X}(\lambda )$. Backscattering by phytoplankton cells is included in $b_{b,X}(\lambda )$.

Phytoplankton concentration is expressed as mass of the pigments chlorophyll-a plus phaeophytin-a per water volume ($\mathrm{mg}{\mathrm{m}}^{-3}$), its specific absorption coefficient is species dependent. Frequently, a mixture of several phytoplankton species is present in the water. The resulting phytoplankton absorption coefficient is the sum of the individual contributions:

Gelbstoff (yellow substance) is the colored dissolved organic matter (CDOM) in the water and is composed of a huge variety of organic molecules. Its absorption coefficient is calculated as $a_Y(\lambda )={Ya}_Y^{*}(\lambda )$, where $a_Y^{*}(\lambda )$ is the specific absorption spectrum, normalized at ${\lambda }_0$, and $Y=a_Y({\lambda }_0)$ is the absorption coefficient at ${\lambda }_0$. The spectrum $a_Y^{*}(\lambda )$ can either be imported from file, or it can be modeled by the frequently used approximation $\exp [-S(\lambda -{\lambda }_0)]$. For the calculations below, the exponential approximation is used with ${\lambda }_0=440\mathrm{nm}$ and $S=0.014{\mathrm{nm}}^{-1}$, which is representative of a great variety of water types [29,30].

Detritus (also known as tripton [31] or bleached particles [32]) is the collective name for all absorbing nonalgal particles in the water. Its absorption spectrum is parameterized as $a_d(\lambda )={Da}_d^{*}(\lambda )$, with $a_d^{*}(\lambda )$ denoting specific absorption, normalized at the same wavelength ${\lambda }_0$ as Gelbstoff, and $D=a_d({\lambda }_0)$ describing the absorption coefficient at ${\lambda }_0$. The spectrum $a_d^{*}(\lambda )$ is imported from file. In this study, detritus is neglected.

Suspended particle concentration is expressed as mass per water volume ($\mathrm{g}{\mathrm{m}}^{-3}$). Its backscattering coefficient is calculated as

The relative path lengths of diffuse and direct radiation were estimated for $z$ in the range from 0.5 to 5 m and $C_0=2\mathrm{mg}{\mathrm{m}}^{-3}$, $X=0.6\mathrm{g}{\mathrm{m}}^{-3}$, $Y=0.3{\mathrm{m}}^{-1}$ as a function of the Sun zenith angle in water using HydroLight simulations. As described in Section 6, the following relationships were obtained:

5. Software Implementation

The described model was implemented into version 4 of the Water Color Simulator WASI [13–15]. The public-domain software WASI can be used to simulate and analyze the most common types of spectral measurements of shipborne instruments.

Simulation of a measurement (forward modeling) is done by attributing a value to each model parameter and then by calculating the corresponding model curve using a set of equations that represents the model. For simulation of an $E_d(\lambda )$ measurement at depth $z$, the spectra $T_{aa}(\lambda )$, $T_{as}(\lambda )$, $T_r(\lambda )$, $T_{\mathrm{oz}}(\lambda )$, $T_o(\lambda )$, $T_{\mathrm{wv}}(\lambda )$, $E_{dd}(\lambda ,z)$, $E_{ds}(\lambda ,z)$, $E_{dr}(\lambda ,z)$, $E_{da}(\lambda ,z)$, $r_d(\lambda ,z)$, $a(\lambda )$, $K_{dd}(\lambda )$, and $K_{ds}(\lambda )$ are calculated as intermediate results. WASI allows visualization and export to file of each of these spectra. Forward modeling is used in Subsection 6.F to compare the developed irradiance model with a reference model.

Data analysis (inverse modeling) aims to determine the values of unknown model parameters, called fit parameters. The fit parameters of the irradiance model are listed in Table 1. Their values are determined iteratively as follows. In the first iteration, a model spectrum $E_d(\lambda )$ is calculated using Eq. (1) for user-defined initial values of the fit parameters. This spectrum is compared with the measured one, $E_d^{\text{meas}}(\lambda )$, by calculating the residuum $\sum {|E_d^{\text{meas}}({\lambda }_i)-E_d({\lambda }_i)|}^2$ as a measure of correspondence. Then, in the further iterations, the fit parameter values are altered using the Downhill Simplex algorithm [33,34], resulting in altered model curves and altered residuals. The procedure is stopped when the calculated and the measured spectrum agree as good as possible, which corresponds to the minimum residuum. The values of the fit parameters of the final step are the fit results. WASI provides different methods to tune this algorithm, in particular, wavelength-dependent weighting and changing the residuum definition; for details, see Gege and Albert [14]. Inverse modeling is used in Section 7 to analyze field measurements.

Table 1. Parameters of the Irradiance Model That Can Be Used in WASI as Fit Parameters

View Table | View all tables in this article

If the initial value of a fit parameter is very different from its correct value, the inversion algorithm may not be able to find the parameter combination for which the model curve has the best correspondence with the measured spectrum. This problem, which increases with the number of fit parameters, can be reduced if a reasonable “first guess” of the fit parameters can be made, in particular, of those parameters that can cause large deviations. Sensor depth is such a parameter. An analytic algorithm has been implemented into WASI, which estimates a first guess, $z_0$, from the ratio $r=E_d({\lambda }_1)/E_d({\lambda }_2)$ of an irradiance measurement at two wavelengths, ${\lambda }_1$ and ${\lambda }_2$, using the following equation:

6. Comparison with a Reference Model

The commercial software HydroLight [10] is taken as reference to determine some model parameters ($l_{dd}$, $l_{ds}$, ${\rho }_{ds}$) and to estimate the accuracy of the analytical model. HydroLight is probably the most widely used numerical model for calculating the radiative transfer in water bodies. The current version is called HE5 (HydroLight-EcoLight version 5.1). It solves the radiative transfer equation in scattering media using the invariant embedding method [35].

A. Consistency of Input Spectra

A comparison of the extraterrestrial solar irradiance spectrum, $H_0(\lambda )$, is shown in Fig. 4. The HE5 and WASI4 spectra of $H_0(\lambda )$ correspond well. Typical differences are of the order of a few percent, except concerning spectral fine structures. Similar results were obtained in Section 3 (see Fig. 2) for the absorption spectra of the atmospheric gases ozone, oxygen, and water vapor. Consequently, all input spectra of WASI for calculating irradiance above the surface can be considered consistent with those of HE5. Remember that spectral fine structures of the order of 1 nm cannot be modeled accurately. To avoid discrepancies resulting from slightly different input spectra, the RADTRAN-X database file of HE5 was exchanged for all further runs by the spectra $H_0(\lambda )$, $a_{\mathrm{oz}}(\lambda )$, $a_o(\lambda )$, and $a_{\mathrm{wv}}(\lambda )$ used in WASI. Also, the optical properties of the water body were made consistent by using the spectra $a_W(\lambda )$ and $a_{\mathrm{pi}}^{*}(\lambda )$ of the WASI database in HE5.

 

Fig. 4. Extraterrestrial solar irradiance.

Download Full Size | PDF

B. Consistency of Above-Water Irradiance Model

Irradiance spectra $E_d(\lambda )$, $E_{dd}(\lambda )$, $E_{ds}(\lambda )$ were compared above the water surface for the spectral range from 350 to 900 nm. HE5 calculates these spectra using RADTRAN-X, which is an implementation of the Gregg and Carder model [17] using an extended database covering the range 300 to 1000 nm. The bands are specified in HE5 by their border wavelengths, i.e., irradiance $E$ of band $k$ is calculated as the average $〈E({\lambda }_i)〉$, where ${\lambda }_i$ runs from the center wavelength of band $k-1$ to the center wavelength of band $k+1$ in steps of 1 nm. For the calculations described below, such spectral averaging was also implemented in WASI4, and a sampling interval of 2 nm was chosen; i.e., irradiance at wavelength $\lambda$ was calculated as $[E(\lambda -1\mathrm{nm})+E(\lambda )+E(\lambda +1\mathrm{nm})]/3.$

The common parameters of WASI4 and HE5 were set to $d=94$ (April 4), $\mathrm{AM}=1$, $\mathrm{RH}=60\%$, $P=1013.25\mathrm{mbar}$, $WV=2.5\mathrm{cm}$, and $H_{\mathrm{oz}}=0.3\mathrm{cm}$. Horizontal visibility was chosen as $V=15\mathrm{km}$ in HE5. The aerosol parameters were set to $\alpha =1.317$, $\beta =0.2606$, and $H_a=1\mathrm{km}$ in WASI4. Because these cannot be selected directly in HE5, their values were obtained as follows. The sourcecode of HE5 was checked to obtain $H_a=1\mathrm{km}$ and to confirm $\beta =3.91H_a/V$. $\alpha$ could not be derived from the sourcecode because it is calculated in a complex way using a marine aerosol model consisting of three components. Hence, $\alpha$ and $\beta$ were tuned in WASI4 until the $E_{ds}$ spectra for ${\theta }_{\text{sun}}=30°$ corresponded to the HE5 spectra. An example for the correspondence of WASI4 and HE5 irradiance spectra is shown in Fig. 5.

 

Fig. 5. Comparison of above water irradiance spectra for ${\theta }_{\text{sun}}=30°$.

Download Full Size | PDF

The correspondence is excellent: differences between WASI4 and HE5 irradiance spectra are less than line thickness for the majority of spectral bands. For this reason only the HE5 spectra are displayed. The ratio plot shows that the differences are always below 0.1%. Correspondence is of similar quality to Fig. 5 also for other Sun zenith angles ranging from 0° to 80° (data not shown), i.e., HE5 and WASI4 provide consistent irradiance spectra $E_d(\lambda )$, $E_{dd}(\lambda )$, and $E_{ds}(\lambda )$above water, indicating that the Gregg and Carder model [17] is implemented identically in both software.

C. Parameterization of ${\rho }_{ds}$

The reflectance of the diffuse component of irradiance at the water surface, ${\rho }_{ds}$, depends on the Sun zenith angle. To obtain a parameterization of ${\rho }_{ds}({\theta }_{\text{sun}})$, pairs of irradiance spectra $E_{ds}(\lambda ,0+)$, $E_{ds}(\lambda ,0-)$ were calculated for Sun zenith angles ranging from 0° to 80°. Because $E_{ds}(\lambda ,0-)$ is no output of HE5, $E_{dd}(\lambda ,0+)=0$ was set in the sourcecode of EcoLight. The calculated irradiance is then just the diffuse component. Furthermore, that part of $E_{ds}(\lambda ,0-)$ that is caused by reflection of upwelling radiation at the water surface has been made zero by setting the scattering coefficients of water and all water constituents to zero. ${\rho }_{ds}$ was then calculated as the average of ${\rho }_{ds}(\lambda )=1-E_{ds}(\lambda ,0-)/E_{ds}(\lambda ,0+)$ for the wavelength range from 400 to 700 nm. The result is shown in Fig. 6.

 

Fig. 6. Dependency of the reflectance factor for diffuse irradiance on the Sun zenith angle.

Download Full Size | PDF

The obtained ${\rho }_{ds}$ values can be approximated by Eq. (5), which was derived by fitting the values shown in Fig. 6 with a second-order polynomial function.

D. Parameterization of $l_{ds}$

A factor $l_{ds}$ was introduced in Eq. (2) as the average path length of diffuse radiation relative to sensor depth $z$. To obtain numerical values for $l_{ds}$, a set of 130 spectra $E_{ds}(\lambda ,z,{\theta }_{\text{sun}})$ was calculated using HE5 for $z$ ranging from 0.5 to 5 m and ${\theta }_{\text{sun}}$ values between 0° and 80°. The concentrations of water constituents were set to $C_0=2\mathrm{mg}{\mathrm{m}}^{-3}$, $X=0.6\mathrm{g}{\mathrm{m}}^{-3}$, and $Y=0.3{\mathrm{m}}^{-1}$. As before, $E_{dd}(\lambda ,0+)=0$ was set in the sourcecode of EcoLight to get just the diffuse component. $l_{ds}$ was then calculated as

 

Fig. 7. Dependency of the relative path length of diffuse radiation on depth and on the Sun zenith angle.

Download Full Size | PDF

E. Parameterization of $l_{dd}$

Similarly, as before, relative path length factors for direct radiation, $l_{dd}$, were obtained by solving Eq. (3) for $l_{dd}$ and generating a number of spectra $E_{dd}(\lambda ,z,{\theta }_{\text{sun}})$ using HE5. $E_{ds}(\lambda ,0+)=0$ was set in the sourcecode to force HE5 to calculate $E_{dd}$ instead of $E_d$. The simulations were performed using a modified version of EcoLight, which has an angular resolution of 2° instead of the standard 10°. In this way, a set of spectra $E_{dd}$ was calculated for the same wavelengths, Sun zenith angles, and depths as above for $E_{ds}$. The derived values of $l_{dd}$ have an average of 1.000 and a standard deviation of 0.004. Hence, the dependency on $\lambda$, $z$, and ${\theta }_{\text{sun}}$ can be neglected for the studied conditions, and $l_{dd}=1$ is set.

F. Validation

The developed model was validated by comparing irradiance spectra obtained from WASI4 using forward modeling with the corresponding spectra from HE5 for 99 different combinations of $z$ and ${\theta }_{\text{sun}}$. The ranges were 0 to 5 m for $z$ and 0° to 80° for ${\theta }_{\text{sun}}$. The atmosphere parameters were chosen as in Section 6.B. The water parameters were as follows: $C_0=2\mathrm{mg}{\mathrm{m}}^{-3}$, $X=0.6\mathrm{g}{\mathrm{m}}^{-3}$, $Y=0.3{\mathrm{m}}^{-1}$, and $n_W=1.33$. For calculating $E_{dd}(\lambda )$ using HE5, $E_{ds}(\lambda ,0+)=0$ was set in the HE5 sourcecode, and $E_{dd}(\lambda ,0+)=0$ was set to obtain $E_{ds}(\lambda )$. The computation time for a single spectrum with 276 spectral bands was on average 120 s for the 2° version of EcoLight and ${10}^{-3}\mathrm{s}$ for WASI4 on a 2 GHz Intel Core 2 Duo processor, i.e., the analytical model is approximately ${10}^5$ times faster. Recently a fast version of EcoLight has been developed, called EcoLight-S, with typical computation times of less than 1 s [36].

A typical matchup is shown in Fig. 8 for $z=4\mathrm{m}$ and ${\theta }_{\text{sun}}=40°$. The HE5 and WASI4 curves are hardly distinguishable by eye. They differ less than $\pm 1.5\%$ in the range from 400 to 700 nm; the average ratio for that range is 0.9991 for $E_{dd}$, 1.0097 for $E_{ds}$, and 1.0057 for $E_d$.

 

Fig. 8. Comparison of spectral shape for WASI4 and HE5 simulations of irradiance.

Download Full Size | PDF

Figure 9 shows the spectrally averaged ratios for all matchups of $E_d$, $E_{dd}$, and $E_{ds}$. Averaging was performed from 400 to 700 nm. In most cases, the deviations are below $\pm 1\%$. Deviations of more than 1.5% occur only for Sun zenith angles above 65° and depths above 3 m. The average ratio is $1.002\pm 0.003$ for $E_d$, $0.998\pm 0.005$ for $E_{dd}$, and $1.006\pm 0.003$ for $E_{ds}$ for $z$ ranging from 0.5 to 5 m and ${\theta }_{\text{sun}}$ from 0° to 70°. Hence, WASI4 and HE5 provide consistent spectra $E_d(\lambda )$, $E_{dd}(\lambda )$, and $E_{ds}(\lambda )$ for a plane water surface under the studied conditions.

 

Fig. 9. Comparison of intensity for WASI4 and HE5 simulations of irradiance.

Download Full Size | PDF

7. Application to Field Measurements

Measurements were performed in 2003 and 2004 in the German lakes Bodensee (Lake Constance), Starnberger See, and Waginger See mostly in shallow waters using a small boat [37]. A data set of 421 $E_d$ measurements was collected using a RAMSES-ACC-VIS irradiance sensor (TriOS, Oldenburg, Germany). Each measurement consisted of a sequence of 4 to 50 irradiance spectra, which cover the range from 350 to 900 nm at a spectral sampling interval of 5 nm. 98% of the data have integration times between 16 and 64 ms, and the average is 34 ms. Measurement time of a sequence varied from 21 to 700 s with an average of 105 s. The $E_d$ sensor was lowered into the water at the sunlit side of the boat at a distance of 2 to 3 m to avoid shadowing. More details can be found in Gege and Pinnel [5].

A representative example of a single spectrum and the corresponding fit curve obtained by inverse modeling using WASI is shown in Fig. 10(a). The measurement was performed on a cloudless day at Bodensee (July 28, 2004, 15:12 h local time). For illustration purposes, the model curve was calculated at higher spectral resolution (1 nm) than the measurement (5 nm). The plot demonstrates that the downwelling irradiance reveals many spectral fine structures that are not resolved by the instrument; thus, the measured spectrum is smooth compared with the model curve. Except this difference of spectral resolution, the calculated spectrum fits well to the measurement. The obtained fit parameters are listed in the figure’s legend. A further result of inverse modeling using WASI is the separation of the direct and diffuse components of irradiance, $E_{dd}(\lambda )$ and $E_{ds}(\lambda )$; see Fig. 10(b). Their spectral shapes differ significantly. Consequently, changes of their relative intensities alter the spectral shape of $E_d$.

 

Fig. 10. Example for inverse modeling of an irradiance measurement. Fit parameters: $X=0.85\mathrm{g}{\mathrm{m}}^{-3}$, $Y=0.21{\mathrm{m}}^{-1}$, $f_{dd}=1.02$, $f_{ds}=1.02$, $z=1.06\mathrm{m}$. (a) Measured spectrum and model curve. (b) Direct and diffuse component obtained from inverse modeling.

Download Full Size | PDF

Inverse modeling was performed for all 4375 individual $E_d$ spectra of the 421 measurements. The model curves were calculated for the spectral range from 400 to 800 nm in 5 nm steps using $X$, $Y$, $f_{dd}$, $f_{ds}$, and $z$ as fit parameters. For ${\theta }_{\text{sun}}$ the actual angle at the time of the measurement was taken, and $C_0=2\mathrm{mg}{\mathrm{m}}^{-3}$ was set, which is between the average concentration of Starnberger See ($1.4\pm 0.3\mathrm{mg}{\mathrm{m}}^{-3}$) and Bodensee ($2.5\pm 0.9\mathrm{mg}{\mathrm{m}}^{-3}$) as derived from simultaneous in situ measurements [38].

The results for $f_{dd}$ and $f_{ds}$ are shown in Fig. 11. These parameters describe the variability of $E_d$ induced by the changing geometry of a wind-roughened and wave-modulated water surface. The histograms of both parameters peak near the expected value for undisturbed geometry, $f_{dd}=1$ and $f_{ds}=1$, suggesting that Eq. (1) is a reasonable model to separate the direct and diffuse components of $E_d$. The asymmetry of the $f_{dd}$ histogram, which favors values $<1(f_{dd}=0.88\pm 0.38)$, is caused by the numerous measurements that were performed under cloud cover. The $f_{ds}$ histogram is symmetrical around $f_{ds}=1.04$ with a standard deviation of 0.46. The $f_{dd}$ and $f_{ds}$ values shown here were analyzed in more detail in Gege and Pinnel [5] to determine quantitatively the depth dependent influence of wave focusing on $E_d$ during a measurement. It was found that $f_{dd}$ and $f_{ds}$ are only weakly correlated ($r^2$ is between 0.02 and 0.16), changes of $f_{dd}$ explain up to 82% of intensity variability and alter the spectral shape of $E_d$ between 400 and 700 nm on average by 5.4%, and changes of $f_{ds}$ have only minor impact on $E_d$.

 

Fig. 11. Frequency distributions of the parameters $f_{dd}$ and $f_{ds}$.

Download Full Size | PDF

Because the spectral shape of $E_d$ depends on the ratio $r_d$ of direct to diffuse irradiance [5], the variability of $r_d$ during a measurement is of interest to estimate wavelength-dependent changes of $E_d$. The $r_d$ values of all 4375 spectra are shown in Fig. 12 as a function of the Sun zenith angle. They were calculated as spectral average of $f_{dd}{E}_{dd}(\lambda )/f_{ds}{E}_{ds}(\lambda )$ in the range from 400 to 700 nm. The solid red line shows for comparison the function $r_d({\theta }_{\text{Sun}})$ as expected for $f_{dd}=f_{ds}=1$; it was calculated using Eq. (17) for the average depth of 0.7 m of all measurements. Obviously, the large and uncorrelated variability of $f_{dd}$ and $f_{ds}$ induces very large variability to $r_d$. The standard deviation of $r_d$ is 4.5 for the actual data set, which is far above the average of 2.6.

 

Fig. 12. Variability of the ratio of direct to diffuse irradiance.

Download Full Size | PDF

Sensor depth $z$ changes during a measurement because of waves that alter the thickness of the water column above the sensor, and because of roll movements of the boat. Figure 13 shows a statistics of these changes. The $z$ values obtained by inverse modeling of the individual spectra forming a measurement were averaged, and then the differences $\mathrm{\Delta }z$ between the individual $z$ values and the mean were calculated. The standard deviation of $\mathrm{\Delta }z$ is 0.043 m. Because a portion of this variability is caused by wave action and boat movement, it can be concluded that the uncertainty of $z$ determination by inverse modeling is below 4 cm.

 

Fig. 13. Sensor depth variability during a measurement.

Download Full Size | PDF

The accuracy of water constituents determination increases with the thickness of the water column above the irradiance sensor. An analysis of the Bodensee and Starnberger See data set showed that the uncertainty of phytoplankton concentration decreases exponentially with depth and reaches class-specific lower limits at depths between 1 and 1.5 m [38]. For this reason, only measurements from depths $z>1.5\mathrm{m}$ were used to calculate lake-specific averages of the parameters $X$ and $Y$. These are summarized in Table 2. The standard deviations reflect the natural concentration variability together with the uncertainties of inverse modeling. The averages for all three lakes, $X=0.6\mathrm{g}{\mathrm{m}}^{-3}$ and $Y=0.3{\mathrm{m}}^{-1}$, were used in Gege and Pinnel [5] as typical concentrations to study the impact of environmental and experimental conditions on the variance of downwelling irradiance in water.

Table 2. Mean and Standard Deviation of Suspended Matter Concentration ($X$, ${\sigma }_X$) and Gelbstoff Concentration ($Y$, ${\sigma }_Y$) Derived by Inverse Modeling of Measurements at Depths$>1.5\mathrm{m}$

View Table | View all tables in this article

Because the in situ measurements indicated only little variability of phytoplankton concentration, the parameter $C_0$ was set to $\text{2\hspace{0.17em}}\mathrm{mg}{\mathrm{m}}^{-3}$ and kept constant during inverse modeling in this study. It was shown in Gege [38] that the described model also can be used to determine the concentration of phytoplankton. By applying an adapted fit strategy to the same data set from Bodensee and Starnberger See, it was possible to determine the concentrations of three phytoplankton classes (diatoms, dinoflagellates, green algae) above class-specific thresholds between 0.4 and $0.9\mathrm{mg}{\mathrm{m}}^{-3}$. The uncertainty of total pigment concentration (sum of chlorophyll-a and phaeophytin-a in all classes) was $0.7\mathrm{mg}{\mathrm{m}}^{-3}$.

8. Summary

An analytic model of the downwelling irradiance in water ($E_d$) was developed that calculates the direct ($E_{dd}$) and diffuse ($E_{ds}$) components separately. This separation allows to account for the different path lengths of the two components and to handle the large variability of $E_d$ at typical field conditions. Because the computation time is of the order of ${10}^{-3}\mathrm{s}$, the model can be used for computationally extensive simulations like inverse modeling.

The intensities and wavelength dependencies of $E_{dd}$ and $E_{ds}$ just beneath the water surface are parameterized using the model of Gregg and Carder for cloudless maritime atmospheres. Their database, which is restricted to the spectral range of 400–700 nm, was extended to a range of 300–1000 nm through radiative transfer simulations using MODTRAN. Changes of $E_{dd}$ and $E_{ds}$ with depth ($z$) are calculated individually to account for the different path lengths of direct ($z/\cos {\theta }_{\text{Sun}}^{\prime }$) and diffuse (${zl}_{ds}$) radiation. The radiative transfer program HydroLight was used to show that using these path lengths, the Lambert—Beer law describes accurately the depth dependency of both irradiance components with a common attenuation coefficient, $a+b_b$, where $a$ is the average absorption coefficient and $b_b$ is the average backscattering coefficient of the water column between surface and depth $z$. The relative path length of diffuse irradiance ($l_{ds}$) was parameterized as a function of the Sun zenith angle in water (${\theta }_{\text{sun}}^{\prime }$) for concentrations of water constituents that are typical for the German lakes Bodensee, Starnberger See, and Waginger See.