vSmartMOM

Appearance

Surface BRDF Models: Theory and Usage ​

vSmartMOM supports a family of Bidirectional Reflectance Distribution Function (BRDF) models for the lower boundary condition in radiative transfer. This tutorial covers the mathematical background, code usage, and guidance for every available surface type – from the simplest isotropic Lambertian to the fully polarized Cox-Munk ocean surface.

Contents

  1. Overview and model comparison

  2. Lambertian surfaces

  3. RPV (Rahman-Pinty-Verstraete) model

  4. Ross-Li kernel-based model

  5. Cox-Munk polarized ocean surface

  6. Canopy-coupled surfaces

  7. Polarization and Mueller matrices

  8. Surface Jacobians (linearized RT)

  9. Choosing the right model

  10. Running a full example

  11. YAML configuration

  12. References

julia
using vSmartMOM
using vSmartMOM.CoreRT
using vSmartMOM.CoreRT: rpvSurfaceScalar, RossLiSurfaceScalar, LambertianSurfaceLegendre, LambertianSurfaceSpline
using CairoMakie

1) Overview ​

Every surface model is a subtype of AbstractSurfaceType. The RT solver dispatches through create_surface_layer!, which computes the Fourier-decomposed surface reflectance matrix for each azimuthal moment .

ModelConstructorParamsPolarized?Jacobians?Use case
Lambertian (scalar)LambertianSurfaceScalar(α)1noyesSimple land, snow, dark ocean
Lambertian (spectral)LambertianSurfaceSpectrum([α₁,…])NnonoWavelength-varying surfaces
RPVrpvSurfaceScalar(ρ₀, ρ_c, k, Θ)4nonoVegetation, soils
Ross-LiRossLiSurfaceScalar(fvol, fgeo, fiso)3nonoMODIS-style vegetation
Cox-MunkCoxMunkSurface(wind_speed=U,…)1-5yesyesOcean, sun glint
CanopyCanopySurface(soil=…, LAI=…,…)manyvia soilpartialVegetation canopy + soil

Two additional Lambertian variants (LambertianSurfaceLegendre, LambertianSurfaceSpline) are available internally for polynomial or spline-based spectral albedo but are not exported.

2) Lambertian Surface ​

The simplest model: reflectance is independent of viewing and illumination angles. A perfect Lambertian surface satisfies

where   is the hemispherical albedo. Because the BRDF has no azimuthal dependence, only the   Fourier moment is nonzero. The code applies a factor of 2 to   for consistency with the Fourier series normalization used in the RT solver.

Scalar albedo ​

Use LambertianSurfaceScalar for a single albedo across all wavelengths in the band:

julia
lam_scalar = LambertianSurfaceScalar(0.15)

Spectral albedo ​

Use LambertianSurfaceSpectrum when the albedo varies per spectral point (one value per grid point in the band):

julia
lam_spectrum = LambertianSurfaceSpectrum([0.12, 0.13, 0.14, 0.15])

Legendre and spline variants ​

Not exported, but accessible from CoreRT:

  • CoreRT.LambertianSurfaceLegendre([c₀, c₁, …]): polynomial spectral albedo    , where is the spectral grid rescaled to  .

  • CoreRT.LambertianSurfaceSpline(interpolator, wlGrid): spline-interpolated albedo from a user-supplied Interpolations.jl object.

3) RPV (Rahman-Pinty-Verstraete) Model ​

The RPV model (Rahman, Pinty & Verstraete, 1993) provides a four-parameter empirical BRDF widely used for vegetation and soils:

where the three factors are:

Minnaert function (angular shape):    

Hot-spot function (Henyey-Greenstein-like):    

Geometric bowl function:     

with the scattering angle cosine     and the geometric factor    .

ParameterSymbolTypical rangePhysical meaning
Overall amplitude0 - 0.5Isotropic reflectance level
Minnaert exponent0.3 - 1.5 : bowl shape,  : bell shape
Hot-spot asymmetry to  : backscatter,  : forward
Geometric amplitude to 1 = no geometric term

RPV is scalar-only: it returns zero for Stokes components beyond .

julia
rpv = rpvSurfaceScalar(0.15, 0.1, 0.7, -0.3)  ## ρ₀, ρ_c, k, Θ

