Page updated: May 3, 2021
Author: Curtis Mobley
View PDF

The Single-Scattering Approximation

As previously noted, exact analytical solutions of the RTE exist only for a few idealized and unphysical situations such as no scattering. There are, however, a few approximate analytic solutions. In pre-computer days these were useful computational tools. These approximate solutions are no longer needed for numerical computation, but they are still useful for isolating the most important processes governing light propagation in the ocean and can provide guidance in interpretation of radiometric data. This page develops one such solution: the single-scattering approximation (SSA). The next page discusses the related quasi-single scattering approximation (QSSA).

The Successive Order of Scattering Solution Technique

We begin with the optical depth form of the time independent, 1D (plane parallel geometry) RTE, which is Eq. (4) of the Scalar Radiative Transfer Equation page:

μdL(ζ,μ,ϕ,λ) dζ = L(ζ,μ,ϕ,λ) + ωo(ζ,λ)02π11L(ζ,μ,ϕ,λ)β̃(ζ; μ,ϕ μ,ϕ; λ)dμdϕ + 1 c(ζ,λ)Σ(ζ,μ,ϕ,λ). We next make a number of simplifications by assuming that
  • The water is homogeneous, so that the IOPs do not depend on depth;
  • The water is infinitely deep;
  • The sea surface is level (zero wind speed);
  • The sun is a point source is a black sky, so that the incident radiance onto the sea surface is collimated;
  • There are no internal sources or inelastic scattering.

The RTE then becomes, for a given wavelength λ, which we henceforth drop for brevity,

μdL(ζ,μ,ϕ) dζ = L(ζ,μ,ϕ) + ωo02π11L(ζ,μ,ϕ)β̃(μ,ϕ μ,ϕ)dμdϕ. (1)

A powerful technique for solving differential equations is to attempt a power series solution in which higher order terms of the series are weighted by a powers of a parameter whose magnitude is less than 1. The higher order terms then contribute less and less to the sum that represents the solution. The albedo of single scattering, ωo, meets the requirement for an expansion parameter. We therefore attempt a solution of Eq. (1) of the form

L(ζ,μ,ϕ) = k=0ω okL(k)(ζ,μ,ϕ) = L(0)(ζ,μ,ϕ) + ω oL(1)(ζ,μ,ϕ) + ω o2L(2)(ζ,μ,ϕ) + (2) The notation L(0) denotes radiance that is unscattered, L(1) is radiance from rays that have been scattered once, L(2) is radiance from rays that have been scattered twice, and so on. This is consistent with the interpretation of ωo as the probability of ray survival in an interaction with matter, i.e., the probability that a ray will be scattered and not absorbed.

We now substitute Eq. (2) for the radiance into Eq. (1) to obtain

μ dL(0) dζ + ωodL(1) dζ + ωo2dL(2) dζ + = [L(0) + ω oL(1) + ω o2L(2) + ] + ωo02π11[L(0) + ω oL(1) + ω o2L(2) + ]β̃(μ,ϕ μ,ϕ)dμdϕ. (3) We next group terms that have the same power of ωo:

μdL(0) dζ + L(0) + ωo μdL(1) dζ + L(1) 02π11L(0)β̃(μ,ϕ μ,ϕ)dμdϕ + ωo2 μdL(2) dζ + L(2) 02π11L(1)β̃(μ,ϕ μ,ϕ)dμdϕ + = 0. This equation must hold true for any value of 0 ωo < 1. Setting ωo = 0 would leave only the first line of the equation, whose terms must sum to 0. Similarly, when ωo0, each group of terms multiplying a given power of ωo must equal zero in order for the entire left side of the equation to sum to zero. We can therefore equate to zero the groups of terms in brackets multiplying each power of ωo. This gives a sequence of equations:

