Fast light calculations for ocean ecosystem and inverse models

Curtis D. Mobley

doi:10.1364/OE.19.018927

1. Introduction

Currently available ocean ecosystem models often use very sophisticated treatments of the hydrodynamics that employ curvilinear coordinate systems and primitive-equation solutions to predict advection and upper-ocean thermodynamics and mixing. The formulations of the biology are becoming increasingly sophisticated with multiple phytoplankton functional groups and complex connections between primary production, nutrient utilization, and grazing. However, most ecosystem models still use grossly oversimplified and often inaccurate treatments of the optics. The optics component of coupled ecosystem models is sometimes just a single equation parameterizing photosynthetically available radiation (PAR) or the scalar irradiance in terms of the chlorophyll concentration and a few parameters such as the solar zenith angle and an assumed Jerlov water type. Such simple models often fail even in Case 1 waters, and they can be wrong by orders of magnitude in Case 2 or optically shallow waters. There is thus a need for a radiative transfer numerical model that can be called as a subroutine within ecosystem models to bring the optics component up to the level of accuracy and sophistication needed for simulations of any water body, including Case 2 and optically shallow waters.

Inverse models seek to recover ocean absorption and scattering properties from in-water optical measurements such as the downwelling irradiance and upwelling radiance, or from the remote-sensing reflectance. Implicit inverse models do this by repeated solution of the radiative transfer equation (RTE) as the input absorption and scattering properties are systematically varied starting with some initial guess. At each iteration the RTE-predicted values of the optical quantities are compared with their measured values, and the current IOP values are adjusted according to the difference in predicted and measured values. This process may require hundreds of RTE solutions before convergence is obtained.

The widely used HydroLight [1–3] radiative transfer software solves the unpolarized RTE with high accuracy to compute in-water and water-leaving radiance distributions L(z, θ, ϕ, λ) as functions of depth z, polar θ and azimuthal ϕ directions, and wavelength λ for any plane-parallel water body. It therefore meets the need for generality in ecosystem and implicit inverse models. However, HydroLight is used primarily as a research tool and its emphasis is on accuracy of the RTE solution, with run time being of secondary importance. Moreover, it computes the full angular distribution of the radiance. It is therefore computationally much too expensive to be practicable for ecosystem models that require light calculations at thousands of grid points and time steps or for inverse models that require hundreds of RTE solutions to achieve convergence.

EcoLight-S(ubroutine) was developed to meet the need for a computationally fast radiative tranfer code that emphasizes extremely fast run times at the expense of accuracy of the RTE solution. EcoLight-S has various user-selectable options such as solving the RTE at only selected wavelengths and to different depths at different wavelengths, with values at unsolved wavelengths and depths being obtained by interpolation or extrapolation. Such options can greatly reduce the run time, but also reduce the accuracy of the computed irradiances and other quantities.

The inputs for EcoLight-S are the same as for HydroLight or any other ocean radiative transfer code: the inherent optical properties (IOPs, namely the absorption and scattering properties) of the water body, the sky radiance incident onto the sea surface, a sea-surface boundary condition parameterized by the wind speed, and a bottom boundary condition parameterized by the bottom reflectance (in finite-depth waters). All inputs needed to solve the RTE must be explicitly defined by the user before calling the EcoLight-S subroutine. Unlike HydroLight, which has a user interface with many options for defining the inputs needed to solve the RTE, EcoLight-S has no user interface or built-in bio-geo-optical IOP, sky, or bottom reflectance models. However, such models can be called by the user’s own program to obtain the inputs needed by EcoLight-S, and examples are provided with the code.

EcoLight-S is entirely new code written in Fortran 95, except for a few public-code Fortran 77 legacy routines for standard mathematical operations. All communication between the user’s code and EcoLight-S is via two Fortran 95 modules and one call to subroutine EcoLight-S. The input module contains all of the inputs that must be defined by the user (on a “fill-in-the-blanks” format) before calling EcoLight-S, and the output module returns all of the outputs from EcoLight-S for use by the user. The user does not need to have any knowledge of the internal workings of the EcoLight-S code. Likewise, the EcoLight-S code is completely independent of the nature and purpose of the user’s code. The programming details of how to specify the EcoLight-S inputs are given in the user’s guide [4].

Section 2 describes the formulation of the RTE solved by EcoLight-S. Section 3 discusses various options for reducing the computation time. Section 4 then presents example simulations showing how EcoLight-S can be optimized to reduce run time while keeping the computed quantities within acceptable error bounds.

2. RTE formulation

HydroLight solves the depth-dependent RTE

