light_model.tex

\chapter{LIGHT MODEL}
\label{chap:light}

Now that the spatial kelp distribution has been modeled, the radiative transfer equation is introduced, which is used to calculate the light field.
An asymptotic series centered at the case of no scattering is developed, which forms the basis for the faster and less memory-intensive of the two solution algorithms presented in Chapter \Rom{\ref{chap:numerical}}.

\section{Optical Definitions}
Before introducing the radiative transfer equation, it is necessary to discuss some basic radiometric quantities of interest which characterize the light field, as well as inherent optical properties which describe the medium of propagation.
It is necessary to begin by saying that the study of light is riddled with different sets of units which vary between disciplines.

From a physics point of view, light is a form of electromagnetic radiation which carries energy determined exclusively by its wavelength.
\textit{Radiometric units} thus describe the transfer of energy, measured in Watts, with quantities such as radiance, irradiance, and radiant flux to describe various types of energy transfer.
These quantities come in frequency-dependent and frequency-integrated varieties, depending on the context.
From the standpoint of human perception, the importance of light is not in the raw energy that it transfers, but the degree to which it facilitates vision from a physiological perspective.

Therefore, \textit{Photometric units} take all of the quantities from radiometry, rename them, and weight them by a \textit{luminosity function} which describes frequency dependence of the human eye's sensitivity to light.
Various luminosity functions exist which describe the eye's response to light in different circumstances.
In photometric units, radiance becomes luminance, irradiance illuminance, and so on.

Meanwhile, plant biologists prefer to measure light by counting photons in the photosynthetic frequency band.
The most common quantity in plant biology is \textit{photosynthetically active radiation} (PAR), which is generally measured in moles of photons per square meter per second.
All photons in the photosynthetic band (\SI{400}{\nm}--\SI{700}{\nm}) may be counted equally or weighted according to a plant's photosynthetic response.
Radiometric quantities are used throughout this thesis.

\subsection{Radiometric Quantities}
One of the most fundamental quantities in optics is radiant flux $\Phi$, which has units of energy per time (Watts).
Considering an element of surface area $A$, the energy \textit{density} moving through the surface is called irradiance, denoted $I$, and has units energy per time.
Note that the angle $\vartheta$ of the surface relative to the light source should also be considered for full detail. Assuming the surface to be flat and the light source to be distant (parallel rays), the flux through the surface is $\Phi=IA\cos\vartheta$.
Further, the angular dependence of the light field must be considered.
The radiance, $L$, expresses this dependence, and is defined as the radiant flux per steradian per projected surface area perpendicular to the direction of propagation of the beam \cite{chandrasekhar_radiative_1960}.
That is,
\begin{equation*}
	L = \frac{d^2\Phi}{dA d\vec{\omega}},
\end{equation*}
where $\vec{\omega}$ is an element of solid angle.
Once the radiance $L$ is calculated everywhere, the irradiance is
\begin{equation*}
  I(\vec{x}) = \int_{4\pi}L(\vec{x},\vec{\omega})\, d\vec{\omega},
\end{equation*}
This quantity is also called \textit{scalar irradiance} since it does not consider light through particular surface, but rather weights radiance from all directions equally \cite{kirk_light_1994}.
For brevity, this quantity is simply called irradiance here.

Irradiance can be approximately converted to PAR units of moles of photons (a mole of photons is also called an Einstein) per second, with the conversion being
\begin{equation}
  \SI{1}{\W} = \SI{4.2}{\micro\mole \,photons\per\meter\squared\per\second}.
  \label{eqn:watts_photons}
\end{equation}
This is not an exact conversion, but has been found to be accurate to roughly $\pm10\%$ across a variety of waters \cite{mobley_light_1994}.

\subsection{Perceived Irradiance}
\label{sec:perceived_irrad}
Assuming that the irradiance $I(\vec{x})$ is known,
the average irradiance at a depth $z$ is calculated as
\begin{equation}
  %\bar{I}_k = \frac{\sum_{ij}I_{ijk}}{n_x n_y},
  \bar{I}(z) = \frac{\iint I(x,y,z)\, dx\, dy}{\iint 1\, dx\, dy}.