4) Ross-Li Model ​

The Ross-Li model (Lucht, Schaaf & Strahler, 2000) decomposes the BRDF into a linear combination of semi-physical kernels:

Isotropic kernel:  .

Ross Thick volumetric kernel (dense canopy scattering):    where    is the scattering angle.

Li Sparse geometric kernel (shadowing by sparsely distributed objects):       where is the overlap function and primed angles are scaled by the crown shape ratio ( ,   RAMI defaults).

ParameterSymbolTypical rangePhysical meaning
Isotropic0 - 0.3Baseline reflectance
Volumetric0 - 0.1Dense-canopy scattering
Geometric0 - 0.1Mutual shadowing

Ross-Li is scalar-only: it returns zero for Stokes components beyond .

By default both kernel zenith angles are capped at 80° (zenith_cap_deg = 80), a vSmartMOM regularization of the grazing quadrature nodes, where the Li-sparse kernel turns strongly negative. zenith_cap_deg = nothing uses the true angles, i.e. the MODIS kernel definition that VLIDORT implements.

julia
rossli = RossLiSurfaceScalar(0.05, 0.03, 0.1)  ## fvol, fgeo, fiso
rossli_modis = RossLiSurfaceScalar(0.05, 0.03, 0.1; zenith_cap_deg = nothing)

5) Cox-Munk Polarized Ocean Surface ​

CoxMunkSurface is the first fully polarized surface model in vSmartMOM. It models the ocean as an ensemble of specular wave facets whose tilt distribution follows a wind-dependent Gaussian. Fresnel reflection from each facet produces a full Mueller-matrix BRDF that couples all four Stokes components. This is critical for remote sensing of ocean scenes (OCO-2/3, PACE, etc.) where sun glint polarization is a dominant signal.

5a) Slope distribution (Cox & Munk, 1954) ​

The ocean surface is approximated as a collection of small flat facets ("capillary waves") whose slopes follow an isotropic Gaussian:

where is the variance of each slope component. Cox & Munk (1954) measured the total mean-square slope   , which depends linearly on the 10-m wind speed (m/s):

Low wind (  m/s) gives a nearly mirror-like surface with a narrow specular glint; high wind (  m/s) spreads the glint over a wide angular range.

5b) Fresnel reflection ​

Each facet reflects specularly according to Fresnel's equations. For a complex relative refractive index     (water/air) and local incidence angle :

where   (Snell's law).

The Fresnel Mueller matrix for reflection is:

It refers to the facet basis   , with normal to the plane of reflection. The off-diagonal (1,2) and (2,1) elements couple the and Stokes parameters – this is the polarizing effect of the ocean surface.

5c) BRDF Mueller matrix ​

For incident cosine (azimuth 0), reflected cosine and relative azimuth (the glint lies at  ) the full BRDF is:

where: ​

$$\beta$$
is the facet tilt angle ($\cos\beta = n_z$ for the facet normal $\hat n \propto \mathbf{k}_r - \mathbf{k}_i$)

  • is the local incidence angle on the tilted facet

  • changes the Stokes basis between a beam's meridian plane and the facet basis . The field components transform as   with      — a reflection of the basis, which is its own inverse — so

belongs to the incident and to the reflected beam. ​

$$S$$
is the Smith (1967) bistatic shadow masking factor

5d) Shadow masking (Smith, 1967) ​

At low grazing angles, parts of the surface are shadowed by neighboring wave crests. The Smith shadowing function corrects for this:

where   with    and  .

5e) Whitecap contribution (Monahan & O'Muircheartaigh, 1980) ​

At higher wind speeds, breaking waves form whitecaps that scatter light diffusely. The fractional whitecap coverage is:

Whitecaps are modeled as unpolarized Lambertian reflectors with albedo (default 0.22, Koepke 1984). The total BRDF is:

where has only a element equal to .

5f) Water refractive index ​

When n_water = nothing (the default), the solver uses the Segelstein (1981) table (34 nm – 10 m) at the wavelength of every spectral point. The function water_refractive_index(λ_nm) returns the complex refractive index at any wavelength, with n and k interpolated linearly in wavelength:

