Radiative Transfer

BING implements radiative transfer models for calculating remote sensing reflectance (Rrs) from inherent optical properties (IOPs). This module provides both standard and advanced formulations of the Gordon semi-analytical model.

Overview

The radiative transfer module (bing.rt) converts optical properties into observable remote sensing reflectance:

\[R_{rs}(\lambda) = f(\lambda) \cdot \frac{b_b(\lambda)}{a(\lambda) + b_b(\lambda)}\]

where the function \(f(\lambda)\) captures the geometry of light propagation through the water column and depends on solar zenith angle, viewing geometry, and inherent optical properties.

Gordon Semi-Analytical Model

The Gordon (1988) model provides a semi-analytical relationship between subsurface remote sensing reflectance and IOPs:

\[r_{rs}(\lambda) = G_1(\lambda) \cdot u(\lambda) + G_2(\lambda) \cdot u^2(\lambda)\]

where:

\(u(\lambda) = \frac{b_b(\lambda)}{a(\lambda) + b_b(\lambda)}\) is the ratio of backscattering to total absorption+backscattering
\(G_1, G_2\) are Gordon coefficients that account for bidirectional light field geometry
\(r_{rs}\) is subsurface remote sensing reflectance

The subsurface reflectance \(r_{rs}\) is related to above-surface \(R_{rs}\) through:

\[R_{rs} = \frac{A \cdot r_{rs}}{1 - B \cdot r_{rs}}\]

where \(A = 0.52\) and \(B = 1.7\) are standard conversion coefficients.

Standard Gordon Coefficients

The classical Gordon model uses constant coefficients:

\[G_1 = 0.0949, \quad G_2 = 0.0794\]

These values were derived from radiative transfer simulations for typical oceanic conditions (clear sky, moderate sun angle).

from bing import rt as bing_rt
import numpy as np

# Define IOPs at measurement wavelengths
wavelengths = np.array([410, 440, 490, 510, 555, 670])  # nm

# Absorption coefficient [m^-1]
a = np.array([0.050, 0.045, 0.035, 0.038, 0.052, 0.350])

# Backscattering coefficient [m^-1]
bb = np.array([0.0025, 0.0020, 0.0015, 0.0013, 0.0010, 0.0006])

# Calculate Rrs with standard Gordon coefficients
Rrs = bing_rt.calc_Rrs(a, bb)

print(f"Rrs at 490 nm: {Rrs[2]:.6f} sr^-1")

Variable Gordon Coefficients

Recent research (Loisel et al. 2023) shows that Gordon coefficients vary with wavelength due to changes in scattering phase functions and light field geometry. BING provides wavelength-dependent Gordon coefficients derived from Hydrolight radiative transfer simulations.

The variable Gordon formulation uses interpolated coefficients:

\[G_1(\lambda) \text{ and } G_2(\lambda) \text{ are wavelength-dependent}\]

Key characteristics of variable Gordon coefficients:

G1 wavelength dependence: - Ranges from ~0.100-0.102 at blue wavelengths (350-450 nm) - Decreases to ~0.085-0.097 at red wavelengths (600-750 nm) - Generally positive across the spectrum
G2 wavelength dependence: - Positive at blue wavelengths (350-490 nm): ~0.05-0.07 - Transitions to negative at green-red wavelengths (>490 nm) - Reaches large negative values at far-red (>700 nm): -0.5 to -2.0 - Negative G2 values are physically valid and represent changes in light field structure

Note

The transition of G2 from positive to negative around 490-500 nm reflects changes in water absorption and the relative importance of single vs. multiple scattering events.

Loading Variable Gordon Coefficients

from bing.rt import rrs
import numpy as np

# Define model wavelengths
wavelengths = np.arange(400, 701, 5)  # 400-700 nm, 5 nm resolution

# Load wavelength-dependent Gordon coefficients
G1, G2 = rrs.wave_dependent_gordon(wavelengths)

