Open Access
Add to BibSonomyAbstract
A numerical model based on the successive order of scattering method has been developed. In this model, the vertical and angular distributions of underwater radiance are computed without using an empirical procedure, by assuming that a water body comprises numerous parallel horizontal layers. The model is validated using the observation data obtained from Lake Pend Oreille and Suruga Bay and for highly scattering water. The model results are in good agreement with the observation data. The model stability is demonstrated by assuming the single-scattering albedo for highly scattering water to be 1.0.
©2010 Optical Society of America
1. Introduction
For underwater light fields, radiance is a key parameter for describing underwater light fields since other parameters such as irradiance can be derived from the radiance distribution. With the development of algorithms for ocean-color remote sensing, the importance of determining the radiance distribution in the ocean is being well recognized. However, in ocean-color remote sensing, only the radiation distribution near the sea surface is considered. In order to understand and compute underwater light fields, it is essential to study the radiance distribution near the sea surface as well as that below the surface [1].
While downward irradiance and upward radiance have been frequently measured in various studies, radiance distribution measurements have not been pursued actively because of the difficulties associated with the measurement method and data calibration. The measurement process becomes tedious if a single radiometer is used, as time is lost in sweeping across the zenith and azimuth planes [2,3]. To solve the abovementioned problem, Voss and Chapin designed an instrument equipped with cooled charge-coupled devices (CCDs) on the basis of the fisheye technique [4]. Nevertheless, there are very few tabulated sets of radiance distribution data reported in the literature [2,3].
Various numerical models based on the Monte Carlo method, invariant imbedding method, and discrete ordinate method have been developed to determine the solutions of radiative transfer equations (RTEs) for underwater radiance distribution [5–7]. Mobley et al. [8] compared the light fields computed using the models based on the abovementioned methods and found that the models based on the invariant imbedding and discrete ordinate methods were stable and afforded good results. The result obtained with HydroLight, which has been developed on the basis of the invariant imbedding and quad-averaging methods, is often used as a reference standard for the development and/or validation of new models [10–12]. Since the models based on the invariant imbedding method and discrete ordinate method have various simplifying assumptions, they are mathematically complex. Such mathematical complexities make the physical interpretation of the RTEs difficult. Further, existing methods are hard t o estimate computation errors.
One of a natural solution is the successive order of scattering (SOS). The SOS method is consists of simple equations based on a phenomenological process of light. Therefore, a numerical model based on the SOS method can make physical interpretation of RTEs easily. Recently, atmosphere-ocean system models based on the SOS method have been developed [13,14]. These models were intended to be used for remote sensing and not for obtaining underwater radiance distribution data, implying that the results computed by these models were based on the radiance distribution measured at the top of the atmosphere or at the sea surface.
The SOS method has also been used to calculate the irradiance reflectance just below the ocean surface [15]. Kishino [16] has developed an SOS model for calculating the underwater radiance distribution; this model includes a single-scattering model proposed by Jerlov and Fukuda [17] that can be applied to a number of horizontal parallel plane thin layers. The SOS model is necessary to derive water body into a numerous number of thin layers in which single scattering is safely assumed. In Kishino’s SOS model, for instance, the layer thickness was determined by experience. Namely, when we apply this model for various optical properties of sea water, it meets serious problem because of no quantitative guide line concerning with appropriate layer thickness.
In this paper, we report the first results obtained with the newly developed model (TRAD) based on the SOS method. TRAD enables us to obtain the angular and vertical distribution of underwater radiance, and the results are in good agreement with the observation data. Using TRAD, we can determine the layer thickness without adopting the empirical approach. For validation of the model, we compared the radiance and irradiance computed by TRAD with the measurement data obtained for Lake Pend Oreille and Suruga Bay as well as with the results obtained with HydroLight [9]. However, it is difficult to use this numerical model to obtain the underwater radiance distribution in the case of highly scattering waters [8]. Hence, in the case of highly scattering waters, we compared the irradiance values computed using TRAD and HydroLight and validated the former by the theoretical knowledge gained from Gershun’s equation. In this study, we do not take the computation time into account since we aim to determine the underwater radiance distribution by using a mathematically simple method that makes the physical interpretation of underwater light propagation possible.
2. Method
2.1 Successive order of scattering (SOS) method
In the SOS method, the total radiance L(z,θ,ϕ) is calculated by summation of scattered radiance as follows:
where z is the depth (positive downward); θ, the zenith angle defined in the downward direction in this paper; ϕ, the azimuth angle; and N, is the maximum order of scatterings. In all the equations presented herein, the wavelength argument λ is omitted for brevity. Ln(z,θ,ϕ) indicates the n th scattered radiance and is given bywhere L p ,n is the n th scattered path radiance; L t ,n, the transmitted radiance. L p ,n scattered from the θ’, ϕ’ direction to the θ, ϕ direction is split into two terms in the following manner:where L d ,n and L u ,n are the radiances scattered from the (downward) θ’, ϕ’ and (upward) θ’, ϕ’ directions to the θ, ϕ, direction, respectively, and are expressed as follows:Here, ω 0 is the single-scattering albedo, which is the ratio of the scattering coefficient b and the beam attenuation coefficient c. , the scattering phase function from the θ’, ϕ’ direction to the θ, ϕ; Ω’ direction, is an infinitesimal solid angle in the direction θ’, ϕ’. Δz is the thickness of the horizontal plane parallel layers. z – and z + indicate the depth above and below the plane at a given z, respectively. L t ,n in Eq. (2) can be expressed as
where z 0 is the depth just below the sea surface; z b, the bottom depth; , the average beam attenuation coefficient between z 0 and z for θ < π/2 or that between z and z b for θ > π/2.When n = 0, L p , 0 is assumed to be zero, and only L t , 0 is taken into account owing to the sky radiance distribution L sky (including the sun radiance) at z 0.
The n th scattered radiance Ln can be derived by solving Eqs. (2)–(7) iteratively after solving Eq. (8). By substituting L n in Eq. (1), the radiance distribution can be obtained.
2.2 Influence of the second order scattering as a function of layer thickness
The influence of second-order scattering on radiant power is shown in Fig. 1 . The radiant power F 0 attenuates with the path length r , as indicated by the beam attenuation coefficient c. The attenuated radiant power F c can be derived
Let us consider the path length l and l – Δl (0 < l – Δl < l < r). The attenuation in the radiant power, F 1, from l – Δl to l can be determined from the relation
If there is no second-order scattering, F 1 remains constant between l and r. If second-order scattering occurs at l, F 1 reduces to F 2, which is given by
Using Eqs. (9)–(11) and by increasing Δl to r and making Δl infinitesimal, the attenuated radiant power F s, which is influenced by the second-order scattering at r, can be obtained from the following expression:
The parameter ε, which shows the influence of second-order scattering on the scattering is defined as
If we want to suppress the influence of second order scattering less than a given level, we can easily estimate it by solving Eq. (13) for cr.
2.3 Configuration of TRAD
TRAD assumes that a water body consists of a number of horizontal parallel layers of thickness Δz. The angular coordinates in TRAD are defined by classifying the directions on the basis of numerous segments bounded by the constants Δθ and Δϕ, similar to the quads in HydroLight, but there are no polar caps and equatorial quads in TRAD. The central angle θ, ϕ of a segment represents the direction of the radiance computed for that segment.
Following are the assumptions made in TRAD: (1) single scattering within a layer, (2) a flat sea surface, (3) no reflection from the sea bottom, (4) homogeneous scattering phase function through water body, (5) homogeneous inherent optical properties (IOPs) within a layer, and (6) homogeneous IOPs in the horizontal plane.
The input parameters for TRAD are absorption coefficient (a), scattering coefficient (b), scattering phase function (), sky radiance distribution including the sun’s radiance (L sky), angular resolutions (Δθ, Δϕ), depth interval (layer thickness, Δz), bottom depth (z b), and maximum order of scattering (N). In this paper, a and b are assumed to be homogeneous through water body. L sky is computed by using the algorithms proposed by Harrison and Coombes [18] and Gregg and Carder [19] as well as by using HydroLight. The angular resolutions Δθ, Δϕ of all the segments are set to 10° × 10°. By using the method described in detail in the following subsection (subsection 2.2), Δz is determined to be less than 1% of the influence of second-order scattering. Since TRAD cannot compute the infinite depth, z b is used as an input parameter. The default value of N is set as 1000. In order to accelerate the computation, the computation process is terminated even before N if the higher order scattering is negligible contribution on radiance, which is evaluated as
where Nz, Nθ, and Nϕ are the number of depths, vertical segments (18, if Δθ is set to 10°), and horizontal segments (36 if Δϕ = 10°).The values for each of the segments and the refraction and reflection at the sea surface were determined before the radiance calculation. Then, the scattering phase functions for the segments are computed for a resolution of 1° and stored.
As shown in Eq. (13), ε is a function of cr and layer thickness Δz ( = rcosθ) is expected to be dependent on c. In TRAD, Δz is adjusted depending on the longest path length, which is the angle θ representing the segment that is closest to the horizontal plane. For example, if ε is 1%, cr is set as 0.02. For the maximum path length, θ = 85° if Δθ = 10°. Then, if c is 0.1 m–1, Δz has to be set to a value less than 1.7 × 10−2 m. The radiance can be overestimated if Δz is assigned a large value, as the radiance does not attenuate, as is expected for the given value of Δz.
2.4 Data for validation
For validating the TRAD computation, the measurement data reported for Lake Pend Oreille by Tyler [2] and the data for Suruga Bay [20] were used. The IOPs data and value for Lake Pend Oreille were obtained from Tyler et al. [21] and Mobley [8], respectively. The a value and volume scattering function (VSF) for Suruga Bay were measured using our own instruments, while the downward irradiance and upward radiance L u were measured using PRR-600 (Biospherical Instruments Inc., San Diego). The VSF was measured for the range 10°–160°. Since the shape of VSF was very similar to that obtained by Petzold for the water in San Diego Harbor [22], the phase function determined from the VSF for San Diego Harbor water was used for our computation. Table 1 lists the values of a, b, Δz, solar zenith angle θ s, and wavelength λ used for the computation.
Table 1. Input parameters for computation
In TRAD, we require a value of the bottom depth from which there is no reflection. Since the influence of reflected radiance by the sea bottom is not considered in this study, the maximum depth is set more 10m than the maximum depth of measurement data, i.e., 80 m for Lake Pend Oreille and 40 m for Suruga Bay.
HydroLight 4.3 has been used as a reference for our numerical computation. This is because TRAD and HydroLight have been so designed that their input parameters—a, b, , and θ s—have the same values, and both these models assume that the IOPS are homogeneous from the sea surface to the sea bottom and that the sea surface is flat.
In general, numerical radiative transfer computations have the difficulty to handle non-absorbing medium with high forward scattering peak. To see the performance of our model, we made the computation for highly scattering water, the input parameters used in TRAD and HydroLight are a = 0.0 and b = 1.0, which implies that ω 0 = 1.0; by Mobley et al. [8] is used; Mobley et al. [8] compared various numerical models by considering ω 0 = 0.9, while Gallegos et al. [23] reported the value of ω 0 to be 0.987 (the expression b/a ( = 78.9) was reported in their article) for an estuary. The irradiance results obtained with TRAD and HydroLight were compared, since Gershun’s equation could be used for the theoretical validation of our model. Gershun’s equation for horizontally stratified water column is
where E d and E u are the downwelling and upwelling irradiances, respectively; E 0d and E 0u are the downwelling and upwelling scalar irradiances, respectively. If there is no absorption, the following relation of irradiance should be held from Eq. (17):3. Results
3.1 Comparison with observation data
A comparison of the relative radiance distribution computed by TRAD and HydroLight and that observed at Lake Pend Oreille is shown in Fig. 2 . The radiance was normalized at a depth of 53.7 m. The normalization angles for the radiance distribution were 170°, 165.5°, and 170° for Lake Pend Oreille, TRAD, and HydroLight, respectively. The radiance distributions obtained by TRAD were in quantitative agreement with the observation data. That is, there was no significant difference in the relative radiances at different depths (up to 53.7 m), as shown in Fig. 2(a). For a depth of 4.24 m, a peak around 25° was corresponding to the sun radiance. Further, the steep curve at a shallow depth and the symmetric curve at a greater depth had similar angular shapes. The peak indicating the influence of the edge of optical manhole, which occurs due to refraction at the water surface, would appear at around ± 48° if the refractive index was 1.341; similarly, in the case of TRAD, an optical manhole peak is accurately re-produced, small shoulder at –45°. The difference is caused by angular resolution of TRAD. The data obtained for Lake Pend Oreille and HydroLight are compared in Fig. 2(b). The relation between the radiance distribution features for HydroLight and Lake Pend Oreille, such as peak angle and angular shape, was almost identical to that between the radiance distribution features for Lake Pend Oreille and TRAD. The radiance distributions obtained by TRAD and HydroLight and those levels were almost the same at all depths [Fig. 2(c)].
Fig. 2 Computed and observed radiance distribution in the plane of the sun. (a) solid line: TRAD, dashed line: Lake Pend Oreille (b) solid line: HydroLight, dashed line: Lake Pend Oreille (c) solid line: TRAD, dashed line: HydroLight
Figure 3 shows a comparison of the relative downward irradiance E d and L u observed at Suruga Bay and those computed by TRAD. All the data were normalized by the downward irradiance just below the sea surface. E d and L u computed by TRAD were confirmed to be consistent with the observation data up to a depth of around 15 m. At depths greater than 15 m, the decrease in E d and L u in the case of TRAD was slightly different from that in the observation data. Figure 3 also shows the E d and L u values only up to 50 orders of scattering, whereas in reality, E d and L u were determined by summation of the scattered radiances up to the 126th order. From the E d and L u profiles, it could be seen that higher order of scattering becomes dominant as a function of deeper. As shown in Fig. 3(a), the scattered E d profiles pronounce peaks of which depth is a function of scattering order. On the other hand, the scattered L u profiles do not show the significant peak. The depth of the peak increased with the scattering order. On the other hand, in the vertical profile of L u [Fig. 3(b)], there was no significant peak at any order.
Fig. 3 Vertical distributions of downward irradiance (a) and upward radiance (b) observed in Suruga Bay (circle) and calculated by TRAD (black line); downward irradiance and upward radiances at each order of scattering up to the 50th order (rainbow-colored line) are also shown.
3.2 Influence of layer thickness setting
The data shown in Fig. 4 was diffuse attenuation coefficient (K d) in depth 5–10 m and, ε corresponding to Δz. All the input parameters except Δz were the same as those used for the computation at Suruga Bay: c = 0.6 m–1; θs = 26° degree. ε is a logistic function of Δz. And, K d showed a decreasing function for Δz.
Fig. 4 Relation between K d, ε, and Δz when c = 0.6 m–1 (ω 0 = 0.75). K d was calculated by E d in the depths of 5m and 10m.
3.3 Case of highly scattering water
The E d and E u profiles computed by HydroLight and TRAD are shown in Fig. 5 . The absorption coefficient profiles are also shown, which are obtained by computed E d, E u, E 0d and E 0u using Eq. (17). The difference between the E d values computed by HydroLight and TRAD was less than 4.5%, and the difference between the E u values computed by these models was less than 1.7% [Fig. 5(a)]. The difference between E d and E u (E d – E u) was constant in depth but this difference was 0.271 for HydroLight and 0.259 for TRAD [Fig. 5(b)]. The absorption coefficients obtained with both these models were less than 1.0 × 10−5, indicating that TRAD could be inversed absorption coefficient from the computed irradiances [Fig. 5(c)]. The peaks in the results obtained for HydroLight were due to significant digits.
Fig. 5 Irradiances computed by HydroLight (red) and TRAD (blue) at assuming no absorption (ω 0 = 1.0) with high forward scattering peak. (a) Vertical profile of relative downward and upward irradiances; (b) difference between E d and E u; (c) absorption coefficient derived by Eq. (17).
4. Discussion
The discrepancies between the observation data and the results obtained with TRAD were negligible in the case of Lake Pend Oreille (Fig. 2). Tyler’s data were used for several corrections since the actual measurement process was time-consuming. The optical properties in the observation data were not strictly homogeneous because of the slight difference in the concentration of particulate matter in the vertical direction, as indicated by the observation data [2]. For the TRAD computation, we used Mobley’s scattering phase function [6], which may be different from the actual value. HydroLight also used the same scattering phase function, and the relation between the HydroLight results and the observation data was the same as that between the TRAD results and the observation data. This confirmed that the radiance distribution computed by TRAD was in good agreement with the observation data.
In the case of Suruga Bay, a marked difference between the observation data and the TRAD results can be seen at depths greater than 15 m. The TRAD computation is carried out under the assumption that there is no reflection from the sea bottom, whereas the bottom reflectance may have some influence on the vertical profiles of irradiance and radiance distribution; this is because the bottom depth is only around 30 m [20]. A similar tendency has been shown by Saruya et al. in their Monte Carlo simulation [20].
For the validation of Suruga Bay data, TRAD computed the downward irradiance and upward radiance with 1.0 × 10−3 m of Δz. From K d (Fig. 4), the Δz setting considered to be enough thin. Also, the Δz setting for Lake Pend Oreille (1.0 × 10−3 m) was sufficiently thin, since ε for Lake Pend Oreille was smaller than that for Suruga Bay. Since the results by TRAD showed a good agreement with the observation data and computed radiance distribution, the determination of Δz based on Eq. (13) is considered to be effective for the SOS method.
In order to see the numerical stability of the model, we made the computation assuming non-absorbing medium (ω 0 = 1) with highly forward scattering peak. The vertical profiles of E d and E u obtained by TRAD and HydroLight are almost identical. As expected from Eq. (18), the vertical profiles indicate that E d – E u is constant. Further, input parameters for the computation can be obtained by converting the irradiance values to absorption coefficients by using Eq. (17). On the basis of these results, we state that the theoretical consistency of TRAD is unaffected even for highly scattering water.
5. Conclusion
In the present study, a numerical model (TRAD) based on the SOS method has been developed. Since TRAD consists of simple equations whose physical interpretation is easy, it is considered a natural solution of RTEs. Despite the use of simple equations, the results obtained with TRAD are in good agreement with the observation data. Further, the results obtained with TRAD are consistent with those obtained with HydroLight. Since in TRAD, the layer thickness is quantitatively evaluated with the help of the error estimator, the model can predict underwater light fields with sufficient accuracy for various oceanic conditions. In addition, since its physical interpretation is simple, TRAD can be applied to a wide range of oceanographic problems. TRAD will be a useful tool for developing algorithm for ocean color remote sensing, for instance, as a forward model of inverse modeling like neural network, I understand the value of the model. In this study, several assumptions have been made in TRAD. One of these assumptions has already been proved, that is, the IOPs are inhomogeneous from the sea surface to the sea bottom. The results obtained with the new TRAD model will be presented in the near future. The other assumptions such as flat sea surface and zero reflection from the bottom will be confirmed, as has been done for numerical models other than TRAD and HydroLight.
Acknowledgments
This study was supported by a grant from the Japan Aerospace Exploration Agency (JAXA), JX-PSPC-227026, and by a Grant-in-Aid for Scientific Research on Priority Areas21014002. Fruitful discussions with Dr. M. Kishino and Dr. T. Oishi of Tokai University are appreciated.
References and links
1. J. R. V. Zaneveld, A. Barnard, and E. Boss, “Theoretical derivation of the depth average of remotely sensed optical parameters,” Opt. Express 13(22), 9052–9061 (2005). [CrossRef] [PubMed]
2. J. E. Tyler, “Radiance distribution as a function of depth in the submarine environment,” SIO Ref. 58–25, pp.37 (1958).
3. E. Aas, N. K. Højerslev, and B. Lundgren, “Spectral irradiance, radiance and polarization data from the Nordic Cruise in the Mediterranean Sea during June-July 1971,” Rep. Dep. of Geophys, Univ. of Oslo 102, 97 (1997).
4. K. J. Voss and A. L. Chapin, “Upwelling radiance distribution camera system, NURADS,” Opt. Express 13(11), 4250–4262 (2005). [CrossRef] [PubMed]
5. C. N. Adams and G. W. Kattawar, “Effect of volume scattering function on the errors induced when polarization is neglected in radiance calculation in an atmosphere-ocean system,” Appl. Opt. 32(24), 4610–4617 (1993). [CrossRef] [PubMed]
6. C. D. Mobley, “A numerical model for the computation of radiance distribution in natural waters with wind-roughened surfaces,” Limnol. Oceanogr. 34(8), 1473–1483 (1989). [CrossRef]
7. Z. Jin and K. Stamnes, “Radiative transfer in nonuniformly refracting layered media: atmosphere-ocean system,” Appl. Opt. 33(3), 431–442 (1994). [CrossRef] [PubMed]
8. C. D. Mobley, B. Gentili, H. R. Gordon, Z. Jin, G. W. Kattawar, A. Morel, P. Reinersman, K. Stamnes, and R. H. Stavn, “Comparison of numerical models for computing underwater light fields,” Appl. Opt. 32(36), 7484–7504 (1993). [CrossRef] [PubMed]
9. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic, San Diego, Calif., 1994).
10. Z. Lee, K. L. Carder, C. D. Mobley, R. G. Steward, and J. S. Patch, “Hyperspectral remote sensing for shallow waters. I. A semianalytical model,” Appl. Opt. 37(27), 6329–6338 (1998). [CrossRef]
11. Z. Lee, K. L. Carder, C. D. Mobley, R. G. Steward, and J. S. Patch, “Hyperspectral remote sensing for shallow waters. 2. Deriving bottom depths and water properties by optimization,” Appl. Opt. 38(18), 3831–3843 (1999). [CrossRef]
12. A. Albert and C. D. Mobley, “An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters,” Opt. Express 11(22), 2873–2890 (2003). [CrossRef] [PubMed]
13. M. Chami, R. Santer, and E. Dilligeard, “Radiative transfer model for the computation of radiance and polarization in an ocean-atmosphere system: polarization properties of suspended matter for remote sensing,” Appl. Opt. 40(15), 2398–2416 (2001). [CrossRef]
14. P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17(4), 2057–2079 (2009). [CrossRef] [PubMed]
15. A. Morel and L. Prieur, “Analysis of variation in ocean color,” Limnol. Oceanogr. 22(4), 709–722 (1977). [CrossRef]
16. M. Kishino, “Numerical calculation of radiative transfer in the sea,” Mer (Paris) 12, 26–33 (1974).
17. N. G. Jerlov and M. Fukuda, “Radiance distribution in the upper layers of the sea,” Tellus 12(3), 348–355 (1960). [CrossRef]
18. A. W. Harrison and C. A. Coombes, “An opaque cloud cover model of sky short wavelength radiance,” Sol. Energy 41(4), 387–392 (1988). [CrossRef]
19. W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35(8), 1657–1675 (1990). [CrossRef]
20. Y. Saruya, T. Oishi, M. Kishino, Y. Jodai, K. Kadokura, and A. Tanaka, “Influence of ship shadow on underwater irradiance fields,” Proc. SPIE 2963, 760–765 (1996). [CrossRef]
21. J. E. Tyler, W. H. Richardson, and R. W. Holmes, “Method for obtaining the optical properties of large bodies of water,” J. Geophys. Res. 64(6), 667–673 (1959). [CrossRef]
22. T. J. Petzold, “Volume scattering functions for selected ocean waters,” SIO Ref. 72–78, 79pp (1972).
23. C. L. Gallegos, D. L. Correll, and J. W. Pierce, “Modeling spectral diffuse attenuation, absorption, and scattering coefficients in a turbid estuary,” Limnol. Oceanogr. 35(7), 1486–1502 (1990). [CrossRef]
