2.1.

Radiative Transfer Simulation

The derivation of model (1) by Lyzenga³ presumes a level surface, homogeneous optical properties of the water, and a level Lambertian bottom with homogeneous reflectance. It also assumes that the effects of polarization, fluorescence, and inelastic scattering are negligible. We simulated radiative transfer in order to estimate the “true” value of $R$ in light of the radiative transfer theory under these assumptions: our simulations are based on the same assumptions as model (1) except for the ignorance of internal reflection and the simplification of multiple scattering. This enables us to examine only the inaccuracy caused by the ignorance or simplification.

Under these assumptions, $R$ depends only on the optical conditions listed in Table 1. Specifically, these conditions are the incident angle ${\theta }_i$ of the incident beam at the surface; the zenith angle ${\theta }_L$ of $L$; the single scattering albedo, defined as $b/c$ ($b$: the scattering coefficient; $c$: the beam attenuation coefficient, defined as $c\equiv a+b$, where $a$ is the absorption coefficient); the scattering phase function $\tilde{\beta }$; the optical depth $c·h$ of water; the irradiance reflectance $r_b$ of the bottom; and the relative refractive index of the air–water interface $n_{aw}$.

Table 1

Optical conditions and values set in radiative transfer simulation.

We calculated $R$ on the basis of the radiative transfer simulation for all the combinations of optical conditions shown in Table 1 ($4\times 1\times 3\times 3\times 15\times 11\times 1=5940$ combinations). ${\theta }_i$ values of 0, 30, and 60 deg represent direct sunlight with different solar elevations, and ${\theta }_i=\text{uniform}$ (the incident radiance is uniformly distributed over the sky) approximates diffuse sky light. ${\theta }_L$ was fixed at zero to simulate observation by quasi-zenith satellites. Three values of $b/c$ (0.2, 0.5, and 0.8) and three types of $\tilde{\beta }$ (“Mob,” “M01,” and “P02”) were used to represent various in-situ measurements.^19–21 Here, “Mob” indicates the scattering phase function described in Refs. 17 and 18 obtained by averaging the measurements of Ref. 19. The probability of backscattering $B$ (the ratio of the backward scattering coefficient to $b$) was 0.0178 in our implementation. “M01” and “P02” are the scattering phase functions presented in Ref. 20 and named in Ref. 21, characterized by small (0.00453) and large (0.0445) $B$, respectively. We set $c·h$ in geometric progression from ${10}^{-2}(=0.01)$ to ${10}^{0.6}(=3.981)$ with a ratio of ${10}^{0.2}(=1.585)$. This is because the scale of $h$ of interest in shallow-water remote sensing using multispectral satellite imagery is diverse: sometimes very shallow ranges such as 0 to 0.3 m¹³ are discussed in a centimeter scale; in other cases, wide ranges of about 1 to 20 m^1,22 are targeted. The $c·h$ range handled in this paper corresponds to 0.05 to 19.91 m when $c=0.2$ (a possible value for clear seawater at blue and green wavelengths).^18,19,23 We set $r_b$ to cover a wide range of 0.1 to 0.6 with a small interval of 0.05. An $r_b$ value of 0.6 is possible for carbonate sand at green and red wavelengths.²⁴ Although an $r_b$ of $<0.1$, such as 0.05, is common for algae and corals,²⁴ this condition was not considered because it is not suitable for the use of model (1) in the form of Eq. (2): $R-R_{\infty }$ becomes nonpositive when $R_b\le R_{\infty }$, and we cannot calculate the logarithm in Eq. (2). Because $n_{aw}$ is rather stable in natural waters,¹⁸ it was fixed at a typical value (1.333).

Our simulation was based on the forward Monte Carlo method described in Ref. 25. For each combination of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $c·h$, we injected ${10}^{10}$ photons downward to the water surface and traced their coordinates until they disappeared by absorption or flew into the air. In this process, we used MT19937²⁶ pseudo-random numbers with the optical conditions to determine the behaviors of each photon when it interacts with the water surface, water body, and bottom: whether the photon is reflected or refracted when it hits the water surface, the distance it travels before interacting with water, the type of interaction (scattering or absorption), and the direction it travels after being scattered in water or reflected at the bottom. Then, we tallied the photons just above the surface to estimate $R$. Here, because we fixed ${\theta }_L$ at zero, only the upwelling photons with zenith angles smaller than 5 deg were counted to calculate $L$.