# Calculate Rrs with variable Gordon coefficients
Rrs = bing_rt.calc_Rrs(a, bb, in_G1=G1, in_G2=G2)

# Examine coefficient variation
print(f"G1 range: {G1.min():.4f} - {G1.max():.4f}")
print(f"G2 range: {G2.min():.4f} - {G2.max():.4f}")
print(f"G2 transitions to negative at: {wavelengths[G2 < 0][0]:.0f} nm")

Gordon Coefficient Data Source

The variable Gordon coefficients are stored in bing/data/RT/gordon_coefficients.csv and include:

wavelength: 350-750 nm at 5 nm resolution
G1, G2: Fitted Gordon coefficients
G1_err, G2_err: Uncertainty estimates
rRMS, RMS: Fit quality metrics

Data derived from Loisel et al. (2023) Hydrolight synthetic dataset with ~20,000 bio-optical profiles representing diverse ocean conditions.

When to Use Variable Gordon Coefficients

Use standard Gordon (constant G1, G2) when:

Working with legacy algorithms or historical comparisons
Processing multispectral data (MODIS, SeaWiFS) with limited spectral information
Computational speed is critical
Results need to match standard products

Use variable Gordon (wavelength-dependent) when:

Processing hyperspectral data (PACE, SBG)
Maximum accuracy is required, especially at red/NIR wavelengths
Modeling oligotrophic (clear) waters where G2 effects are strongest
Developing new algorithms or validating models

Comparison of Standard vs. Variable Gordon

Spectral Differences

import matplotlib.pyplot as plt

# Compare standard and variable Gordon
wavelengths = np.arange(400, 751, 1)

# Variable Gordon coefficients
G1_var, G2_var = rrs.wave_dependent_gordon(wavelengths)

# Standard Gordon coefficients
G1_std = rrs.G1_STANDARD * np.ones_like(wavelengths)
G2_std = rrs.G2_STANDARD * np.ones_like(wavelengths)

# Plot comparison
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# G1 comparison
ax1.plot(wavelengths, G1_std, 'k--', label='Standard G1 (0.0949)', linewidth=2)
ax1.plot(wavelengths, G1_var, 'b-', label='Variable G1', linewidth=2)
ax1.set_ylabel('G1')
ax1.set_title('Gordon Coefficient G1')
ax1.legend()
ax1.grid(True, alpha=0.3)

# G2 comparison
ax2.plot(wavelengths, G2_std, 'k--', label='Standard G2 (0.0794)', linewidth=2)
ax2.plot(wavelengths, G2_var, 'r-', label='Variable G2', linewidth=2)
ax2.axhline(0, color='gray', linestyle=':', linewidth=1)
ax2.set_xlabel('Wavelength (nm)')
ax2.set_ylabel('G2')
ax2.set_title('Gordon Coefficient G2 (note negative values at λ > 490 nm)')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('gordon_coefficients_comparison.png', dpi=150)

Impact on Retrieved Rrs

# Calculate Rrs difference using both formulations
def compare_gordon_formulations(a, bb, wavelengths):
    """Compare Rrs calculated with standard vs. variable Gordon."""

    # Standard Gordon
    Rrs_std = bing_rt.calc_Rrs(a, bb)

    # Variable Gordon
    G1, G2 = rrs.wave_dependent_gordon(wavelengths)
    Rrs_var = bing_rt.calc_Rrs(a, bb, in_G1=G1, in_G2=G2)

    # Calculate differences
    abs_diff = Rrs_var - Rrs_std
    rel_diff = 100 * (Rrs_var - Rrs_std) / Rrs_std

    return Rrs_std, Rrs_var, abs_diff, rel_diff

# Example: oligotrophic water
wavelengths = np.arange(400, 701, 10)
a_oligo = np.array([0.03, 0.025, 0.02, 0.02, 0.025, 0.035, 0.05,
                    0.08, 0.12, 0.18, 0.25, 0.30, 0.35, 0.40,
                    0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75,
                    0.80, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10,
                    1.15, 1.20, 1.25])