\end{equation}
More relevant, however, is the average irradiance perceived by the kelp.
To calculate this value, the irradiance is weighted by the
normalized spatial kelp distribution before taking the mean.
Then, the average perceived irradiance at each depth is
\newcommand{\Iperk}{I_{\rm{perc}}}
\begin{equation}
   %\Iperk = \frac{\sum_{ij}P_{ijk}I_{ijk}}{\sum_{ij}P_{ijk}}.
   \Iperk(z) = \frac{\iint P_k(x,y,z)I(x,y,z)\, dx\, dy}{\iint P_k(x,y,z)\, dx\, dy}.
\end{equation}
The irradiance perceived by the kelp is expected to be lower than the average irradiance,
since the kelp is more densely located at the center of the domain where the light field is reduced,
whereas the simple average is influenced by regions of higher irradiance at the edges of the domain where kelp is not present.

\subsection{Inherent Optical Properties}
\label{sec:iops}
We now define a few inherent optical properties (IOPs) which depend only on the medium of propagation.
The absorption coefficient $a(\vec{x})$ (units m$^{-1}$) defines the
proportional loss of radiance per unit length due to absorption by the medium.
For example, this includes radiant energy which is converted to heat.
The scattering coefficient $b$ (units m$^{-1}$), defines the proportional loss
of radiance per unit length due to scattering, and is assumed to be constant over space.
Scattered light is not lost from the light field, it simply changes direction.