The simulation code was developed in the Fortran 90 language and validated using the approaches described in Refs. 27 and 28. The former approach is based on a rigorous relationship among the vector and scalar irradiances and $b/c$ (the first formula of Ref. 27). The latter approach is based on an analytical solution (shown in Table 1 of Ref. 28) for a collimated beam normally incident on a homogeneous single-layer slab with isotropic scattering. For further validation, we also developed a backward Monte Carlo code and confirmed that its output was consistent with that of the forward Monte Carlo code.

Gordon and Brown²⁵ showed that any radiometric quantity $Q$ for an arbitrary bottom reflectance ($r_b$) can be exactly computed from the contribution to $Q$ from the photons that strike the bottom zero, one, and two times when $r_b=1$. Therefore, we performed the simulation with $r_b=1$ for each combination of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $c·h$, and the $R$ values for various $r_b$ listed in Table 1 were then calculated using this time-saving technique. In addition to $R$, the remote-sensing reflectance without the contribution of photons internally reflected at least once ($R^{\prime }$) was also tallied.

The upward reflection at the water surface was not simulated directly except for the loss of incident photons in the case of ${\theta }_i=\text{uniform}$. Instead, the effect was corrected afterward based on Fresnel equation.

2.2.

Accuracy Assessment

First, we evaluated how accurately model (1) can describe the relationship between $c·h$ and $R$ calculated by the radiative transfer simulation for each combination of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$. Model (1) can be transformed to the following form via Eq. (2):

For each combination, we prepared a dataset consisting of the $c·h$ values listed in Table 1 (except for $\infty$; 14 in total) and the corresponding values of the simulated $\log [R-R_{\infty }]$. Here, the simulated $R$ for $c·h=\infty$ was used as $R_{\infty }$. Next, Eq. (3) was fitted to the dataset by least squares fitting, treating $\log [R_b-R_{\infty }]$ and $k/c$ as free parameters. Then, the determination coefficient ($R^2$) and root mean square residual (RMSR) of the fitting were evaluated. We also performed the same evaluation for $R^{\prime }$ (by replacing $R$ and $R_{\infty }$ with $R^{\prime }$ and $R_{\infty }^{\prime }$ in the above procedure) in order to examine the effect of internal reflection on the model’s accuracy.

In most of the applications of model (1) where $c$ is assumed to be a constant, the relationship between $R$ (or $R^{\prime }$) and $c·h$ is essentially equivalent to that between $R$ (or $R^{\prime }$) and $h$. Therefore, the above procedure is equivalent to an accuracy evaluation of model (1) in terms of the relationship between $R$ (or $R^{\prime }$) and $h$ for arbitrary $c$.

Second, we investigated the error caused by using model (1) in remote sensing of water depth when the depth range of the calibration data is limited. We considered the simplest case in which the bottom type is uniform and thus only one visible band is used. First, Eq. (3) was calibrated with (fitted by least squares to) the data for limited $c·h$ ranges selected from the dataset described above. Then, the calibrated Eq. (3) was used to estimate $c·h$ from the simulated $\log [R-R_{\infty }]$ values for the entire $c·h$ range, and the estimation error was evaluated.

3.1.

Accuracy of Model (1)

Figure 2 demonstrates the relationships between $R$ or $R^{\prime }$ and $c·h$ calculated by the radiative transfer simulation for several combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$. Obviously, $R$ or $R^{\prime }$ increases with decreasing $c·h$ and increasing $r_b$. Naturally, this is due to the increase in the bottom reflection component of $R$ or $R^{\prime }$.

Fig. 2

Relationship between $R$ or $R^{\prime }$ and $c·h$ calculated by radiative transfer simulation for several combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$. The curves of model (1) fitted [in the form of Eq. (3)] to the plotted data and their $R^2$ values are also shown.

Figure 2 also shows the curves of model (1) as expressed in Eq. (3) fitted to the plotted data as described above. Overall, the fitting is fairly good for both $R$ and $R^{\prime }$. The model appears to underestimate small $c·h$ values when $r_b$ is large, but the error (residual of the fitting) is actually small: the error in the region is exaggerated because the $c·h$ axis is in log scale.