bb_oligo = 0.001 * np.ones_like(a_oligo)

Rrs_std, Rrs_var, abs_diff, rel_diff = compare_gordon_formulations(
    a_oligo, bb_oligo, wavelengths)

# Differences typically largest at red/NIR wavelengths
print(f"Maximum relative difference: {np.max(np.abs(rel_diff)):.2f}%")
print(f"Wavelength of max difference: {wavelengths[np.argmax(np.abs(rel_diff))]:.0f} nm")

Typical relative differences:

Blue wavelengths (400-500 nm): 0.5-2%
Green wavelengths (500-600 nm): 1-5%
Red wavelengths (600-700 nm): 3-15%

Differences are largest in oligotrophic waters and at longer wavelengths.

Raman Scattering Correction

BING supports Raman scattering corrections to account for the inelastic scattering contribution to Rrs. This is particularly important for clear (oligotrophic) waters at wavelengths > 500 nm.

Integrated Raman Correction in calc_Rrs

The main calc_Rrs function accepts optional Raman parameters:

from bing.rt import rrs, raman
import numpy as np

# Define wavelengths
wave = np.arange(400, 701, 5)

# IOPs at emission (model) wavelengths
a = ...   # Total absorption [m^-1]
bb = ...  # Total backscattering [m^-1]

# Compute excitation wavelengths (~3400 cm^-1 Raman shift)
wave_ex = raman.emission_to_excitation_wavelength(wave)

# IOPs at excitation wavelengths
a_ex = ...   # Absorption at excitation wavelengths
bb_ex = ...  # Backscattering at excitation wavelengths

# Raman backscattering coefficient
bb_R = raman.raman_backscattering_coeff(wave_ex)

# Calculate Rrs with Raman correction
Rrs = rrs.calc_Rrs(a, bb, a_ex=a_ex, bb_ex=bb_ex, bb_R=bb_R)

The Raman correction factor is computed using Sathyendranath & Platt (1998):

# Calculate correction factor separately
corr = rrs.calc_raman_correction_factor(a, bb, a_ex, bb_ex, bb_R)
print(f"Raman correction range: {corr.min():.3f} - {corr.max():.3f}")

Raman Correction in Fitting Workflows

When fitting with BING, enable Raman correction in the parameter configuration:

from bing.parameters import standard

params = standard.expb_pow(
    satellite='PACE',
    include_Raman=True,      # Enable Raman correction
    variable_Gordon=True,    # Recommended for best accuracy
    nsteps=40000
)

The fitting module automatically:

Initializes excitation wavelengths for each model
Computes IOPs at both emission and excitation wavelengths
Applies the Raman correction factor during forward modeling

Note

For detailed Raman scattering physics and additional functions, see Raman Scattering.

Chlorophyll Fluorescence

BING also supports a chlorophyll fluorescence inelastic-emission contribution to Rrs. The top-level method bing.rt.calc_Rrs_fluorescence() integrates the fluorescence reflectance from a band of excitation wavelengths (370-690 nm) to emission wavelengths near the 685 nm primary peak (and optionally the 730 nm secondary peak), following Gordon (1979) and Sathyendranath & Platt (1998).

The fluorescence calculation depends on a downwelling irradiance spectrum at the excitation wavelengths. This is provided by the companion package correct_atmosphere:

Repository: https://github.com/ocean-colour/correct-atmosphere
Install: pip install git+https://github.com/ocean-colour/correct-atmosphere.git
Import: from correct_atmosphere import downwelling

import numpy as np
from bing.rt import calc_Rrs_fluorescence
from correct_atmosphere import downwelling

wave    = np.linspace(650, 750, 101)        # Emission wavelengths
wave_ex = np.arange(400.0, 681.0, 1.0)      # Excitation wavelengths