μdL(0) dζ = L(0) (S0) μdL(1) dζ = L(1) +02π11L(0)β̃(μ,ϕ μ,ϕ)dμdϕ (S1) μdL(2) dζ = L(2) +02π11L(1)β̃(μ,ϕ μ,ϕ)dμdϕ (S2) and so on. Note that because ωo multiples the path integral term in Eq. (3), the path integrals in this sequence of equations always involve the radiance at one order of scattering less than the derivative term. We first solve Eq. (S0), which governs the unscattered radiance. The solution for L(0) then can be used in Eq. (S1) to evaluate the path integral, which becomes a source function for singly scattered radiance. After solving Eq. (S1) for singly-scattered radiance, L(1) can be used to evaluate the path function in Eq. (S2), and so on. This process constitutes the successive-order-of-scattering (SOS) solution technique.

Solution of Eq. (S0) for the unscattered radiance

To solve (S0) we need boundary conditions at the sea surface and bottom. Figure 1 reminds us that the incident unscattered radiance onto the sea surface, and transmitted into the water, is perfectly collimated because we have assumed that the sun is a point source in a black sky and the surface is level. In that figure, E(0) denotes the irradiance measured just below the sea surface on a plane that is perpendicular to the direction of photon travel (denoted by the red dashed line), and 𝜃sw is the Sun’s zenith angle in the water after refraction by the level surface.


Figure 1: The collimated incident radiance onto the sea surface and the refracted radiance just below the surface.

Recalling the Dirac delta function, we can write the unscattered radiance just below the surface as

L(0)(0,μ,ϕ) = E (0)δ(μ μsw)δ(ϕ ϕsw), (BC1)

where (μsw,ϕsw) is the direction of the Sun’s beam in the water. The two delta functions, which together have units of sr1, “pick out” the direction of the Sun’s beam; the unscattered radiance is zero in all other directions. Note that integrating this radiance over all downward directions to compute the downwelling plane irradiance gives

Ed(0) = 0102πE (0)δ(μ μsw)δ(ϕ ϕsw)μdμdϕ = E(0)μsw as expected.

It is assumed that the incident solar irradiance is given, so Eq. (BC1) is the boundary condition on L(0)(ζ,μ,ϕ) at the sea surface (i.e., in the water at depth ζ = 0). We are assuming that the water is infinitely deep and source free, so the radiance must approach 0 at great depth. The boundary condition at the bottom is thus

L(0)(ζ,μ,ϕ) 0asζ . (BC2)

We can now solve Eq. (S0) subject to boundary conditions (BC1) and (BC2). Rewriting (S0) as

dL(0)(ζ) L(0)(ζ) = dζ μ

and integrating from depth 0 to ζ, corresponding to radiances L(0)(0) and L(0)(ζ) respectively, gives

ln L(0) L(0)(0)L(0)(ζ) = ζ μ0ζ


L(0)(ζ,μ,ϕ) = L(0)(0,μ,ϕ)eζμ (4a) = E(0)δ(μ μsw)δ(ϕ ϕsw)eζμ. (4b)

Solution (4a) is simply the Lambert-Beer law: the initial unscattered radiance decays exponentially with optical depth. Using (BC1) to rewrite the radiance at the surface gives (4b), which will be the convenient form for solution of (S1) below. Equation (4b) also shows explicitly that the unscattered radiance is 0 except in direction (μsw,ϕsw). The exponential forces the radiance to 0 as the depth increases, so that (BC2) is satisfied. Thus our solution satisfies both the surface and bottom boundary conditions and thus constitutes a complete solution of the two-point boundary value problem for unscattered radiance. This solution gives the contribution of unscattered radiance to the total radiance.

Solution of Eq. (S1) for the singly scattered radiance

The first step in solving (S1) is to evaluate the scattering term using the solution for L(0). To do this we use (4b) to get

02π11L(0)(ζ,μ,ϕ)β̃(μ,ϕ μ,ϕ)dμdϕ = 02π11E (0)δ(μ μ sw)δ(ϕ ϕ sw)eζμ β̃(μ,ϕ μ,ϕ)dμdϕ = E(0)eζμsw β̃(μsw,ϕsw μ,ϕ). (5) This result shows how much of the unscattered radiance reaches depth ζ and then gets scattered into the direction of interest (μ,ϕ). In other words, the unscattered radiance is a local (at depth ζ) source term for singly scattered radiance.