Table 2 lists the statistics of the $R^2$ and RMSR values of the fitting for all the combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$. According to this table, the overall mean $R^2$ is 0.9991 for $R$ and 0.9997 for $R^{\prime }$. The overall mean RMSR is 0.02861 for $R$ and 0.01274 for $R^{\prime }$. A RMSR of 0.02861 is just 0.72% of the target $c·h$ range (${10}^{0.6}-{10}^{-2}=3.971$). The overall minimum $R^2$ for $R$ is 0.9935, whereas the overall maximum RMSR for $R$ is 0.09158. Even this maximum RMSR is just 2.3% of the target $c·h$ range. These results indicate that when model (1) is fitted to data that cover the entire target range of $c·h$, it is reasonably accurate in describing the relationship between $R$ and $c·h$ although the internal reflection slightly degrades the accuracy on average. Figure 3 shows $R-R^{\prime }$ (the contribution of internally reflected photons to $R$) as a function of $c·h$ for the same combinations of optical conditions as in Fig. 2. We can observe that $R-R^{\prime }$ increases with decreasing $c·h$ and increasing $r_b$. In fact, $R-R^{\prime }$ increased monotonically with decreasing $c·h$ for all the combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$ ($4\times 3\times 3\times 11=396$ combinations) except for two combinations with $r_b=0.1$, in which $R-R^{\prime }$ for $c·h={10}^{0.6}$ was slightly larger than that for ${10}^{0.4}$. $R-R^{\prime }$ increased monotonically with increasing $r_b$ for all the combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $c·h$ ($4\times 3\times 3\times 14=504$ combinations) without exception.

Table 2

Statistics of R2 and RMSR for fitting model (1) [in the form of Eq. (3)] to the dataset of c·h and simulated log[R−R∞] for all the combinations of θi, b/c, β˜, and rb. The values for R′ (results obtained when log[R′−R∞′] is used instead of log[R−R∞]) are also shown.

Fig. 3

$R-R^{\prime }$ (contribution of internally reflected photons to $R$) as a function of $c·h$ for the same combinations of optical conditions as in Fig. 2. The curve of exponential function (4) fitted to the plotted data is also shown.

Figure 4 is a contour plot of the averaged ratio of $R-R^{\prime }$ to $R$ versus $r_b$ and $c·h$. We can see from this figure that internal reflection significantly increases $R$ even in common conditions. For example, the averaged ratio of $R-R^{\prime }$ to $R$ exceeds 0.05 (5%) when $c·h={10}^{-0.2}\cong 0.63$ and $r_b=0.4$, as indicated in the figure by the symbol “$\times$.” Here, a $c·h$ value of 0.63 corresponds to an $h$ value of 2.1 m when $c=0.3$, and both of these values are not too small to be common: $c<0.3$ has been observed in various locations in oceans and the Mediterranean for green wavelengths^19,23 and in a coral reef even at 660 nm.²⁹ An $r_b$ value of 0.4 is not a large value for carbonate sand in coral reefs: the average $r_b$ of 670 measurements by Hochberg²⁴ falls within 0.4 to 0.6 at wavelengths of 475 to 700 nm. The maximum value of the ratio of $R-R^{\prime }$ to $R$ for all the combinations of optical conditions listed in Table 1 was 0.28 (28%).

Fig. 4

Contour plot of averaged ratio of $R-R^{\prime }$ to $R$ against $r_b$ and $c·h$. For each combination of $r_b$ and $c·h$, there are $4\times 3\times 3=36$ combinations of ${\theta }_i$, $b/c$, and $\tilde{\beta }$. Each plotted value is the average of the ratio of $R-R^{\prime }$ to $R$ for the 36 combinations. The symbol “$\times$” indicates the point for $r_b$ =0.4 and $c·h={10}^{-0.2}\cong 0.63$.

Figure 3 also shows the curve of the following exponential function fitted to the plotted data:

3.2.

Error Caused by Model (1) When Calibration Data are Limited

Figure 5 shows the relationship between $R$ and $c·h$ calculated by the radiative transfer simulation for the same combinations of optical conditions as in Fig. 2. However, the curves in Fig. 5 are not those of model (1) calibrated using all the data plotted. The red curve shows the model as expressed in Eq. (3) fitted to only the five data points with the smallest $c·h$ ($c·h={10}^{-2},{10}^{-1.8},{10}^{-1.6},{10}^{-1.4},{10}^{-1.2}$), and the blue curve shows model (1) fitted to only the five data points with the largest $c·h$ ($c·h={10}^{-0.2},{10}^0,{10}^{0.2},{10}^{0.4},{10}^{0.6}$). This figure shows what happens when the depth range of the calibration data does not fully cover the target depth range for depth estimation using model (1). For medium and large $r_b$ ($r_b=0.3$, 0.5), the red curve underestimates $c·h$ in the region with large $c·h$, whereas the blue curve underestimates $c·h$ in the region with small $c·h$. This shows that when model (1) is used to estimate $h$ from $R$, it significantly underestimates $h$ if it is calibrated using data for a different depth range. In short, model (1) cannot accurately extrapolate $h$ using $R$.