Ed_ex = downwelling.solar_irradiance(wave_ex)
Ed_em = downwelling.solar_irradiance(685.)

Rrs_fl = calc_Rrs_fluorescence(
    wave, a_em, bb_em,
    a_ex, bb_ex, aph_ex,
    wave_ex, Ed_ex, Ed_em,
    phi_C=0.02,
    double_gaussian=True,
)

For the full set of low-level fluorescence routines (emission line shapes, quantum yields, FLH, …) see Chlorophyll Fluorescence.

The Radiative-Transfer Dictionary (`rt_dict_from_p`)

The fitting routines (bing.fitting.chisq_fit, bing.fitting.inference, bing.fitting.l23) receive their radiative-transfer options through a small dictionary, the rt dict. The convenience method bing.rt.defs.rt_dict_from_p() builds this dictionary from a parameter named-tuple (see Parameters):

from bing.parameters import standard
from bing.rt import defs as rt_defs

# Build a parameter named-tuple (e.g. via standard.expb_pow(...))
p = standard.expb_pow(satellite='PACE')

# Or generate one explicitly with the RT options of interest
from bing.parameters import p_ntuple
p = p_ntuple.gen(
    model_names=['ExpBricaud', 'Pow'],
    variable_Gordon=True,
    include_Raman=True,
    include_Chl_fl=True,
    phi_C=0.02,
    double_gaussian=True,
)

# Convert to the radiative-transfer dictionary
rt_dict = rt_defs.rt_dict_from_p(p)
# rt_dict == {
#     'variable_Gordon':  True,
#     'include_Raman':    True,
#     'include_Chl_fl':   True,
#     'phi_C':            0.02,
#     'double_gaussian':  True,
# }

# Pass through to a fit
from bing.fitting import chisq_fit
ans, cov, idx = chisq_fit.fit(items[0], models, rt_dict, ...)

bing.rt.defs.rt_dict_from_p(p)

Build a radiative-transfer options dictionary from a BING parameter named-tuple.

Parameters:: p – Parameter named-tuple (typically produced by bing.parameters.p_ntuple.gen() or one of the standard.* helpers).
Returns:: dict with keys variable_Gordon, include_Raman, include_Chl_fl, phi_C, double_gaussian. Any attribute that is missing on p is set to None.
Return type:: dict

The returned dictionary is the canonical way to pass radiative-transfer options into the forward model and the fitting routines.

API Reference

Core Functions

calc_Rrs(a, bb, in_G1=None, in_G2=None, a_ex=None, bb_ex=None, bb_R=None)

Calculate remote sensing reflectance from IOPs using Gordon model, with optional Raman scattering correction.

Parameters:

a (float or ndarray) – Total absorption coefficient at emission wavelength [m^-1]
bb (float or ndarray) – Total backscattering coefficient at emission wavelength [m^-1]
in_G1 (float or ndarray, optional) – Gordon coefficient G1 (default: 0.0949)
in_G2 (float or ndarray, optional) – Gordon coefficient G2 (default: 0.0794)
a_ex (float or ndarray, optional) – Absorption at Raman excitation wavelength [m^-1]. Required for Raman correction.
bb_ex (float or ndarray, optional) – Backscattering at Raman excitation wavelength [m^-1]. Required for Raman correction.
bb_R (float or ndarray, optional) – Raman backscattering coefficient [m^-1]. Required for Raman correction.

Returns:

Remote sensing reflectance [sr^-1]

Return type:

float or ndarray

Raises:

IOError – If a_ex is provided but bb_ex or bb_R are not.

Examples:

# Elastic-only (standard Gordon)
Rrs = calc_Rrs(a=0.05, bb=0.002)

# Variable Gordon
wavelengths = np.arange(400, 701, 5)
G1, G2 = wave_dependent_gordon(wavelengths)
Rrs = calc_Rrs(a, bb, in_G1=G1, in_G2=G2)