The volume scattering function (VSF) $\beta(\Delta): [-1, 1] \to \RR^+$ (units sr$^{-1}$) defines the probability of light scattering at any given angle from its source, where $\Delta=\cos\vartheta$ is the cosine of the angle between the initial and final directions.
Formally, given two directions $\vec{\omega}$ and $\vec{\omega}'$, $\beta(\vec{\omega} \cdot \vec{\omega}')$ is the probability density of light scattering from $\vec{\omega}$ into $\vec{\omega}'$ (or vice-versa).
Now, since a single direction subtends no solid angle, the probability of scattering occurring exactly from $\vec{\omega}$ to $\vec{\omega}'$ is 0.
Rather, we say that the probability of radiance being scattered from a direction $\omega$ into an element of solid angle $\Omega$ is $\int_\Omega \beta(\vec{\omega} \cdot \vec{\omega}')\, d\vec{\omega}'$.

The VSF is normalized such that
\begin{equation*}
  \int_{-1}^1\beta(\Delta)\, d\Delta=\frac{1}{2\pi},
\end{equation*}
so that for any $\omega$,
\begin{equation*}
  \int_{4\pi}\beta(\vec{\omega}\cdot\vec{\omega}')\, d\vec{\omega}' = 1,
\end{equation*}
i.e., the probability of light being scattered to some direction on the unit sphere is 1.

\section{The Radiative Transfer Equation}
We are now prepared to present the full details of radiative transfer equation, whose solution is the radiance in the medium as a function of position $\vec{x}$ and angle $\vec{\omega}$.

\subsection{Ray Notation}
\label{sec:ray_notation}
Consider a fixed position $\vec{x}$ and direction $\vec{\omega}$ such that
$\vec{\omega} \cdot \hat{z} \neq 0$ (the ray is not horizontal).
Let $\vec{l}(\vec{x}, \vec{\omega}, s)$ denote the linear path from one domain boundary to another containing $\vec{x}$ in the direction $\vec{\omega}$.
Since the ray is not horizontal, it originates either at the surface or bottom of the domain, with initial $z$ coordinate given by
\begin{equation*}
  z_0 =
   \begin{cases}
    0, & \vec{\omega} \cdot \hat{z} > 0 \\
    \zmax, & \vec{\omega} \cdot \hat{z} < 0.
  \end{cases}
\end{equation*}
Hence, the ray path is parameterized as
\begin{equation}
  \vec{l}(\vec{x}, \vec{\omega}, s) = \frac{1}{\tilde{s}} (s\vec{x} + (\tilde{s} - s)\vec{x_0}(\vec{x}, \vec{\omega})),
  \label{eqn:ray_path}
\end{equation}
where
\begin{equation}
  \vec{x_0}(\vec{x}, \vec{\omega}) = \vec{x} - \tilde{s} \vec{\omega}
  \label{eqn:x_0}
\end{equation}
is the origin of the ray, and
\begin{equation*}
  \tilde{s} = \frac{\vec{x} \cdot \hat{z} - z_0}{\vec{\omega} \cdot \hat{z}}
\end{equation*}
is the path length from $\vec{x_0}(\vec{x}, \vec{\omega})$ to $\vec{x}$.

\subsection{Colloquial Description}
Denote the radiance at $\vec{x}$ in the direction $\vec{\omega}$ by $L(\vec{x}, \vec{\omega})$.
As light travels along $\vec{l}(\vec{x}, \vec{\omega}, s)$, interaction with the
medium produces four phenomena of interest:
\begin{enumerate}
  \item Radiance is decreased due to absorption.
  \item Radiance is decreased due to scattering out of the path to other
    directions.
  \item Radiance is increased due to scattering into the path from other
      directions.
  \item Radiance is increased or decreased due to light sources or sinks.
\end{enumerate}

\subsection{Equation of Transfer}
Combining these phenomena yields the radiative transfer equation along
$\vec{l}(\vec{x}, \vec{\omega}, s)$ evaluated at $(\vec{x}, \vec{\omega})$,
\begin{align}
  \label{eqn:rte1d}
  \frac{dL\left(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}\right)}{ds}\Bigg|_{s=\tilde{s}}
  = &-(a(\vec{x}) + b)L(\vec{x}, \vec{\omega}) \nonumber \\
  &+ b \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L(\vec{x})\, d\vec{\omega}' + \sigma(\vec{x}, \vec{\omega}),
\end{align}
where $\int_{4\pi}$ denotes integration over the unit sphere.
The derivative of $L$ over the path can be rewritten as
\begin{align*}
  \frac{dL\left(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}\right)}{ds}\Bigg|_{s=\tilde{s}}
    &= \left[\frac{d\vec{l}(\vec{x}, \vec{\omega}, s)}{ds} \cdot \nabla L(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}', \vec{\omega})\right]\Bigg|_{s=\tilde{s}} \\
    &= \vec{\omega} \cdot \nabla L(\vec{x}, \vec{\omega}),
\end{align*}
which reveals the vector form of the radiative transfer equation,
\begin{equation*}
  \vec{\omega} \cdot \nabla L(\vec{x}, \vec{\omega})
  = -(a(\vec{x}) + b)L(\vec{x}, \vec{\omega})
  + b \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L(\vec{x}, \vec{\omega}')\, d\vec{\omega}' + \sigma(\vec{x}, \vec{\omega}),
\end{equation*}
or equivalently,
\begin{equation}
  \vec{\omega} \cdot \nabla L(\vec{x}, \vec{\omega})
  + a(\vec{x})L(\vec{x}, \vec{\omega})
  = b \left(
    \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L(\vec{x}, \vec{\omega}')\, d\vec{\omega}'
    - L(\vec{x}, \vec{\omega})
  \right)+ \sigma(\vec{x}, \vec{\omega}).
  \label{eqn:rte}
\end{equation}

\subsection{Boundary Conditions}

We use periodic boundary conditions in the $x$ and $y$ directions,
\begin{align*}
  L\left((\xmin, y, z), \vec{\omega}\right) &= L\left((\xmax, y, z), \vec{\omega}\right), \\
  L\left((x, \ymin, z), \vec{\omega}\right) &= L\left((x, \ymax, z), \vec{\omega}\right).
\end{align*}
In the $z$ direction, we specify a spatially uniform downwelling light just
under the surface of the water by a function $f(\vec{\omega})$.
Or if $\zmin>0$, then the radiance at $z=\zmin$ should be specified instead (as opposed to the radiance at the first grid cell center).
Further, we assume that no upwelling light enters the domain from the bottom.
Letting $\vec{x_s}$ be a point on the surface of the domain and $\vec{x_b}$ a point on the bottom, we have
\begin{align*}
  L(\vec{x_s}, \vec{\omega}) &= f(\vec{\omega}) \mbox{ if } \vec{\omega} \cdot \hat{z} > 0, \\
  L(\vec{x_b}, \vec{\omega}) &= 0 \mbox { if } \vec{\omega} \cdot \hat{z} < 0.
\end{align*}

\section{Low-Scattering Approximation}
\label{sec:symbolic_asymptotics}
In waters where absorption dominates scattering, an asymptotic series in terms of the scattering coefficient $b$ can be constructed.
The physical interpretation of the asymptotic series is that each term represents a discrete scattering event.
With the addition of each term, light from the previous term is scattered and attenuated from each point along the ray path.
In reality, the scattering cannot be considered to occur in discrete events, but rather all scattering occurs simultaneously (on a macroscopic timescale).

Since this is only an approximation, it is important to note that while the asymptotic series converges as $b \to 0$, it is not necessarily true that the series converges as the number of terms increases, although it may occur in certain cases.
Especially in cases of large scattering, the asymptotic series diverges rapidly.
The convergence properties of the algorithm are discussed in detail in Chapter \Rom{\ref{chap:model_analysis}}.

\subsection{Asymptotic Expansion}
Taking $b$ to be small, we introduce the asymptotic series
\newcommand{\Lasym}{L_0(\vec{x},\vec{\omega}) + b L_1(\vec{x},\vec{\omega}) + b^2 L_2(\vec{x},\vec{\omega}) + \cdots}
\newcommand{\Lasyms}{L_0(\vec{x_s},\vec{\omega}) + b L_1(\vec{x_s},\vec{\omega}) + b^2 L_2(\vec{x_s},\vec{\omega}) + \cdots}
\newcommand{\Lasymb}{L_0(\vec{x_b},\vec{\omega}) + b L_1(\vec{x_b},\vec{\omega}) + b^2 L_2(\vec{x_b},\vec{\omega}) + \cdots}
\newcommand{\Lasymp}{L_0(\vec{x},\vec{\omega}') + b L_1(\vec{x},\vec{\omega}') + b^2 L_2(\vec{x},\vec{\omega}') + \cdots}
\newcommand{\sigasym}{\sigma_0(\vec{x},\vec{\omega}) + b \sigma_1(\vec{x},\vec{\omega}) + b^2 \sigma_2(\vec{x},\vec{\omega}) + \cdots}
\begin{equation*}
  \label{eqn:Lasym}
  L(\vec{x},\vec{\omega}) = \Lasym.
\end{equation*}
Since the source $\sigma$ may also depend on $b$, it deserves a similar expansion,
\begin{equation*}
  \label{eqn:sigasym}
  \sigma(\vec{x},\vec{\omega}) = \sigasym.
\end{equation*}
Substituting the above into \eqref{eqn:rte} yields
\begin{align*}
    &\vec{\omega} \cdot \nabla \left[ \Lasym \right] \\
    &+ a(\vec{x}) \left[ \Lasym \right] \\
    &= b\Bigg(
      \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}')
      \left[ \Lasymp \right] \, d\vec{\omega}' \\
    &- \left[ \Lasym \right]
      \Bigg) \\
    &+ \left[ \sigasym \right].
\end{align*}
Grouping like powers of $b$, we have the decoupled set of equations
\begin{align}
  \vec{\omega} \cdot \nabla L_0(\vec{x}, \vec{\omega}) + a(\vec{x})L_0(\vec{x}) &= \sigma_0(\vec{x}, \vec{\omega}),
  \label{eqn:asymptotics_0}\\
  \vec{\omega} \cdot \nabla L_1(\vec{x}, \vec{\omega}) + a(\vec{x})L_1(\vec{x})
  &= \sigma_1(\vec{x}, \vec{\omega})
  + \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L_0(\vec{x}, \vec{\omega}')\,d\vec{\omega}' - L_0(\vec{x}, \vec{\omega}), \nonumber\\
  \vec{\omega} \cdot \nabla L_2(\vec{x}, \vec{\omega}) + a(\vec{x})L_2(\vec{x})
  &= \sigma_2(\vec{x}, \vec{\omega})
  + \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L_1(\vec{x}, \vec{\omega}')\,d\vec{\omega}' - L_1(\vec{x}, \vec{\omega}). \nonumber \\
  &\vdots \nonumber
\end{align}
In general, for $n \geq 1$,
\begin{equation}
  \vec{\omega} \cdot \nabla L_n(\vec{x}, \vec{\omega}) + a(\vec{x})L_n(\vec{x})
  = \sigma_n(\vec{x}, \vec{\omega})
  + \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}') L_{n-1}(\vec{x}, \vec{\omega}')\,d\vec{\omega}' - L_{n-1}(\vec{x}, \vec{\omega}).
  \label{eqn:asymptotics_n}
\end{equation}

For boundary conditions, let $\vec{x_s}$ be a point on the surface of the domain and $\vec{x_b}$ a point on the bottom.
Then,
\begin{equation*}
  \begin{cases}
    \Lasyms = f(\vec{\omega}) &\mbox{ if } \hat{z}\cdot\vec{\omega} > 0, \\
    \Lasymb = 0 &\mbox{ if } \hat{z}\cdot\vec{\omega} < 0.
    \end{cases}
\end{equation*}
Grouping by powers of $b$, we have
\begin{equation}
  \begin{cases}
    L_0(\vec{x_s}, \vec{\omega}) = f(\vec{\omega}) &\mbox{ if } \hat{z}\cdot\vec{\omega} > 0, \\
    L_0(\vec{x_b}, \vec{\omega}) = 0 &\mbox{ if } \hat{z}\cdot\vec{\omega} < 0,
  \end{cases}
  \label{eqn:asymptotics_bc_0}
\end{equation}
for the first term, and
\begin{equation}
  \begin{cases}
    L_n(\vec{x_s}, \vec{\omega}) = 0 &\mbox{ if } \hat{z}\cdot\vec{\omega} > 0, \\
    L_n(\vec{x_b}, \vec{\omega}) = 0 &\mbox{ if } \hat{z}\cdot\vec{\omega} < 0,
  \end{cases}
  \label{eqn:asymptotics_bc_n}
\end{equation}
for $n > 0$.

\subsection{Analytical Solution}
\label{sec:asymptotic_sol}

Given $\vec{x}, \vec{\omega}$, consider the path $\vec{l}(\vec{x}, \vec{\omega}, s)$ from \eqref{eqn:ray_path} for $s \in [0, \tilde{s}]$.
Denote the leading order radiance, absorption coefficient, and source term along the path by
\begin{align*}
  u_0(s) = L_0(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}), \\
  \tilde{a}(s) = a(\vec{l}(\vec{x}, \vec{\omega}, s)), \\
  \tilde\sigma_0(s) = \sigma_0(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}).
\end{align*}
Then, the first equation from the asymptotic expansion, \eqref{eqn:asymptotics_0} and its associated boundary condition, \eqref{eqn:asymptotics_bc_0}, can be rewritten as the first order linear ordinary differential equation and inhomogeneous boundary condition,
\begin{equation}
  \left\{
  \begin{aligned}
  \tilde\sigma_0(s) &= \frac{du_0}{ds}(s) + \tilde{a}(s) u_0(s), \\
  u_0(0) &= f(\vec{\omega})H(\vec{\omega}\cdot\hat{z}),
  \end{aligned}
  \right.
  \label{eqn:asymptotics_ode_0}
\end{equation}
where $H(x)$ is the Heaviside step function.
This equation is solved by multiplying by the appropriate integrating factor, as follows.
\begin{align*}
  \exp\left(\int_0^s \tilde{a}(s')\, ds'\right)\tilde\sigma_0(s) &= \exp\left(\int_0^s \tilde{a}(s')\, ds'\right) \frac{du_0}{ds}(s) + \exp\left(\int_0^s \tilde{a}(s')\, ds'\right) \tilde{a}(s) u_0(s) \\
  &= \frac{d}{ds}\left[\exp\left(\int_0^s \tilde{a}(s')\, ds'\right) u_0(s)\right].
\end{align*}
Then, integrating both sides yields
\begin{align*}
  \int_0^s \exp\left(\int_0^{s'} \tilde{a}(s'')\, ds''\right)\tilde\sigma_0(s')\, ds' &= \int_0^s \frac{d}{ds'}\left[\exp\left(\int_0^{s'} \tilde{a}(s'')\, ds''\right) u_0(s')\right]\, ds' \\
  &= \exp\left(\int_0^s \tilde{a}(s')\, ds'\right) u_0(s) - f(\vec{\omega})H(\vec{\omega}\cdot\hat{z}).
\end{align*}
Hence,
\begin{align}
  u_0(s) &= \left[f(\vec{\omega})H(\vec{\omega}\cdot\hat{z}) + \int_0^s \exp\left(\int_0^{s'} \tilde{a}(s'')\, ds''\right) \tilde\sigma_0(s')\, ds'\right] \exp\left(-\int_0^s \tilde{a}(s')\, ds'\right) \nonumber\\
  &= \exp\left(-\int_0^s \tilde{a}(s')\, ds'\right) f(\vec{\omega})H(\vec{\omega}\cdot\hat{z}) \nonumber\\
    &+ \int_0^s \exp\left(-\int_{s'}^s \tilde{a}(s'')\, ds''\right) \tilde\sigma_0(s')\, ds'.
  \label{eqn:asymptotics_soln_0}
\end{align}
Then, we convert back from path length $s$ to the spatial coordinate $\vec{x}$ by evaluating the one-dimensional solution at the end of the ray path.
That is,
\begin{equation*}
  L_0(\vec{x}, \vec{\omega}) = u_0(\tilde{s}).
\end{equation*}

In addition to the explicit source term, the  $n \geq 1$ equations also have a scattering term, which is an integral of the previous term in the series. The sum of these two sources is called the effective source,
\begin{equation*}
  g_n(\vec{x}, \vec{\omega}) = \sigma(\vec{x}, \vec{\omega}) + \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}')
  L_{n-1}(\vec{x}, \vec{\omega}')\,d\vec{\omega}' - L_{n-1}(\vec{x}, \vec{\omega}).
\end{equation*}
This can be similarly extracted along a ray path as
\begin{equation*}
  \tilde{g}_n(s) = \tilde{\sigma}(s) + \int_{4\pi} \beta(\vec{\omega}\cdot\vec{\omega}')
  L_{n-1}(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}')\,d\vec{\omega}' - L_{n-1}(\vec{l}(\vec{x}, \vec{\omega}, s), \vec{\omega}).
\end{equation*}
Then, since $\tilde{g}_n$ depends only on $L_{n-1}$ and is therefore independent of $u_n$, \eqref{eqn:asymptotics_n} and its boundary condition \eqref{eqn:asymptotics_bc_n} can be written as the first order linear ordinary differential equation along the ray path,
\begin{equation}
  \left\{
  \begin{aligned}
    \tilde{g}_n(s) &= \frac{du_n}{ds}(s) + \tilde{a}(s)u_n(s) \\
    u_n(0) &= 0
  \end{aligned}
  \right.
  \label{eqn:asymptotics_ode_n}
\end{equation}
This is exactly \eqref{eqn:asymptotics_ode_0} with $\tilde{\sigma}_0 \to \tilde{g}_n$ and $f(\vec{\omega}) \to 0$.
Hence,
\begin{equation}
  u_n(s) = \int_0^s\tilde{g}_n(s')\exp\left( -\int_{s'}^{s}\tilde{a}(s'')\,ds'' \right)\, ds'.
  \label{eqn:asymptotics_soln_n}
\end{equation}
As before, the conversion back to full spatial and angular coordinates is
\begin{equation*}
  L_n(\vec{x}, \vec{\omega}) = u_n(\tilde{s}).
\end{equation*}