julia
n_water_550 = water_refractive_index(550.0)
println("Water RI at 550 nm:  n = ", round(real(n_water_550), digits=4),
        ", k = ", round(imag(n_water_550), sigdigits=3))

You can also provide a fixed value for the whole band, or a vector with one value per spectral point (used element-wise), via the n_water keyword of CoxMunkSurface.

5g) TMS single-scattering correction ​

The specular sun-glint peak can be extremely narrow (especially at low wind speed), requiring many Fourier moments to resolve. Rather than increasing max_m, the code applies a Truncated Multiple Scattering (TMS) correction after the Fourier loop: it adds the difference between the exact single-scattering surface contribution and the truncated Fourier reconstruction at each viewing geometry. This is handled automatically by apply_ss_correction! inside rt_run, for the TOA radiance only: upwelling radiances at interior observer heights (obs_alt) keep the Fourier-series glint, in the forward and the linearized solvers alike.

5h) Construction ​

julia
ocean = CoxMunkSurface(wind_speed = 5.0)  ## uses built-in water RI, whitecaps on

ocean_custom = CoxMunkSurface(
    wind_speed      = 7.5,
    n_water         = complex(1.34, 1e-7),  ## override refractive index
    whitecap_albedo = 0.22,
    include_whitecaps = true,
    shadowing       = true,
)

6) Canopy-Coupled Surface ​

CanopySurface represents a vegetation canopy backed by a soil BRDF. It internally solves canopy sub-layers via the adding-doubling method before presenting an effective reflectance to the atmospheric RT. This is covered in detail in the Canopy Tutorial; here is a brief example:

julia
canopy = CanopySurface(
    soil = LambertianSurfaceScalar(0.1),
    LAI  = 3.0,
    n_layers = 1,
    leaf_reflectance   = 0.45,
    leaf_transmittance = 0.05,
)

CanopySurface_from_prospect(leaf, wl_grid; ...) provides a convenience constructor that computes spectral leaf optics from the PROSPECT model. See the Canopy Tutorial for multi-layer setups, spectral leaf optics, and within-canopy atmospheric absorption.

7) Polarization and Mueller Matrices ​

The Lambertian, RPV, and Ross-Li models are scalar-only: they populate only the Stokes block of the reflectance matrix and return zero for all off-diagonal blocks. This is sufficient for intensity-only () RT.

The Cox-Munk surface is the first model to fill the full Mueller matrix at the surface boundary. The Fresnel reflection off tilted wave facets couples the - and - Stokes components. Concretely:

  • For Stokes_I(): the reflectance matrix is  .

  • For Stokes_IQUV(): it becomes   with all Stokes cross-coupling blocks populated.

The Fourier decomposition uses polarization-aware azimuthal kernels: for the - / - blocks and for the cross blocks, matching the postprocessing reconstruction convention.

8) Surface Jacobians (Linearized RT) ​

Several surface types support analytical Jacobians for use with rt_run_lin:

Lambertian (analytical) ​

The derivative of the surface reflectance matrix w.r.t. albedo is trivially:       for  , where is the identity-like Stokes selector and are the quadrature weights.

Cox-Munk (analytical) ​

The derivative of the surface reflectance w.r.t. wind speed is computed analytically via the chain rule:         

where  . The geometry, Fresnel coefficients, and Stokes basis changes are all independent of and shared with the forward evaluation.

Canopy ​

Jacobians with respect to the soil albedo use finite-difference perturbation of the internal adding-doubling solve.

RPV, Ross-Li ​

No built-in linearization. Use external finite differences or ForwardDiff.

9) Choosing the Right Model ​

ScenarioRecommendedReason
Quick test, isotropic surfaceLambertianSurfaceScalarFastest, simplest
Ocean, sun glint, polarimetryCoxMunkSurfaceFull polarization, wind-speed Jacobians
Vegetation (MODIS/RAMI)RossLiSurfaceScalarStandard kernel-based, 3 params
Vegetation (empirical)rpvSurfaceScalar4-param empirical with hotspot
Physical canopy modelCanopySurfaceLeaf-level scattering, spectral optics
Spectral albedo from dataLambertianSurfaceSpectrumPer-wavelength albedo vector
Polarized atmosphere, scalar surfaceLambertianSurfaceScalarOnly block filled