# With Raman correction
from bing.rt import raman
wave_ex = raman.emission_to_excitation_wavelength(wavelengths)
bb_R = raman.raman_backscattering_coeff(wave_ex)
Rrs = calc_Rrs(a, bb, in_G1=G1, in_G2=G2,
               a_ex=a_ex, bb_ex=bb_ex, bb_R=bb_R)

calc_Rrs_fluorescence(wavelength, a_em, bb_em, a_ex, bb_ex, aph_ex, wavelength_ex, Ed_ex, Ed_em, mu_d=None, mu_f=None, phi_C=0.02, double_gaussian=True)

Calculate the Rrs contribution from chlorophyll fluorescence, integrating over excitation wavelengths.

Parameters:

wavelength (float or ndarray) – Emission wavelength(s) [nm]. Typically 650-750 nm.
a_em (float or ndarray) – Total absorption at emission wavelength(s) [m^-1]
bb_em (float or ndarray) – Total backscattering at emission wavelength(s) [m^-1]
a_ex (ndarray) – Total absorption at excitation wavelengths [m^-1]
bb_ex (ndarray) – Total backscattering at excitation wavelengths [m^-1]
aph_ex (ndarray) – Phytoplankton absorption at excitation wavelengths [m^-1]
wavelength_ex (ndarray) – Excitation wavelengths [nm]
Ed_ex (ndarray) – Downwelling irradiance at excitation wavelengths (from correct_atmosphere.downwelling)
Ed_em (float or ndarray) – Downwelling irradiance at the (peak) emission wavelength
phi_C (float, optional) – Quantum yield (default 0.02)
double_gaussian (bool, optional) – Use double-Gaussian emission line (default True)

Returns:

Fluorescence contribution to Rrs [sr^-1]

Return type:

ndarray

See Chlorophyll Fluorescence for full physics, low-level methods, and FLH calculations.

Note

The upwelling attenuation \(\kappa^F(\lambda) = (a(\lambda) + b_b(\lambda))/\mu_f\) is evaluated at every emission wavelength, not pinned to the 685 nm peak. This is what allows the double-Gaussian model to reproduce the correct (suppressed) 730 nm secondary shoulder, where pure-water absorption is ~4x larger than at 685 nm.

calc_raman_correction_factor(a_em, bb_em, a_ex, bb_ex, bb_R, Ed_ratio=1.0, include_second_order=True)

Calculate the multiplicative correction factor for Raman scattering.

Parameters:

a_em (float or ndarray) – Absorption at emission wavelength [m^-1]
bb_em (float or ndarray) – Backscattering at emission wavelength [m^-1]
a_ex (float or ndarray) – Absorption at excitation wavelength [m^-1]
bb_ex (float or ndarray) – Backscattering at excitation wavelength [m^-1]
bb_R (float or ndarray) – Raman backscattering coefficient [m^-1]
Ed_ratio (float, optional) – Ratio Ed(λ’)/Ed(λ). Default is 1.0
include_second_order (bool, optional) – Include second-order Raman terms. Default is True

Returns:

Correction factor (R_elastic + R_raman) / R_elastic

Return type:

float or ndarray

Typical correction factors range from 1.0 to ~1.25, with largest corrections in clear (oligotrophic) waters at longer wavelengths.

wave_dependent_gordon(wave, bounds_error=True)

Load wavelength-dependent Gordon coefficients.

Reads from bing/data/RT/gordon_coefficients.csv and interpolates to requested wavelengths using cubic splines.

Parameters:

wave (ndarray) – Wavelength array [nm]
bounds_error (bool, optional) – Raise error if wavelengths outside 350-750 nm range

Returns:

G1 and G2 arrays at requested wavelengths

Return type:

tuple(ndarray, ndarray)

Examples:

# Get coefficients for PACE wavelengths
pace_waves = np.arange(340, 701, 5)
G1, G2 = wave_dependent_gordon(pace_waves, bounds_error=False)