All quantities on the right hand side of Eq. (5) are known from the given IOPs and surface boundary condition. We can therefore proceed with the solution of (S1) for the singly scattered radiance L(1). The equation to be solved is

μdL(1) dζ + L(1) = E (0)eζμsw β̃(μsw,ϕsw μ,ϕ), (6)

where the right hand side is now a known function of depth. There is no incident scattered radiance from the sky because the sun’s collimated beam is all unscattered light. Thus the boundary conditions for Eq. (6) are

L(1)(0,μ,ϕ) = 0andL(1)(ζ,μ,ϕ) 0asζ . (7)

Figure 2 shows that the singly scattered downwelling radiance at depth ζ comes only from above depth ζ, and that the upwelling radiance at ζ comes only from depths below ζ. We can thus consider the downwelling, Ld(1)(ζ,μ,ϕ), and upwelling, Lu(1)(ζ,μ,ϕ), radiances separately. We can integrate from the surface down to ζ to compute Ld(1), and we can integrate from ζ to to compute Lu(1).


Figure 2: Single-scattering contributions to Ld(1) and Lu(1).

If you were paying attention in your undergraduate differential equations class, you recognize Eq. (6) as an ordinary differential equation with constant coefficients, which can be solved by means of an integrating factor. Multiplying Eq. (6) for downwelling radiance by 1 μeζμ (the integrating factor) gives

1 μeζμ μdLd(1)(ζ) dζ + Ld(1)(ζ) = 1 μeζμ E β̃eζμsw d dζ Ld(1)(ζ)eζμ = Eβ̃ μ exp 1 μ 1 μsw ζ. (8) Now integrating from depth 0 to ζ, where the radiances are Ld(1)(0) and Ld(1)(ζ), respectively, and recalling that Ld(1)(0) = 0 by the upper boundary condition (7) gives
Ld(1)(ζ)eζμ = Eβ̃ μ 1 1 μ 1 μsw exp 1 μ 1 μsw ζ 1 ,

provided that μμsw. Recalling that Ed(0) = E(0)μsw, the preceding equation can be rewritten as

Ld(1)(ζ,μ,ϕ) = E d(0)β̃(μsw,ϕsw μ,ϕ) 1 μsw μ eζμsw eζμ . (9)

For the special case of μ = μsw but ϕϕsw, so that the scattering angle is nonzero, Eq. (8) reduces to

d dζ Ld(1)(ζ)eζμsw = Eβ̃ μsw

which integrates to

Ld(1)(ζ,μ sw,ϕ) = Eβ̃(μsw,ϕsw μsw,ϕ) ζ μsweζμsw = Ed(0)β̃(μsw,ϕsw μsw,ϕ) ζ μsw2eζμsw . (10) The second form results from Ed(0) = E(0)μsw, which was derived above.

The direction of μ = μsw and ϕ = ϕsw is the case of no scattering, so there is no singly scattered radiance.

We next compute the upwelling radiance at ζ by integrating Eq. (6) from ζ to , keeping in mind that now μ = cos 𝜃 < 0 since 𝜃 is measured from 0 in the nadir direction. The integration gives (writing μ = |μ| to emphasize the negativity of μ)

Lu(1)(ζ)eζ|μ| ζ Lu(1)(ζ)eζ|μ| = Eβ̃ |μ| 1 1 |μ| 1 μsw exp 1 |μ| 1 μsw ζ ζ exp 1 |μ| 1 μsw ζ Both limits as ζ are zero. The result can be rewritten as
Lu(1)(ζ) = E d(0)β̃(μsw,ϕsw μ,ϕ) 1 μsw μeζμsw . (11)

Assembling the SSA solution

Recalling from Eq. (2) that the SSA is given by

L(SSA)(ζ,μ,ϕ) = L(0)(ζ,μ,ϕ) + ω oL(1)(ζ,μ,ϕ),

we can assemble L(SSA) from the pieces computed in Eqs. (4a) and (9-11):