10) Running a Full Example ​

Load a parameter file, override the surface, and run the RT:

julia
yaml_path = joinpath(pkgdir(vSmartMOM),
                     "test", "test_parameters", "PureRayleighParameters.yaml")
params = read_parameters(yaml_path)
params.brdf[1] = CoxMunkSurface(wind_speed = 5.0)
model = model_from_parameters(params)
R, T = rt_run(model)
println("R shape: ", size(R))
println("R(VZA 60°, I): ", R[1, 1, 1])

Compare with a Lambertian surface at the same geometry:

julia
params2 = read_parameters(yaml_path)
params2.brdf[1] = LambertianSurfaceScalar(0.06)
model2 = model_from_parameters(params2)
R2, T2 = rt_run(model2)
println("R(VZA 60°, I) Lambertian 0.06: ", R2[1, 1, 1])

Compare the radiance spectra from the two surface types:

julia
fig = Figure(size=(700, 450))
ax = Axis(fig[1,1],
    xlabel = "Spectral index",
    ylabel = "TOA radiance (Stokes I)")
lines!(ax, R[1, 1, :],  label="Cox-Munk (U=5 m/s)")
lines!(ax, R2[1, 1, :], label="Lambertian (α=0.06)")
axislegend(ax, position=:rt)
fig

The rendered docs include a Plotly comparison of the angular signatures for the main surface families:

The interactive panel uses synthetic angular shapes to illustrate these surface families. It is not output from the BRDF implementations and is not a validation reference.

11) YAML Configuration ​

In a YAML parameter file, the surface is specified as a constructor string under radiative_transfer.surface (one entry per spectral band):

yaml
radiative_transfer:
  surface:
    - LambertianSurfaceScalar(0.15)
    - rpvSurfaceScalar(0.12, 0.08, 0.75, -0.25)
    - RossLiSurfaceScalar(0.04, 0.02, 0.08)
    - CoxMunkSurface(wind_speed=5.0)

For CoxMunkSurface, optional keyword arguments can be specified inline:

yaml
    - CoxMunkSurface(wind_speed=7.5, include_whitecaps=true, shadowing=true)

For canopy surfaces, use the dedicated canopy YAML section (see the Canopy Tutorial), which wraps the surface entry as the soil BRDF.

12) References ​

  • Cox, C. & Munk, W. (1954). Measurement of the roughness of the sea surface from photographs of the sun's glitter. JOSA 44(11), 838-850. doi:10.1364/JOSA.44.000838

  • Koepke, P. (1984). Effective reflectance of oceanic whitecaps. Applied Optics 23(11), 1816-1824. doi:10.1364/AO.23.001816

  • Lucht, W., Schaaf, C.B. & Strahler, A.H. (2000). An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Trans. Geosci. Remote Sens. 38(2), 977-998. doi:10.1109/36.841980

  • Mishchenko, M.I. & Travis, L.D. (1997). Satellite retrieval of aerosol properties over the ocean using polarization as well as intensity of reflected sunlight. JGR 102(D14), 16989-17013. doi:10.1029/97JD01084

  • Monahan, E.C. & O'Muircheartaigh, I.G. (1980). Optimal power-law description of oceanic whitecap coverage dependence on wind speed. J. Phys. Oceanogr. 10, 2094-2099.

  • Rahman, H., Pinty, B. & Verstraete, M.M. (1993). Coupled surface-atmosphere reflectance (CSAR) model: 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data. JGR 98(D11), 20791-20801. doi:10.1029/93JD00564

  • Segelstein, D.J. (1981). The complex refractive index of water. M.S. Thesis, University of Missouri-Kansas City.

  • Smith, B.G. (1967). Geometrical shadowing of a random rough surface. IEEE Trans. Antennas Propag. 15(5), 668-671. doi:10.1109/TAP.1967.1138991

  • Zhai, P.-W., Hu, Y., Chowdhary, J., Trepte, C.R., Lucker, P.L. & Josset, D.B. (2010). A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface. JQSRT 111(7-8), 1025-1040. doi:10.1016/j.jqsrt.2009.12.005

This page was generated using Literate.jl.