# Examine spectral structure
import matplotlib.pyplot as plt
plt.plot(pace_waves, G1, label='G1')
plt.plot(pace_waves, G2, label='G2')
plt.xlabel('Wavelength (nm)')
plt.ylabel('Gordon Coefficient')
plt.legend()

Conversion Functions

Rrs_to_rrs(Rrs, A=0.52, B=1.7)

Convert above-surface Rrs to subsurface rrs.

Parameters:

Rrs (ndarray) – Remote sensing reflectance (above surface) [sr^-1]
A (float, optional) – Conversion coefficient (default: 0.52)
B (float, optional) – Conversion coefficient (default: 1.7)

Returns:

Subsurface remote sensing reflectance [sr^-1]

Return type:

ndarray

\[r_{rs} = \frac{R_{rs}}{A + B \cdot R_{rs}}\]

rrs_to_Rrs(rrs, A=0.52, B=1.7)

Convert subsurface rrs to above-surface Rrs.

Parameters:

rrs (ndarray) – Subsurface remote sensing reflectance [sr^-1]
A (float, optional) – Conversion coefficient (default: 0.52)
B (float, optional) – Conversion coefficient (default: 1.7)

Returns:

Remote sensing reflectance (above surface) [sr^-1]

Return type:

ndarray

\[R_{rs} = \frac{A \cdot r_{rs}}{1 - B \cdot r_{rs}}\]

Constants

G1_STANDARD = 0.0949: Standard (constant) Gordon coefficient G1 from Gordon et al. (1988).

G2_STANDARD = 0.0794: Standard (constant) Gordon coefficient G2 from Gordon et al. (1988).

A_Rrs = 0.52: Conversion coefficient from subsurface to above-surface reflectance.

B_Rrs = 1.7: Conversion coefficient from subsurface to above-surface reflectance.

Complete Workflow Example

Forward Modeling with Variable Gordon

This example demonstrates a complete forward modeling workflow using the variable Gordon formulation:

import numpy as np
from bing import rt as bing_rt
from bing.models import utils as model_utils
from bing.parameters import standard
import matplotlib.pyplot as plt

# ===== 1. Setup =====
# Define hyperspectral wavelengths (PACE-like)
wavelengths = np.arange(400, 701, 5)

# Initialize bio-optical models
params = standard.expb_pow(satellite='PACE', variable_Gordon=True)
models = model_utils.init(params.model_names, wavelengths)

# ===== 2. Set bio-optical parameters =====
# Example: oligotrophic water
Chl = 0.1  # mg/m³
Sdg = 0.015  # CDOM spectral slope [nm^-1]
Adg_440 = 0.01  # CDOM absorption at 440 nm [m^-1]
Bnw_ref = 0.0008  # Particle backscatter at 600 nm [m^-1]
beta = 1.5  # Spectral slope for backscatter

# ===== 3. Configure models =====
# Set chlorophyll for phytoplankton absorption
models[0].set_aph(Chl=Chl)

# Model parameters (in log10 space for amplitudes)
model_params = np.array([
    np.log10(Adg_440),  # log10(A_dg)
    Sdg,                 # S_dg
    np.log10(0.05582 * Chl),  # log10(A_ph) from Bricaud
    np.log10(Bnw_ref),   # log10(B_nw)
    beta                 # spectral slope
])

# ===== 4. Evaluate IOPs =====
# Non-water absorption
anw = models[0].eval_anw(model_params[:3])

# Non-water backscattering
bbnw = models[1].eval_bbnw(model_params[3:])

# Add water components
a_total = models[0].a_w + anw[0]  # Total absorption
bb_total = models[1].bb_w + bbnw[0]  # Total backscattering

# ===== 5. Calculate Rrs with variable Gordon =====
G1, G2 = bing_rt.rrs.wave_dependent_gordon(wavelengths)
Rrs = bing_rt.calc_Rrs(a_total, bb_total, in_G1=G1, in_G2=G2)