Ld(SSA)(ζ,μ,ϕ) = L(0)(0,μ sw,ϕsw)eζμsw (12) ifμ = μswandϕ = ϕsw Ld(SSA)(ζ,μ,ϕ) = ω oEd(0)β̃(μsw,ϕsw μsw,ϕ) ζ μsw2eζμsw (13) ifμ = μswbutϕϕsw Ld(SSA)(ζ,μ,ϕ) = ω oEd(0)β̃(μsw,ϕsw μ,ϕ) 1 μsw μ eζμsw eζμ (14) ifμ > 0andμμsw Lu(SSA)(ζ,μ,ϕ) = ω oEd(0)β̃(μsw,ϕsw μ,ϕ) 1 E(0)μsw μeζμsw (15) ifμ 0 Equations (12)-(15) constitute the SSA solution to the RTE. This solution is seen, for example, in Gordon (1994), where it is presented without derivation.

It is easy to show that

lim μμsw 1 μsw μ eζμsw eζμ = ζ μsw2eζμsw , (16)

in which case Eq. (14) reduces to Eq. (13), which was derived independently as a special case of the depth integration.

It must be remembered that the SSA rests upon a number of simplifying assumptions. In particular, the input sky radiance was collimated. The delta functions in direction then made evaluation of the scattering path function in Eq. (5) easy. This would not be the case for any other sky radiance distribution, or for a non-level sea surface. Likewise, the assumption of infinitely deep water removed any bottom effect.

The SSA will be a good approximation to actual radiances only if the higher order terms in the Eq. (2) are negligible. This means that ωo must be sufficiently small, but how small? Figure 3 compares Lu(SSA) and Ld(SSA) with radiances computed by HydroLight for nadir- and zenith-viewing radiances. The sun was at 42 deg, which gives the in-water solar zenith angle of 𝜃sw = 30 deg or μsw = 0.866. This gives a scattering angle of ψ = 30 deg for Ld(SSA) and ψ = 150 deg for Lu(SSA). The Petzold “average-particle” phase function) was used, for which β̃(ψ = 30) = 0.08609 and β̃(ψ = 150) = 0.002365sr1. HydroLight includes all orders of multiple scattering, so comparison of its radiances with the SSA values shows the importance of multiple scattering. The HydroLight runs modeled the SSA conditions as closely as possible, the difference being that the SSA is for one exact direction and HydroLight computes nadir and zenith radiances as averages over polar caps with a 5 deg half angle, and the sun’s direct beam in water is spread out over a quad from 𝜃 = 25 to 35 deg. The HydroLight runs set Ed(inair) = 1.028Wm2sr1 so that Ed(0) = 1.0Wm2sr1.


Figure 3: Comparison of the SSA and HydroLight solutions for nadir-viewing (Lu; dashed lines) and zenith-viewing (Ld; solid lines) radiances for four ωo values; Ed(0) = 1Wm2sr1.

Figure 3 shows that for ωo = 0.01 the agreement between the SSA and HydroLight is very good. The HydroLight values are slightly higher than the SSA values because there is still a small multiple scattering contribution to the total radiance even at very small ωo values. For ωo = 0.1 the SSA still gives good results near the sea surface but differs from the multiple scattering solution by a factor of 3 at 10 optical depths. For ωo = 0.85, which is typical of blue and green wavelengths in ocean waters, the SSA upwelling radiance is a factor of five too small even at the surface, and the SSA radiances are off by orders of magnitude at large optical depths. Thus, as expected, we see that the SSA is of little use in optical oceanography because multiple scattering almost always dominates underwater radiance distributions at visible wavelengths.

We end the SSA discussion by noting that Walker (1994) has carried the SOS solution through second order scattering. His development requires a good bit of mathematical masochism and results in a much more complicated set of equations, which can be seen in his Section 2-6. There is little need for such approximations given the ease of numerical solution of the RTE to include all (in HydroLight) or at least many (in Monte Carlo models) orders of multiple scattering, and without any of the assumptions required for the analytic evaluation of the path integrals in the SSA. Perhaps the greatest value of the SSA solution is that it can be used to check numerical models when ωo is small.

Comments for The Single-Scattering Approximation:

Loading Conversation