Fig. 5

Relationship between $R$ and $c·h$ calculated by radiative transfer simulation for the same combinations of optical conditions as in Fig. 2. The curves shown are model (1) [in the form of Eq. (3)] fitted to the five plotted data with largest $c·h$ (blue curve) and smallest $c·h$ (red curve). The maximum error (m.e.) of estimating the plotted $c·h$ from $R$ using each curve is also shown. The m.e. for the red curve is the error for the maximum $c·h({10}^{0.6})$, and the m.e. for the blue curve is that for the minimum $c·h({10}^{-2})$.

Figure 5 also shows the maximum error (“m.e.” in the figure) of estimating the plotted $c·h$ from $R$ using each curve. In other words, the maximum error is the maximum difference in $c·h$ between the curve and the plotted data points. Naturally, for the red curve, which was calibrated using the data points with the smallest $c·h$, the maximum error is the error for the data point with the largest $c·h({10}^{0.6})$. Similarly, for the blue curve, which was calibrated using the data points with the largest $c·h$, the maximum error is the error for the data point with the smallest $c·h({10}^{-2})$. For example, when ${\theta }_i=30$, $b/c=0.2$, $\tilde{\beta }$=“Mob,” and $r_b=0.5$, the maximum error for the red curve is 1.382. This means that the red curve underestimated the maximum $c·h({10}^{0.6}=3.981)$ by 1.382. When $c=0.3$, this error is as large as $1.382/0.3=4.61\mathrm{m}$, which is 35% of the target depth range (13.24 m).

Table 3 lists the top three maximum errors for all the combinations of ${\theta }_i$, $b/c$, $\tilde{\beta }$, and $r_b$. The largest error, 1.705, was observed when ${\theta }_i=0$, $b/c=0.5$, $\tilde{\beta }$=“M01,” and $r_b=0.6$. This error corresponds to $1.705/0.3=5.68\mathrm{m}$ of depth estimation error when $c=0.3$ and accounts for as much as 43% of the target depth range (13.24 m).

Table 3

The top three maximum errors in c·h estimation using model (1) [in the form of Eq. (3)] fitted to only five data points with smallest or largest c·h. The corresponding optical conditions and the maximum error for R′ (the maximum error obtained when log[R′−R∞′] is used instead of log[R−R∞]) are also shown.

Table 3 also shows the maximum error for $R^{\prime }$ (the maximum error obtained when $\log [R^{\prime }-R_{\infty }^{\prime }]$ is used instead of $\log [R-R_{\infty }]$). Considering that the maximum errors for $R^{\prime }$ are $<10\%$ of those for $R$, the large errors for $R$ can mostly be attributed to the fact that model (1) ignores internal reflection. Although both $R^{\prime }$ and $R-R^{\prime }$ decrease exponentially with $c·h$ as discussed above, their sum $R$ is only approximate and not strictly an exponential function of $c·h$, and hence the exponential model (1) cannot be used for extrapolation.

3.3.

Mathematical Presentation of Simulation Results

To enable readers to reproduce the simulation results obtained in this study for arbitrary conditions, we empirically modeled $R^{\prime }$ and $R-R^{\prime }$ as follows:

Here, $R_s$ is the contribution of specular reflection of incident light at the surface, which is nonzero only for ${\theta }_i=0$ and ${\theta }_i=\text{uniform}$. Based on Fresnel equation, $R_s$ is $0.02037\pi$ for ${\theta }_i=0$ and 0.02037 for ${\theta }_i=\text{uniform}$. $R_{w\infty }^{\prime }$ is the contribution of in-water scattering for infinitely deep water, which was calculated on the basis of Ref. 30 as

The decimal coefficients in models (5) and (6) were determined by least squares fitting of the models to the simulated $R^{\prime }$ and $R-R^{\prime }$ for all the combinations of optical conditions shown in Table 1. The $R^2$ values of the fitting of models (5) and (6) were as large as 0.999854 and 0.999152, respectively. The RMSR values were as small as 0.001132 and 0.000895, respectively. These results indicate that the models can accurately reproduce the simulated $R^{\prime }$ and $R-R^{\prime }$.