# ===== 6. Compare with standard Gordon =====
Rrs_standard = bing_rt.calc_Rrs(a_total, bb_total)

# ===== 7. Visualize results =====
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot IOPs
axes[0, 0].semilogy(wavelengths, a_total, 'b-', label='Total absorption', linewidth=2)
axes[0, 0].semilogy(wavelengths, models[0].a_w, 'b--', label='Water absorption', alpha=0.5)
axes[0, 0].set_ylabel('Absorption [m$^{-1}$]')
axes[0, 0].set_title('Absorption Coefficient')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].semilogy(wavelengths, bb_total, 'r-', label='Total backscatter', linewidth=2)
axes[0, 1].semilogy(wavelengths, models[1].bb_w, 'r--', label='Water backscatter', alpha=0.5)
axes[0, 1].set_ylabel('Backscattering [m$^{-1}$]')
axes[0, 1].set_title('Backscattering Coefficient')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot Rrs comparison
axes[1, 0].plot(wavelengths, Rrs_standard, 'k--', label='Standard Gordon', linewidth=2)
axes[1, 0].plot(wavelengths, Rrs, 'g-', label='Variable Gordon', linewidth=2)
axes[1, 0].set_xlabel('Wavelength (nm)')
axes[1, 0].set_ylabel('R$_{rs}$ [sr$^{-1}$]')
axes[1, 0].set_title('Remote Sensing Reflectance')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Plot relative difference
rel_diff = 100 * (Rrs - Rrs_standard) / Rrs_standard
axes[1, 1].plot(wavelengths, rel_diff, 'purple', linewidth=2)
axes[1, 1].axhline(0, color='gray', linestyle=':', linewidth=1)
axes[1, 1].set_xlabel('Wavelength (nm)')
axes[1, 1].set_ylabel('Relative Difference (%)')
axes[1, 1].set_title('Variable - Standard Gordon')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('forward_model_variable_gordon.png', dpi=150)

# Print summary statistics
print(f"Chlorophyll: {Chl:.2f} mg/m³")
print(f"CDOM a_dg(440): {Adg_440:.4f} m⁻¹")
print(f"Mean Rrs difference: {np.mean(np.abs(rel_diff)):.2f}%")
print(f"Max Rrs difference: {np.max(np.abs(rel_diff)):.2f}% at {wavelengths[np.argmax(np.abs(rel_diff))]} nm")

Best Practices

Parameter Configuration

When using variable Gordon coefficients in fitting workflows:

from bing.parameters import standard

# Enable variable Gordon in parameter configuration
params = standard.expb_pow(
    satellite='PACE',
    variable_Gordon=True,  # Use wavelength-dependent coefficients
    add_noise=True,
    nsteps=40000
)

# The fitting module will automatically use G1(λ), G2(λ)

Error Handling

from bing.rt import rrs

# Handle wavelengths outside valid range
wavelengths = np.arange(300, 801, 5)  # Extends beyond 350-750 nm range

try:
    G1, G2 = rrs.wave_dependent_gordon(wavelengths, bounds_error=True)
except ValueError as e:
    print(f"Wavelength out of range: {e}")
    # Use bounds_error=False to extrapolate
    G1, G2 = rrs.wave_dependent_gordon(wavelengths, bounds_error=False)
    print("Warning: Extrapolating Gordon coefficients beyond valid range")

Physical Constraints

Gordon coefficients must satisfy physical constraints:

G1 > 0 (always positive)
|G2| < 1 (bounded, can be negative)
G1 >> |G2| at short wavelengths
Rrs >= 0 for physical solutions

# Verify physical constraints
def check_gordon_physics(G1, G2, wavelengths):
    """Validate Gordon coefficients."""

    checks = {
        'G1_positive': np.all(G1 > 0),
        'G2_bounded': np.all(np.abs(G2) < 1.0),
        'G1_magnitude': np.all(G1 > np.abs(G2))
    }

    for check, passed in checks.items():
        status = "✓" if passed else "✗"
        print(f"{status} {check}: {passed}")

    return all(checks.values())