\begin{array}{l} cos θ \frac{d L (z, θ, ϕ, λ)}{d z} & = & - c (z, λ) L (z, θ, ϕ, λ) \\ + & \int_{0}^{2 π} \int_{0}^{π} β (z, θ^{'}, ϕ^{'}, θ, ϕ, λ) L (z, θ^{'}, ϕ^{'}, λ) sin θ^{'} d θ^{'} d ϕ^{'} + S (z, θ, ϕ, λ) . \end{array}

Here c is the beam attenuation coefficient (the sum of the absorption and scattering coefficients); β is the volume scattering function (VSF), which describes how radiance is scattered from direction (θ′, ϕ′) into direction (θ, ϕ); and S is a source function representing the contributions of inelastic scattering at wavelength λ. Depth z is measured positive downward from the mean sea surface, θ = 0 corresponds to light heading straight down, and ϕ = 0 can be chosen for convenience, e.g. in the sun’s azimuthal direction.

The full radiance distribution L(z, θ, ϕ, λ) obtained from Eq. (1) gives far more information than is needed by ecosystem models. The spectral total scalar irradiance

E_{o} (z, λ) = \int_{0}^{2 π} \int_{0}^{π} L (z, θ, ϕ, λ) sin θ d θ d ϕ

is the quantity needed for calculations of photosynthesis and heating of the water. Indeed, most ecosystem biological models parameterize the light in terms of PAR,

P A R (z) = \int_{400}^{700} \frac{λ}{h c} E_{o} (z, λ) d λ .

The λ/hc factor, where h is Planck’s constant and c is the speed of light, converts the energy units of E _o (W m⁻² nm⁻¹) to quantum units (photons m⁻² s⁻¹) as needed for photosynthesis. Water heating rates by short-wave radiation are usually computed by

\frac{\partial T (z, t)}{\partial t} = \frac{1}{ρ c_{p}} \frac{\partial [{\bar{E}}_{d} (z, t) - {\bar{E}}_{u} (z, t)]}{\partial z},

where T is the temperature, t is time, ρ is the density of sea water, c_p is the specific heat of sea water, and Ē_d and Ē_u are respectively the downwelling and upwelling plane irradiances integrated over visible and near-infrared wavelengths, typically 400 to 1000 nm. These irradiances are computed from the radiance by

{\bar{E}}_{d} (z) = \int_{400}^{1000} \int_{0}^{2 π} \int_{0}^{π / 2} L (z, θ, ϕ, λ) cos θ sin θ d θ d ϕ d λ,

with a similar equation over the upward directions for Ē_u.

It is thus clear that only irradiances are necessary for ecosystem modeling. These irradiances are obtained by integration of the radiance over all azimuthal directions ϕ, as seen in Eqs. (2) and (5). The azimuthal dependence of the radiance, obtained at considerable computational expense in HydroLight, is thus lost when computing irradiances. In most waters Raman scatter by water and fluorescence by chlorophyll and colored dissolved organic matter (CDOM) have little effect on PAR (see Table 2 below and the discussion thereof). Including these inelastic scatter processes requires additional inputs and significantly increases run times. It is therefore reasonable to solve an azimuthally integrated, source-free version of the RTE,

cos θ \frac{d L (z, θ, λ)}{d z} = - c (z, λ) L (z, θ, λ) + \int_{0}^{π} β (z, θ^{'}, θ, λ) L (z, θ^{'}, λ) sin θ^{'} d θ^{'},

where

β (z, θ^{'}, θ, λ) = \int_{0}^{2 π} β (z, θ^{'}, ϕ^{'}, θ, ϕ = 0, λ) d ϕ^{'}

is the azimuthally integrated VSF. Equation (6) yields the azimuthally integrated radiance,

L (z, θ, λ) = \int_{0}^{2 π} L (z, θ, ϕ, λ) d ϕ,

which then can be integrated over polar angle θ to obtain the irradiances. Equation (6) is the form of the RTE solved by EcoLight-S (and also by the standard version of EcoLight that is included with HydroLight, except that the standard EcoLight retains the source term for inelastic effects). The azimuthally integrated RTE yields the same nadir and zenith radiances as HydroLight because in both cases the radiances for θ = 0 and π have no azimuthal dependence. EcoLight-S thus computes the same nadir-viewing water-leaving radiance, remote-sensing reflectance, and in-water zenith and nadir radiances as HydroLight. These are the radiometric quantities most often used, along with the plane irradiances, in inverse models.

Table 1. Effect of Layer Thickness on PAR Values at 15 m for the Environmental Conditions of Fig. 1 ^*

View Table | View all tables in this article

Table 2. Simulations of Pure Water and Turbid Case 2 Water

View Table | View all tables in this article

In HydroLight, the set of all (θ, ϕ) directions is discretized into M θ bins and N ϕ bins. The run time is proportional to (MN)² because radiance can be scattered from any one bin into any other. Removing the azimuthal dependence thus reduces the run time by a factor of N ². Compared to the standard HydroLight angular resolution of 15 deg in ϕ, i.e. N = 24, this is a reduction of (24)² in run time. The mathematical details of directional discretization are given in [1] for HydroLight and in [4] for EcoLight-S. Both codes use invariant imbedding techniques to solve the RTE.

It should be noted that EcoLight-S simply solves Eq. (6) for whatever inputs it is given, regardless of what physical environment those inputs represent. The RTE solution code is not, for example, restricted to a particular wavelength range. The association of inputs with wavelengths takes place in the user’s own code before calling EcoLight-S.

2.1. Depth specification

Most ecosystem models treat the water column as a stack of homogeneous layers, and EcoLight-S does the same. The IOPs are therefore constant with depth within each layer. This is different from HydroLight, which can model arbitrary depth variation in the IOPs. Modeling the water column as homogeneous layers makes it easy to match the EcoLight-S input and output to the depth structure of most ecosystem models, and it also allows the EcoLight-S code to run faster. As HydroLight solves the RTE by integrating Riccati equations as functions of depth (Ref. [1], Chapter 8), the RTE solver repeatedly calls a subroutine that computes and returns the IOPs at the current depth and wavelength. This can be computationally expensive if the IOP subroutine must interpolate between IOPs specified at discrete depths or if bio-optical models must be called to compute the absorption and scattering properties from other inputs. EcoLight-S on the other hand receives its IOP input as 2-D arrays of absorption, scattering, and backscatter coefficients, with the two array indices labeling the depth layer and wavelength band. Thus EcoLight-S needs only index these arrays for the current depth layer and wavelength band to obtain the IOPs. Computing the array depth index corresponding to any physical depth is computationally very fast. The user can define the layer depths by giving either the layer boundary depths or the layer midpoint depths as inputs to EcoLight-S.

Finite-depth bottom boundaries are treated as Lambertian reflectors with an irradiance reflectance specified by the user. If the water is infinitely deep, EcoLight-S computes a non-Lambertian bottom boundary condition corresponding to the IOPs of the bottom-most layer of the user’s ecosystem model. The appropriate boundary condition is then applied at the maximum depth where the RTE it so be solved.

Computed in-water irradiances, radiances, and other quantities are returned at both the layer boundaries and layer mid-points. Diffuse attenuation functions are returned as layer-averaged values.

2.2. IOP specification

The IOPs are specified by the total absorption, scattering, and backscattering coefficients, which in general are functions of depth and wavelength. Only the total IOPs—the sum of all contributions by water, phytoplankton, dissolved matter, mineral particles, etc.—are required to solve the RTE. The user can define the total IOPs in any desired manner, e.g., using the concentrations of the various constituents of the user’s biological model and bio-geo-optical models to convert concentrations to IOPs. Those calculations are done by the user before calling the EcoLight-S subroutine.

EcoLight-S uses the backscatter fraction, the ratio of the total backscatter coefficient to the total scattering coefficient, to determine a depth- and wavelength-dependent Fournier-Forand [5] scattering phase function (the ratio of the VSF to the scattering coefficient) in the manner described in [6]. The same technique is also used by HydroLight as one option for defining phase functions. This parameterization of Fournier-Forand phase functions by the backscatter fraction is admittedly crude. However, it works well for applications such as remote sensing because, to first order, it is the backscatter fraction that is most important in determining quantities such as upwelling radiance and remote-sensing reflectance. The exact shape of the phase function is less important. However, as shown in [6], better results are obtained if the shape of the phase function is determined as accurately as possible. Freda and Piskozub [7] developed a two-parameter technique that uses both the absorption and backscatter coefficients to select a Fournier-Forand phase function. Their technique is based on limited data, but it appears to work well and warrants further investigation for use in codes such as HydroLight and EcoLight-S.

In all be the very clearest waters the contribution of scattering by water itself to the total scattering is negligible. Therefore, the use of a Fournier-Forand phase function, which is based on Mie calculations for particles, is justified for use as the phase function for the total of water plus particles. The errors caused by this approximation are small compared to those resulting from the other approximations used in optimizing EcoLight-S for speed.

3. RTE solution options

EcoLight-S takes the following philosophy. It is necessary to solve the RTE in order to incorporate the effects of the sea-surface boundary and to account for all IOP effects. However, once an accurate value of the scalar irradiance E _o(z,λ) has been computed to some depth z _o deep enough to be free of surface boundary effects, it is not necessary to continue solving the RTE to greater depths, which is computationally expensive. As shown below, in many cases of practical interest it is possible to extrapolate the accurately computed upper-water-column irradiances to greater depths and still obtain irradiances that are acceptably accurate for ecosystem predictions. Likewise, it may not be necessary to solve the RTE at every wavelength used by a particular ecosystem model in order to obtain acceptably accurate irradiances at the needed wavelength resolution. In models requiring high wavelength resolution (e.g., the EcoSim biological model [8, 9], which requires E _o(z,λ) at 5 nm resolution from 400 to 700 nm.), omitting every other wavelength cuts the EcoLight-S run time by roughly one half. Models using only PAR can solve the RTE at relatively few wavelengths and still obtain PAR profiles accurate to within a few percent.

3.1. Dynamic depth solutions

If the option to use dynamic-depth solutions is not selected, then EcoLight-S solves the RTE to the greatest user-defined depth, regardless of the wavelength (just as HydroLight does). This is computationally expensive if the bottom is optically very deep, which is always the case at wavelengths greater than 700 nm for bottoms more than a few meters deep, owing to water absorption.

If the option to use dynamic-depth solutions is selected, EcoLight-S dynamically (i.e., as the program runs) determines the depth to which the RTE is solved at each wavelength. The goal is to solve the RTE to the shallowest depth possible at each wavelength, and then to extrapolate the scalar irradiance to greater depths. To determine the depth z _o to which the RTE will be solved at a particular wavelength, note that the scalar irradiance can be written as

E_{o} (z_{o}, λ) = E_{o} (0, λ) exp [- \int_{0}^{z_{o}} K_{o} (z, λ) d z] .

This equation defines K _o, the diffuse attenuation coefficient for scalar irradiance. Note that K _o is not known until after the RTE has been solved. However, except very near the sea surface where boundary effects are important, K _o is approximately equal to K _d, the diffuse attenuation coefficient for downwelling plane irradiance (K _o becomes exactly equal to K _d at great depths in homogeneous water). To first order (e.g., Ref [1] Eq. (5.65)) K _d ≈ a/ $\bar{μ}$ _d, where a is the absorption coefficient and $\bar{μ}$ _d is the mean cosine of the downwelling radiance distribution. In typical waters, $\bar{μ}$ _d ≈ 3/4. Thus, in Eq. (7) K _o can be roughly approximated by the absorption coefficient a(z,λ). Equation (7) can then be rewritten as

F_{o} = \frac{E_{o} (z_{o}, λ)}{E_{o} (0, λ)} \approx exp [- \int_{0}^{z_{o}} a (z, λ) d z] .

Here F _o is the fraction of the surface scalar irradiance that penetrates to depth z _o at the given wavelength.

Equation (8) can be used to estimate the depth to which the RTE will be solved because the absorption coefficient is an input to EcoLight-S and is known before the RTE is solved. The user selects a value of F _o, 0 ≤ F _o ≤ 1, as one of the inputs to EcoLight-S. Equation (8) is then solved for the depth z _o corresponding to the chosen F _o and the known absorption coefficient. In practice, because of the use of homogeneous layers, this amounts to recursively computing f _o = 1 and then f_k = f _k−1 exp[−a(k, λ) Δz_k], k = 1, 2, 3,... until f_k < F _o. Here a(k, λ) is the absorption coefficient for depth layer k, which has thickness Δz_k. The next deeper layer mid-point or boundary depth is then taken to be z _o. Thus the RTE is always solved to a grid output depth greater than or equal to the actual depth corresponding to F _o and the IOPs. Moreover, this estimate of z _o will always be greater than the actual z _o value because the mean cosine factor is omitted. Omitting the mean cosine is equivalent to having too little absorption; hence the solution will go too deep. After the RTE is solved at the first wavelength, the actual value of z _o can be obtained from the computed irradiance. This value is then used to adjust the estimated z _o at the next wavelength, and so on for subsequent wavelenths. After the first wavelength, which always goes too deep, this algorithm results in z _o estimates that are very close to the actual F _o depths (Fig. 1 below).

Fig. 1 Example simulation showing the use of dynamic solution depths for RTE solutions down to the depth where the scalar irradiance has decreased to 10% of the surface value.

Download Full Size | PDF

Using F _o = 0.1 would result in solving the RTE to the 10% irradiance level at each wavelength, i.e., to the depth where the irradiance has decreased to 0.1 or 10% of its value at the sea surface. There is also an option in the code to let the user-selected value of F _o be the spectral scalar irradiance E _o(z _o, λ) in W m⁻² nm⁻¹. With that option, using F _o = 0.1 would result in solving the RTE to a depth where the irradiance has decreased to 0.1 W m⁻² nm⁻¹, regardless of the surface value. That option is more reasonable for calculations of photosynthesis because it is irradiance magnitudes, not percentages, that determine photosynthesis. Regardless of which option is used to estimate depth z _o, the bottom boundary condition is then applied at the next layer midpoint or boundary depth z_k below z _o, and the RTE is solved only between the surface and z_k. The resulting radiances and derived quantities are then accurate down to depth z_k.

3.2. Extrapolation below the solution depth

After the RTE has been solved to some depth z_k at a particular wavelength, the computed values of E _o(z_k) and E _d(z_k) are extrapolated to greater depths as follows. The extrapolation is based on Eq. (7), except that the mean cosine factors can now be included in the approximations for K _d ≈ a/ $\bar{μ}$ _d and K _o ≈ a/ $\bar{μ}$ . In addition to E _o(z_k) and E _d(z_k), the solution of the RTE gives all irradiances, in particular the downwelling scalar irradiance E _od(z_k) and the upwelling plane irradiance E _u(z_k). These irradiances are used to compute the mean cosines

{\bar{μ}}_{d} (z_{k}) = \frac{E_{d} (z_{k})}{E_{od} (z_{k})} and \bar{μ} (z_{k}) = \frac{E_{d} (z_{k}) - E_{u} (z_{k})}{E_{o} (z_{k})}

at the last solved depth. These values of $\bar{μ}$ _d(z_k) and $\bar{μ}$ (z_k) are then used at all lower depths. Thus Eq. (7) becomes

E_{o} (z, λ) = E_{o} (z_{k}, λ) exp [- \int_{z_{k}}^{z} \frac{a (z, λ)}{\bar{μ} (z_{k}, λ)} d z] .

A similar equation using $\bar{μ}$ _d(z_k) holds for E _d. Equation (10) is then applied (in layer summation form) to the homogeneous layers, beginning at the last computed depth z_k and extending to the maximum depth. The upwelling quantities E _u and L _u are not extrapolated below the maximum solution depth because they can be strongly influenced by bottom reflectance in shallow waters. In simulations with optically shallow bottoms the maximum solution depth determined by the absorption coefficient is usually deeper at visible wavelengths than the actual bottom depth, in which case all quantities are accurately computed down to the physical bottom.

Note that E _o(z_k,λ) and the other irradiances incorporate all of the effects of the surface boundary and of the water IOPs above the maximum depth z_k to which the RTE was solved. The extrapolations based on Eqs. (9) and (10) will be reasonably accurate if the variability in the mean cosines is not great below depth z_k and if the IOPs covary with the absorption coefficient. This is often a good approximation, but might not be the case, for example, if there were a layer of highly scattering but non-absorbing particles below depth z_k. In such a case, it would be necessary to solve the RTE to a depth deeper than the scattering layer if high accuracy is required for the computed irradiances below the scattering layer.

3.3. Wavelength skipping

Some ecosystem models require the scalar irradiance at many wavelengths; others require only PAR. In either case, it is often sufficient to solve the RTE at relatively few wavelengths and then obtain the irradiances at other wavelengths by interpolation between the computed wavelengths. For example, in the EcoSim model [8,9], which requires 5 nm resolution, it might be sufficient to solve the RTE at every fourth EcoSim wavelength band (i.e. in 5 nm wide bands spaced 20 nm apart) and then interpolate to obtain the irradiances at the required 5 nm resolution. This would cut the run time by roughly a factor of four, because each wavelength solution requires approximately the same run time (when using dynamic depth solutions).

If the wavelength-skipping option is not chosen, then the RTE is solved at each of the user’s input wavelengths. If wavelength skipping is chosen, then the user inputs the number of wavelengths to skip, nwskip, between RTE solutions. The RTE is always solved at the first and last user wavelengths. Thus if the user’s input wavelengths (where the IOPs and other inputs are defined) are 400 to 700 nm by 10 nm, and nwskip = 0, the RTE will be solved for 10 nm wide bands centered at 405, 415, ...685, 695 nm. If nwskip = 1, the RTE will be solved at 405, 425, ..., 665, 685, and 695 nm. Values at 415, 435, ..., 675 nm will be computed by linear interpolation at each depth between the computed values. If nwskip = 2, the RTE will be solved at 405, 435, ..., 665, and 695 nm, with values at 415 and 425 being obtained by interpolation between the computed 405 and 435 nm values, and so on.

3.4. Other optimizations

The dynamic depth and wavelength skipping options can be used separately or together. The run-time minimization that can be achieved by use of dynamic depths and wavelength skipping in a particular application of EcoLight-S depends on the level of accuracy of the computed PAR or irradiances required by the user’s particular application. There are no simple guidelines as to what value of F _o or what wavelength resolution are adequate because different applications have different accuracy requirements, and the acceptable optimizations for a given accuracy may depend on water IOPs, bottom conditions, and the like. Some experimentation may be necessary to determine the acceptable amount of depth and wavelength optimization for a given application, before production runs are begun.

Just as with HydroLight, EcoLight-S can compute the asymptotic K functions and reflectances using the IOPs at the deepest input depth. However, those computations add to the run time and can be omitted if the asymptotic values are not of interest. Likewise, the computation of PAR can be omitted if the user’s application does not require PAR.

4. Example simulations

A few example simulations suffice to show the possible decreases in run time and resulting decreases in accuracy of the computed variables for various optimizations. The timing was done on a computer with a 2.4 GHz, 32 bit, Intel Core i5 CPU, 4 Gbytes of RAM, and Microsoft Windows 7 operating system. Since the main anticipated application of EcoLight-S is computation of PAR in ecosystem models, PAR will be used as the variable of primary interest for comparing optimizations and run times.

In the first example, Case 1 water was modeled using a background chlorophyll value of Chl = 0.5 mg m⁻³ plus a gaussian profile with its maximum value of 2.0 mg m⁻³ at 15 m depth. This continuous Chl(z) profile is shown by the red line in the left panel of Fig. 1. The EcoLight-S water column was modeled as 5 m thick homogeneous layers, which is typical of the depth resolution used in ocean ecosystem models. The layer chlorophyll values are the average of the continuous profile within each layer. The layer chlorophyll values were converted to layer absorption, scatter, and backscatter coefficients using a bio-optical model [10]. The sun was at a 30 deg zenith angle in a clear sky; the wind speed was 5 m s⁻¹. The bottom boundary condition was for infinitely deep water below 50 m. The dynamic depth option was selected with F _o = 0.1, interpreted as a percent of the surface irradiance at each wavelength. The greatest depth requested for output was 50 m. The solution wavelengths were 5 nm bands from 400 to 700 nm; no wavelength skipping was done.

The right panel of Fig. 1 shows the estimated and actual, physical and optical, depths corresponding to the F _o = 0.1 or 10% irradiance level at each wavelength. Note that for the first wavelength band (400-405 nm, plotted at 402.5 nm) the RTE was solved considerably deeper than necessary. The estimated solution depths at subsequent wavelengths were only slightly deeper than the actual 10% irradiance depths as determined after the RTE was solved. The five-meter layer thicknesses make it easy to see (from the red dots) that the RTE was always solved to either a layer boundary or layer midpoint depth (2.5 m increments), where the output was saved. The main computational savings in this run come at wavelengths greater than 600 nm, where the optical depth of the bottom at 50 m increases from 30 (at 600 nm) to 48 (at 700 nm), but the solution optical depths decrease from 8 to 4. The corresponding run without dynamic depths, for which the RTE was solved to 50 m at each wavelength, took 1.42 s, but the run with dynamic depths and F _o = 0.1 took only 0.53 s. The two PAR profiles agreed to within 0.4% at all depths down to 50 m.

Figure 2 shows the errors in PAR profiles for the same bio-physical and environmental conditions as Fig. 1 but for runs with different F _o values and no wavelength skipping. The inset labels in the second panel from the left show the values of F _o, the wavelength resolution (5 nm in all cases, corresponding to nwskip = 0), and the run times in seconds. The label colors correspond to the curve colors in all panels. The third panel shows the relative error in percent,

relative error = 100 \frac{P A R (optimized) - P A R (unoptimized)}{P A R (unoptimized)},

of the optimized PAR profiles compared to the unoptimized profile, and the last panel shows the actual error = PAR(optimized) – PAR(unoptimized). The colored dots on the PAR profiles in the second panel show the maximum depth to which the RTE was solved for the various F _o values, after the initial solution. These runs show for example that, for this chlorophyll profile, the RTE can be solved to only the 20% irradiance depth at each wavelength (green curves) with resulting errors in PAR of less than 3% or 1 μmol photons m⁻² s⁻¹ down to 50 m, which is roughly the depth at which PAR has decreased to 1% of its surface value. Solving the RTE to just the 50% irradiance depth (red curves), which is above the depth of the cholorphyll maximum, still gives PAR values within 9% or 7 μmol photons m⁻² s⁻¹. The optimized run times are all less than 1 second.

Fig. 2 Errors in PAR for different dynamic solution depths for the same environmental conditions as Fig. 1. The color-coded inset in the second panel shows the value of F _o, the wavelength resolution (always 5 nm for this figure, corresponding to nwskip = 0), and the run time in seconds. The dots along the PAR profiles show the greatest RTE solution depth (after the solution for the first wavelength) for each of the simulations, e.g., 32.5 m for the simulation of Fig. 1, which is the purple curve in this figure. The third panel shows the relative errors in PAR compared to the unoptimized run, and the rightmost panel shows the actual PAR errors.

Download Full Size | PDF

The compensation depth, where phytoplankton photosynthesis equals respiration, depends on phytoplankton species and other factors such as temperature. A typical value is about 5 μmol photons m⁻² s⁻¹ for microplankton. The common rule of thumb for the depth of the euphotic zone is the depth where PAR has decreased to 1% of its mid-day surface value, which is about 50 m for the water properties of Fig. 2. It is reasonable to assume that errors in PAR below the depth where PAR is roughly 10 μmol photons m⁻² s⁻¹ will not greatly affect primary production calculations in ecosystem models. In the simulation of Fig. 2 that depth is about 60 m.

Figure 3 shows the results for wavelength-skipping optimization, with the RTE being solved to 50 m and the same environmental conditions as the simulation of Fig. 1. The labels in the second panel show the wavelength resolution for solving the RTE: λ = 10 (nwskip = 1) means that the RTE was solved at every other one of the 5 nm bands; λ = 50 (nwskip = 9) corresponds to RTE solutions only at every tenth band, centered at 402.5, 452.5, ..., 697.5 nm. Even for just 50 nm resolution, PAR is computed within 3% of the unoptimized 5 nm values. The largest magnitude errors are now near the sea surface, with the largest error of 33 μmol photons m⁻² s⁻¹ (2885 vs. 2852 μmol photons m⁻² s⁻¹) being at the surface for the simulation at 20 nm resolution.

Fig. 3 Errors in PAR for different wavelength resolutions (different values of nwskip) compared to 5 nm bands from 400 to 700 nm, for the same environmental conditions as Fig. 1.

Download Full Size | PDF

A run with an unoptimized wavelength resolution of 5 nm and nwskip = 3 solves the RTE at 20 nm resolution. Strictly speaking this is not the same radiative transfer problem as using a 10 nm resolution and skipping one wavelength between solutions, or using a 20 nm resolution with no wavelength skipping. However, the results and run times for these cases are very similar. For these three cases and the current chlorophyll profile, the PAR profiles are the same to within ±2% and the run times are the same to within a few hundredths of a second.

Dynamic-depth and wavelength-skipping optimizations can of course be combined. Figure 4 shows the PAR errors for various combinations of F _o and wavelength resolution for the same conditions as Fig. 1. Solving the RTE to only the 50% irradiance level at 25 nm resolution still gives PAR errors of less than 10% down to 50 m with a run time of only 0.05 s. The less-optimized runs all have PAR errors of less than 4% and run times of 0.2 s or less.

Fig. 4 Errors in PAR for combinations of solution depths (F _o values) and wavelength skipping (nwskip values; the labels show the resulting wavelength resolution) compared to solutions to 50 m for 5 nm bands from 400 to 700 nm. The environmental conditions are the same as for Fig. 1.

Download Full Size | PDF

Figure 5 shows the scalar irradiances E _o(z, λ) for the unoptimized solution of Fig. 4, and for the optimized solution with F _o = 0.2 and 25 nm resolution (nwskip = 4), which is the orange curve of that figure. The black squares in the upper right panel show the wavelengths and F _o depths to which the RTE was solved in the optimized run. Visual comparison of the contour plots of the upper panels of Fig. 5 show that the unoptimized and optimized irradiances are similar, although the effects of interpolation between the solved wavelengths are apparent. In particular, the scalloped appearance of the contours in the upper right panel at red wavelengths and large depths occurs because the linear interpolation between the solution wavelengths at a given depth does not accurately describe the non-linear change with wavelength in E _o(z, λ), which results from the rapid increase in water absorption with wavelength. The lower left panel shows the relative errors in percent computed as in Eq. (11). Most of the water column between 400 and 600 nm has errors of less than ±10% (blue colors), but the relative errors are large for the unsolved wavelengths in the red spectral region at large depths (red colors). These large relative errors occur only where the irradiances are several order of magnitude less than the maximum irradiances at the same depth. The small irradiances with the large relative errors make little contribution to the PAR value. The lower right panel shows that the errors in the irradiance magnitudes themselves are usually less than ±0.02, and almost always less than ±0.05 W m⁻² nm⁻¹.

Fig. 5 Scalar irradiances for unoptimized (E _o(un), upper left panel) and optimized (E _o(op), upper right panel) solutions. The optimized solution corresponds to the orange curve of Fig. 4. The lower left panel shows the relative errors E _o computed as in Eq. (11), and the lower right panel shows the actual errors in E _o.

Download Full Size | PDF

Figures 1–5 have all been for 5 m thick layers, which were chosen as typical of the resolution in ecosystem models. However, the depth resolution affects the computed PAR values, just as do the optimizations of solution depth and wavelength resolution. Figure 6 repeats the simulations of Fig. 4 for a one-meter depth resolution. Comparison of Figs. 4 and 6 shows that the errors in PAR relative to the unoptimized values are determined mostly by the solution depth and wavelength optimizations. What is not easily seen in the scale of the PAR(z) profiles of those figures is that the PAR profiles depend on the depth resolution of the IOPs in regions where the water is inhomogeneous. Table 1 shows the values of PAR at 15 m, the depth of the chlorophyll maximum, for the unoptimized and optimized red curves, which have F _o = 0.5 and λ = 25 (nwskip = 4), of Figs. 4 and 6. These values illustrate that for either depth resolution the relative and actual errors are almost the same, but that the PAR values are different for the different layer thicknesses. For this example, the PAR values are about 8% less for the 5 m layer thickness than for the 1 m layers, and these differences are almost the same for the unoptimized and optimized solutions.

Fig. 6 Errors in PAR for the same simulations as Fig. 4, but with a one-meter depth resolution.

Download Full Size | PDF

Thinner layers can be used to get better depth resolution of the continuous chlorophyll profile seen in Fig. 1. A layer thickness of 0.5 m, corresponding for example to data from a profiling instrument with measurements binned into 0.5 m depth bins, gives a good approximation to the continuous chlorophyll profile. Figure 7 compares this high-depth-resolution case with layer thickness of 1, 2, and 5 m. The inset labels now color code the curves by layer thickness Δz and run time. These runs all gave the same PAR profiles to within 1% near the surface and below 25 m. However, near the chlorophyll maximum the difference between the high-resolution (Δz = 0.5 m) and the low-resolution (5 m) cases was as much as 11% or 65 μmol photons m⁻² s⁻¹. The unoptimized run times increased from 1.42 s for the 5 m layers to 2.09 s for the 0.5 m layers. The output files for 0.5 m layers are ten times as large as the files for 5 m layers because output is returned at all layer boundaries and mid-points. Corresponding runs optimized with F _o = 0.2 and 25 nm wavelength resolution (as shown in Table 1) are the same as the unoptimized runs to within 2% down to 30 m, and differ by less than 4% at 50 m. The optimized run times were 0.09 s (5 m layers) to 0.16 s (0.5 m layers). It should be remembered that different depth resolutions of the continuous profile are different radiative transfer problems even though the column-integrated chlorophyll values are the same. Thus the differences seen in Table 1 and Fig. 7 are not errors in the same sense as those due to solution depth and wavelength optimizations, they simply show that different IOP profiles have different PAR profiles.

Fig. 7 Differences in PAR for different depth resolutions of the continuous chlorophyll profile of Fig. 1, for unoptimized solutions.

Download Full Size | PDF

These results show that the accuracy of computed PAR profiles depends on a complicated interplay of depth resolution, solution depth, and wavelength resolution. In many EcoLight-S applications, the layer thicknesses will be determined by the ecosystem model. However, if high depth resolution can be used, then more accurate PAR profiles can be computed. If PAR profiles within 5% accuracy of the continuous profile are required, then these can be obtained for the present chlorophyll profile with layer thicknesses of 2 m, with little increase in run time compared to the 5 m layers. Very sharp gradients in IOP profiles could require higher depth resolution for a given accuracy of PAR, but nearly homogeneous waters can be well modeled with even thicker depth layers. EcoLight-S allows for unequal layer thicknesses.

Table 2 compares standard HydroLight and EcoLight version 5.1 and unoptimized and optimized EcoLight-S runs for simulations of pure sea water and homogeneous, turbid Case 2 water. The Case 2 water simulates phytoplankton with a chlorophyll concentration of 2 mg m⁻³, CDOM absorption of a _CDOM(440) = 0.2 m⁻¹, and “brown earth” minerals with a concentration of 1 gm m⁻³. These values were converted to IOPs using HydroLight’s generic bio-geo-optical model for Case 2 water. The backscatter fraction by phytoplankton was taken to be 0.01 at all wavelengths. The backscatter fraction of the minerals was wavelength dependent and given by 0.03(550/λ)^0.5. These concentrations and spectral functions give roughly equal contributions to absorption by the phytoplankton and minerals, with CDOM absorption being greater than either at blue wavelengths. Scattering by phytoplankton and minerals is roughly equal, but the backscatter coefficient for minerals is 2-3 times that of the phytoplankton. CDOM was assumed to be nonscattering. The sun was placed at the zenith in a clear sky; the resulting above-surface PAR value was 2882 μmol photons m⁻² s⁻¹. The sea surface was level. The HydroLight and EcoLight runs were made from 350–700 nm with 10 nm bandwidths so that inelastic scatter effects (Raman scatter by the water and chlorophyll and CDOM fluorescence) from below 400 nm would be accounted for in the PAR wavelengths of 400-700 nm. The EcoLight-S unoptimized runs were from 400 to 700 nm by 10 nm. The pure water runs were made to a depth of 400 m, and the Case 2 runs to a depth of 20 m.

The HydroLight simulations with inelastic scattering were taken to be the standards for comparison with the other runs. The deviations of PAR from this standard increase with depth for these homogeneous waters, and the second column of the table gives the PAR value at the greatest simulation depth, either 400 m or 20 m, where the difference was greatest. The first time in the third column for the HydroLight and EcoLight runs is the time for the 350-700 nm simulation; the time in parentheses is the time for a run from 400-700 nm. Comparing the 350-700 run times with and without inelastic scatter shows the cost of the extra computations required to include inelastic processes. The run times for 400-700 are directly comparable with the EcoLight-S run times. The last column shows the percent differences in the PAR values of column 2 between the HydroLight run with inelastic effects and the other runs.

The large decreases in run times between corresponding HydroLight 5.1 and EcoLight 5.1 runs are due to the difference in solving the full RTE of Eq. (1) vs. the azimuthally integrated RTE of Eq. (6). The differences in run times between HydroLight or EcoLight runs with and without inelastic effects show the computational costs of including inelastic scatter in the solution of the RTE. The inelastic scatter calculations are a larger fraction of the total run time for EcoLight than for HydroLight. The timing difference between EcoLight 5.1 without inelastic effects and the unoptimized EcoLight-S is due to the greater efficiency of the EcoLight-S code, in particular obtaining the IOPs from indexed arrays rather than from subroutine calls.

The differences in the various EcoLight-S optimizations show that there is no simple rule for determining how much error in computed PAR values results from a particular optimization. For the pure water simulations, increasing F _o from 0.1 to 0.2 (solving the RTE to a shallower depth) for a given wavelength resolution has little effect on the PAR difference at 400 m, but decreasing the wavelength resolution from 10 to 20 nm for a given F _o significantly increases the error. However, the oppose occurs with the Case 2 water: changing the wavelength resolution from 10 to 20 nm has little effect, but changing F _o from 0.1 to 0.2 triples the error at 20 m. However, regardless of which EcoLight-S optimization is used, these runs show that even in extreme cases of pure water or very turbid water, it is possible to compute PAR values to within roughly 10% in a few tenths of a second of computer time. The optimized EcoLight-S runs are usually more than 1,000 times faster than the HydroLight run with inelastic effects.

In addition to the in-water irradiances needed by ecosystem models, EcoLight-S also computes other useful quantities, including the nadir-viewing remote-sensing reflectance R _rs. Figure 8 shows the R _rs spectra for the pure and Case 2 water simulations of Table 2. The solid lines are the spectra computed by HydroLight 5.1 including inelastic effects. The red open circles are the unoptimized EcoLight-S values, and the black dots are the EcoLight-S values with F _o = 0.2 and 10 nm resolution. The EcoLight-S values are always less than the HydroLight values because of the omission of inelastic effects, but the difference is noticeable only in extreme cases such as the pure water simulation, for which Raman scatter makes its maximum contribution to the water-leaving radiance. The contributions of Raman scatter and chlorophyll and CDOM fluorescence to R _rs are insignificant in the Case 2 water case except in the chlorophyll fluorescence band at 685 nm where there is a small fluorescence peak (the HydroLight runs assumed a 2% quantum efficiency for chlorophyll fluorescence). The optimized and unoptimized EcoLight-S runs are indistinguishable in this figure because R _rs is determined by the near-surface light field, which is being accurately computed in all cases, regardless of the effect of the F _o value on the light field at depth.

Fig. 8 Remote-sensing reflectances R _rs for the pure water and Case 2 runs of Table 2. The solid curves are the HydroLight 5.1 runs with inelastic scattering. The open circles are the unoptimized EcoLight-S runs, and the filled dots are EcoLight-S with F _o = 0.2.

Download Full Size | PDF

5. Discussion

Only a few representative simulations have been shown here, but many others have been performed in the course of EcoLight-S development. These include a variety of chlorophyll concentrations and profiles for Case 1 waters and inhomogeneous Case 2 waters with various combinations of non-covarying chlorophyll, CDOM, and mineral particles. These all indicate that regardless of the water type or complexity, EcoLight-S can compute PAR to within 10% accuracy throughout the euphotic zone in a few tenths of a second of time on a moderately fast computer.

An initial (and slower) version of EcoLight-S was used [11] as the optical component of a coupled ROMS (Regional Ocean Modeling System [12])–EcoSim ( [8, 9]) ecosystem model and compared with the EcoSim analytical model for spectral scalar irradiance. In idealized five-year simulations of ecosystem annual cycles in Case 1 water, EcoLight-S was able to compute the spectral scalar irradiance E _o(z,λ) at the 5 nm resolution needed by EcoSim with less than a 30% increase in the total simulation time. EcoLight-S is currently being used as the optical component of a coupled ROMS-CoSiNE (Carbon Silicon Nitrogen Ecosystem [13, 14]) model for more complex simulations of ocean upwelling systems. The standard ROMS-CoSiNE model uses one analytical model to obtain PAR for biological primary production calculations within CoSiNE and another model based on short-wave irradiances and an assumed Jerlov water type for heating calculations within the ROMS hydrodynamic model. Initial simulations show [15] that EcoLight-S runs over 400-1000 nm can replace both analytical light models in ROMS-CoSiNE with less than a 20% increase in total ecosystem simulation times. Including the 700-1000 nm wavelength range roughly doubles the run times for optimized EcoLight-S solutions. This is because the high absorption by water beyond 700 nm quickly attenuates the irradiance with depth, and the dynamic-depth solutions therefore stop at very shallow depths. The run times for 400-1000 nm are still a fraction of a second.

EcoLight-S has also been used [16] as the RTE forward model within an implicit inverse model for recovery of IOPs from in-water measurements of E _d and L _u obtained from a glider. The effect of errors in the computed E _d and L _u, i.e., the acceptable level of EcoLight-S optimization, on the convergence of implicit inverse models and on the accuracy of their retrievals must be evaluated by the user for each particular inversion algorithm. Most users of EcoLight-S in an inversion model would probably choose to solve the RTE at the wavelengths and depths corresponding to the observational data used by the inversion algorithm, in which case the errors in the computed light quantities would be minimal.

Unlike simple analytical light models for PAR or spectral irradiance, EcoLight-S can account for depth- and wavelength-dependent IOPs, variable sky conditions, sea surface boundary effects, and reflectance by shallow bottoms. It can simulate any water body from pure water to the most complex of Case 2 waters, for which there are no analytic light models. EcoLight-S also computes optical quantities such as the nadir-viewing remote-sensing reflectance, in-water upwelling radiance, plane irradiances, and diffuse attenuation functions corresponding to the bio-optical state of the ecosystem. The remote-sensing reflectance allows for validation of ecosystem model predictions using satellite ocean color radiometry without an intervening step to convert a satellite-measured radiance to a chlorophyll concentration via an imperfect chlorophyll algorithm. The in-water quantities enable validation from optical measurements made by moorings, gliders, or autonomous underwater vehicles. The availability of these ancillary quantities alone argues for the use of EcoLight-S as the optical component of ecosystem models, even for situations for which analytical PAR or spectral irradiance models give satisfactory results.

Acknowledgments

The development of EcoLight-S was supported by the Environmental Optics (now Littoral Sciences and Optics) Program of the U. S. Office of Naval Research, and by HydroLight revenues. The author thanks Jacek Piskozub and an anonymous reviewer for their helpful comments.

References and links

1. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic, 1994), http://www.curtismobley.com/lightandwater.zip.

2. C. D. Mobley and L. K. Sundman, HydroLight User’s Guide (Sequoia Scientific, Inc., 2008), http://www.hydrolight.info.

3. C. D. Mobley and L. K. Sundman, HydroLight Technical Documentation (Sequoia Scientific, Inc., 2008), http://www.hydrolight.info.

4. C. D. Mobley, EcoLight-S 1.0 User’s Guide and Technical Documentation (Sequoia Scientific, Inc., 2011), http://www.hydrolight.info. [PubMed]

5. G. R. Fournier and J. L. Forand, “Analytic phase function for ocean water,” in Ocean Optics XII, J. Jaffe, ed., Proc. SPIE2258 (with corrections), 194–201 (1994). [CrossRef]

6. C. D. Mobley, L. K. Sundaman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. 41(6), 1036–1050 (2002). [CrossRef]

7. W. Freda and J. Piskozub, “Improved method of Fournier-Forand marine phase function parameterization,” Opt. Express 15(20), 12763–21768 (2007). [CrossRef] [PubMed]

8. W. P. Bissett, J. J. Walsh, D. A. Dieterle, and K. L. Carder, “Carbon cycling in the upper waters of the Sargasso Sea: I. numerical simulation of differential carbon and nitrogen fluxes,” Deep-Sea Res. 46, 205–269 (1999a). [CrossRef]

9. W. P. Bissett, K. L. Carder, J. J. Walsh, and D. A. Dieterle, “Carbon cycling in the upper waters of the Sargasso Sea: II. numerical simulation of apparent and inherent optical properties,” Deep-Sea Res. 46, 271–317 (1999b). [CrossRef]

10. C. D. Mobley, “A new IOP model for Case 1 water,” in Ocean Optics Web Book, http://www.oceanopticsbook.info/view/optical_constituents_of_the_ocean/__level_2/a_new_iop_model_for_case_1_water.

11. C. D. Mobley, L. K. Sundman, W. P. Bissett, and B. Cahill, “Fast and accurate irradiance calculations for ecosystem models,” Biogeosci. Discuss. 6, 10625–10662 (2009), http://www.biogeosciences-discuss.net/6/10625/2009/bgd-6-10625-2009.pdf. [CrossRef]

12. A. F. Shchepetkin and J. C. McWilliams, “The regional ocean modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model,” Ocean Model. 9, 347–404 (2005), https://www.myroms.org/wiki/index.php/Documentation_Portal. [CrossRef]

13. F. Chai, R. C. Dugdale, T.-H. Peng, F. P. Wilkerson, and R. T. Barber, “One dimensional ecosystem model of the equatorial Pacific upwelling system, part I: model development and silicon and nitrogen cycle,” Deep-Sea Res. , Part II 49, (13–14), 2713–2745 (2002).

14. M. Fujii, E. Boss, and F. Chai, “The value of adding optics to ecosystem models: a case study,” Biogeosciences 4, 817–835 (2007), http://www.biogeosciences.net/4/817/2007/. [CrossRef]

15. F. Chai, “Incorporating optics into physical and ecosystem modeling,” presented at the Gordon Research Conference on Coastal Ocean Modeling, South Hadley, MA, 26 June–1 July, 2011.

16. E. Rehm, C. D. Mobley, and J. Smart, “Inverting light with constraints,” presented at Ocean Optics XIX, Barga, Italy, 6–10 Oct. 2008, http://staff.washington.edu/erehm/OOXIX-Rehm-final.pdf.

layer thickness	unoptimized PAR [μ mol m⁻² s⁻¹]	optimized PAR [μ mol m⁻² s⁻¹]	relative error [percent]	actual error [μ mol m⁻² s⁻¹]
1 m	366.08	359.77	−1.72	−6.31
5 m	335.75	330.22	−1.65	−5.53
1 m – 5 m	30.33	29.55	−2.57	−0.78

Pure water: z _max = 400 m; Secchi depth = 120 m

model	PAR(z _max) [μ mol m⁻² s⁻¹]	time [seconds]	difference [percent]
HydroLight 5.1 with inelastic	10.1330	427.1 (376.9)	—
HydroLight 5.1 without inelastic	9.6223	283.6 (268.4)	−5.0
EcoLight 5.1 with inelastic	10.1090	15.3 (12.8)	−0.2
EcoLight 5.1 without inelastic	9.1690	5.1 (4.6)	−5.1
EcoLight-S unoptimized	9.6578	2.71	−4.7
EcoLight-S with F _o = 0.1, 10 nm	9.6342	0.17	−4.9
EcoLight-S with F _o = 0.2, 10 nm	9.9650	0.13	−4.3
EcoLight-S with F _o = 0.1, 20 nm	8.7881	0.08	−13.3
EcoLight-S with F _o = 0.2, 20 nm	8.8332	0.07	−12.8

Case 2 water: z _max = 20 m; Secchi depth = 3.7 m

HydroLight 5.1 with inelastic	1.6128	198.6 (162.3)	—
HydroLight 5.1 without inelastic	1.6032	108.3 (89.4)	−0.6
EcoLight 5.1 with inelastic	1.6068	6.7 (5.5)	−0.4
EcoLight 5.1 without inelastic	1.5968	2.3 (1.9)	−1.0
EcoLight-S unoptimized	1.5975	0.98	3.0
EcoLight-S F _o = 0.1, 10 nm	1.6606	0.28	3.0
EcoLight-S F _o = 0.2, 10 nm	1.7696	0.21	9.7
EcoLight-S F _o = 0.1, 20 nm	1.6594	0.15	2.9
EcoLight-S F _o = 0.2, 20 nm	1.7702	0.11	9.8

layer thickness	unoptimized PAR [μ mol m⁻² s⁻¹]	optimized PAR [μ mol m⁻² s⁻¹]	relative error [percent]	actual error [μ mol m⁻² s⁻¹]
1 m	366.08	359.77	−1.72	−6.31
5 m	335.75	330.22	−1.65	−5.53
1 m – 5 m	30.33	29.55	−2.57	−0.78

Pure water: z _max = 400 m; Secchi depth = 120 m

model	PAR(z _max) [μ mol m⁻² s⁻¹]	time [seconds]	difference [percent]
HydroLight 5.1 with inelastic	10.1330	427.1 (376.9)	—
HydroLight 5.1 without inelastic	9.6223	283.6 (268.4)	−5.0
EcoLight 5.1 with inelastic	10.1090	15.3 (12.8)	−0.2
EcoLight 5.1 without inelastic	9.1690	5.1 (4.6)	−5.1
EcoLight-S unoptimized	9.6578	2.71	−4.7
EcoLight-S with F _o = 0.1, 10 nm	9.6342	0.17	−4.9
EcoLight-S with F _o = 0.2, 10 nm	9.9650	0.13	−4.3
EcoLight-S with F _o = 0.1, 20 nm	8.7881	0.08	−13.3
EcoLight-S with F _o = 0.2, 20 nm	8.8332	0.07	−12.8

Case 2 water: z _max = 20 m; Secchi depth = 3.7 m

HydroLight 5.1 with inelastic	1.6128	198.6 (162.3)	—
HydroLight 5.1 without inelastic	1.6032	108.3 (89.4)	−0.6
EcoLight 5.1 with inelastic	1.6068	6.7 (5.5)	−0.4
EcoLight 5.1 without inelastic	1.5968	2.3 (1.9)	−1.0
EcoLight-S unoptimized	1.5975	0.98	3.0
EcoLight-S F _o = 0.1, 10 nm	1.6606	0.28	3.0
EcoLight-S F _o = 0.2, 10 nm	1.7696	0.21	9.7
EcoLight-S F _o = 0.1, 20 nm	1.6594	0.15	2.9
EcoLight-S F _o = 0.2, 20 nm	1.7702	0.11	9.8

Fast light calculations for ocean ecosystem and inverse models

Abstract

1. Introduction

2. RTE formulation

2.1. Depth specification

2.2. IOP specification

3. RTE solution options

3.1. Dynamic depth solutions

3.2. Extrapolation below the solution depth

3.3. Wavelength skipping

3.4. Other optimizations

4. Example simulations

5. Discussion

Acknowledgments

References and links

Cited By

Figures (8)

Tables (2)

Equations (13)

Optics Express