# Run validation
wavelengths = np.arange(400, 701, 5)
G1, G2 = rrs.wave_dependent_gordon(wavelengths)
is_valid = check_gordon_physics(G1, G2, wavelengths)

Performance Considerations

Computation Speed

Loading variable Gordon coefficients requires:

Reading CSV file (~100 rows)
Cubic spline interpolation to model wavelengths

This adds negligible overhead (~1 ms) and should be done once per model initialization.

import time

# Benchmark
wavelengths = np.arange(400, 701, 1)  # 301 wavelengths

# Standard Gordon (constant)
t0 = time.time()
for _ in range(1000):
    Rrs = bing_rt.calc_Rrs(a, bb)
t_standard = time.time() - t0

# Variable Gordon (precomputed)
G1, G2 = rrs.wave_dependent_gordon(wavelengths)
t0 = time.time()
for _ in range(1000):
    Rrs = bing_rt.calc_Rrs(a, bb, in_G1=G1, in_G2=G2)
t_variable = time.time() - t0

print(f"Standard Gordon: {t_standard:.3f} s")
print(f"Variable Gordon: {t_variable:.3f} s")
print(f"Overhead: {100*(t_variable-t_standard)/t_standard:.1f}%")

Typical overhead: <1% when G1, G2 are precomputed.

Memory Usage

Variable Gordon coefficients add:

2 arrays × n_wavelengths × 8 bytes
Example: 301 wavelengths = ~5 KB

Negligible compared to MCMC chains or model data.

References

Gordon Semi-Analytical Model:

Gordon, H.R. et al. (1988). “A semianalytic radiance model of ocean color,” Journal of Geophysical Research 93(D9), 10909-10924. https://doi.org/10.1029/JD093iD09p10909

Variable Gordon Coefficients:

Loisel, H., Vantrepotte, V., et al. (2023). “Ocean color remote sensing of coastal waters: Are we ready for hyperspectral satellites?” Journal of Geophysical Research: Oceans (in preparation).

Radiative Transfer Theory:

Mobley, C.D. (1994). Light and Water: Radiative Transfer in Natural Waters. Academic Press.
Morel, A. and Gentili, B. (1996). “Diffuse reflectance of oceanic waters. III. Implication of bidirectionality for the remote-sensing problem,” Applied Optics 35(24), 4850-4862.

Raman Scattering:

Sathyendranath, S. and Platt, T. (1998). “Ocean-colour model incorporating transspectral processes,” Applied Optics 37, 2216-2227.
Bartlett, J.S. et al. (1998). “Raman scattering by pure water and seawater,” Applied Optics 37, 3324-3332.

Related BING Documentation:

Bio-optical Models - Bio-optical model formulations
Fitting Algorithms - Parameter estimation algorithms
Data Processing - End-to-end workflow examples

Radiative Transfer

Overview

Gordon Semi-Analytical Model

Standard Gordon Coefficients

Variable Gordon Coefficients

Loading Variable Gordon Coefficients

Gordon Coefficient Data Source

When to Use Variable Gordon Coefficients

Comparison of Standard vs. Variable Gordon

Spectral Differences

Impact on Retrieved Rrs

Raman Scattering Correction

Integrated Raman Correction in calc_Rrs

Raman Correction in Fitting Workflows

Chlorophyll Fluorescence

The Radiative-Transfer Dictionary (rt_dict_from_p)

API Reference

Core Functions

Conversion Functions

Constants

Complete Workflow Example

Forward Modeling with Variable Gordon

Best Practices

Parameter Configuration

Error Handling

Physical Constraints

Performance Considerations

Computation Speed

Memory Usage

References

See Also

The Radiative-Transfer Dictionary (`rt_dict_from_p`)