The AOMC Model Implementation

The AOMC model uses the Monte Carlo method to simulate the individual life histories of a large number of photons. By tracking the fate of each photon, the model reconstructs the light field (e.g., radiance and irradiance) throughout the water body.

The simulation is driven by the main program mc.f90, which iterates through a predefined number of photons for each specified wavelength.

Main Program (mc.f90) Structure

The mc.f90 program is the central orchestrator of the AOMC model. Its structure is divided into several logical blocks, each responsible for a specific phase of the simulation:

  1. Initialization: Sets up initial variables and seeds the random number generator.
  2. Reading General Parameters: Reads primary simulation parameters from the amc.inp file.
  3. Allocation of Memory: Dynamically allocates memory for arrays and variables based on the input parameters.
  4. Reading of Input Files: Reads data from various input files, such as conc.inp, abs.inp, scat.inp, spf.inp, lambbot.inp, and difcol.inp or skydist.inp.
  5. Calculation of IOPs: Computes inherent optical properties like beam attenuation and scattering albedo.
  6. Run the Model: Executes the core Monte Carlo simulation loop, tracing individual photon paths.
  7. Compute AOPs: Processes the tallied photon data to calculate apparent optical properties.
  8. Free Up Allocated Memory: Deallocates all dynamically allocated memory.

Photon Life Cycle: Initiation and Termination

The core of the simulation is a loop where each iteration represents the entire life of a single photon. A new photon is initiated at the start of each iteration, and its path is traced until it is terminated.

  • Initiation: A new photon’s life begins inside the INTENSITY_LOOP in mc.f90. The light subroutine is called to define its initial properties (e.g., direction), which are determined by the sky radiance distribution. This can be a simple model with direct and diffuse components (difcol.inp) or a tabulated distribution from a file (skydist.inp). The photon is placed at the air-water interface.

  • Termination: A photon’s path ends, and the simulation for that photon ceases, when one of the following events occurs, as handled by the water subroutine in water.f90:

    1. It is absorbed within the water column. This is determined by comparing a random number to the single-scattering albedo.
    2. It scatters upward and passes back out through the air-water interface.
    3. It hits a side boundary of the model domain.
    4. It is absorbed by the bottom boundary.

Once a photon is terminated, the program begins the next iteration of the loop, initiating a new photon.

Photon Path Simulation

The Air-Water Interface

  • When the photon strikes the surface, Fresnel’s equation is used with a random number to determine if it is reflected or transmitted into the water.
  • If the photon is transmitted, Snell’s Law is used to calculate its new, refracted angle of travel within the water.

Path Length Determination

  • The distance the photon travels between interactions is its optical path length. This step size, \(\Delta l\), is determined probabilistically using the beam attenuation coefficient \(c(\lambda)\) and a random number, \(\xi_1\), uniformly distributed between 0 and 1: \[ \Delta l = -\frac{\ln(\xi_1)}{c(\lambda)} \]

Absorption or Scattering

  • After traveling the distance \(\Delta l\), the photon interacts with the medium. A decision is made to either absorb or scatter the photon.
  • A new random number, \(\xi_2\), is generated and compared to the single-scattering albedo, \(\omega_0(\lambda)\).
    • If \(\xi_2 \le \omega_0(\lambda)\), the photon is scattered.
    • If \(\xi_2 > \omega_0(\lambda)\), the photon is absorbed, and its simulation is terminated.

Scattering Process

If the photon is scattered, its direction must change.

  1. Select Constituent: The model first determines which constituent particle scattered the photon. This is a random choice weighted by the relative scattering contribution of each constituent at that wavelength.
  2. Select New Angle: The model then uses the cumulative scattering phase function for that specific particle (from spf.inp) as a probability table. A new random number, \(\xi_3\), is generated to select a new polar scattering angle, \(\theta\). A uniformly distributed random number determines the azimuthal scattering angle, \(\phi\).

For a detailed review of the scattering phase function see here

Boundary Interactions

  • If a photon strikes the bottom boundary, a random number is compared to the bottom reflectance (\(R_b\)) from lambbot.inp.
    • If absorbed, the photon is terminated.
    • If reflected, another random number determines if the reflection is specular (mirror-like) or diffuse (Lambertian), and the photon’s journey continues in an upward direction.

Photon Logging and AOP Calculation

As photons travel through the simulated environment, their crossings of predefined horizontal layers are recorded. This information is used to calculate the Apparent Optical Properties (AOPs).

Logging Photon Crossings

  • The logbin subroutine (in logbin.f90) is called whenever a photon crosses a logging layer.
  • logbin records the event by incrementing a counter in a large multi-dimensional array n(lambda, i, j, k, t, p). This array acts as a histogram, binning photons by:
    • lambda: Wavelength.
    • i, j: Horizontal grid cell location.
    • k: Vertical layer index.
    • t, p: Bins for the polar (\(\theta\)) and azimuthal (\(\phi\)) angles of the photon’s direction.

Polar Angle (Theta) Binning

The AOMC model supports two different methods for defining the polar angle bins, controlled by the angint parameter in amc.inp:

  1. Equal Angle Intervals (angint = 1): The polar angle range [0, \(\pi\)] is divided into alphaint bins of equal angular width.
  2. Equal Cosine Intervals (angint = 0): The cosine of the polar angle range [1, -1] is divided into alphaint equal intervals. This provides higher angular resolution near the horizontal plane (\(\theta = \pi/2\)) compared to the equal angle method.

Conversion to AOPs

  • After all photon simulations are complete, the n array contains the raw counts of photon crossings.
  • In the final section of mc.f90, these counts are converted into physical quantities:
    • Radiance (\(L\)): The number of photons in each angular bin is normalized by the solid angle of that bin to calculate radiance. \[ L = \frac{N_{photons}}{\Delta\Omega} \]
    • Irradiance (\(E_d, E_u\)): Planar irradiances are calculated by summing the radiance over the downward or upward hemisphere, respectively, and weighting by the cosine of the polar angle. Scalar irradiances (\(E_{od}, E_{ou}\)) are also computed by summing radiance over the hemispheres without the cosine weighting.

This process transforms the statistical data from the Monte Carlo simulation into meaningful optical properties of the light field in the water.

The following table summarizes the conversion of the tallied photons in the AOMC model to the final radiometric quantities. \(N(\theta,\phi)\) represents the number of photon crossings, across a horizontal plane, tallied within the angular bin corresponding to polar angle \(\theta\) and azimuth angle \(\phi\).

Radiometric Quantity Conventional Definition Derivation from the AOMC model
Field Radiance, \[L\] \[\frac{d^2\Phi}{dA |\cos(\theta)| d\omega}\] \[\frac{N(\theta,\phi)}{|\cos(\theta)|\sin(\theta)d\theta d\phi}\]
Downwelling Planar Irradiance \[E_d\] \[\int_0^{2\pi}\int_0^{\pi/2} L(\theta,\phi)\cos(\theta)\sin(\theta)d\theta d\phi\] \[\sum_0^{2\pi}\sum_0^{\pi/2} N(\theta,\phi)\]
Upwelling Planar Irradiance, \[E_u\] \[-\int_0^{2\pi}\int_{\pi/2}^{\pi} L(\theta,\phi)\cos(\theta)\sin(\theta)d\theta d\phi\] \[\sum_0^{2\pi}\sum_{\pi/2}^{\pi} N(\theta,\phi)\]
Downwelling Scalar Irradiance, \[E_od\] \[\int_0^{2\pi}\int_0^{\pi/2} L(\theta,\phi)\sin(\theta)d\theta d\phi\] \[\sum_0^{2\pi}\sum_0^{\pi/2} \frac{N(\theta,\phi)}{|\cos(\theta)|}\]
Upwelling Scalar Irradiance, \[E_ou\] \[\int_0^{2\pi}\int_{\pi/2}^{\pi} L(\theta,\phi)\sin(\theta)d\theta d\phi\] \[\sum_0^{2\pi}\sum_{\pi/2}^{\pi} \frac{N(\theta,\phi)}{|\cos(\theta)|}\]
Polar Radiance, \[L_d, L_u\] \[\frac{d^2\Phi}{dA |\cos(\theta)| d\omega}, at\ \theta=0, \pi\] \[\frac{\sum_{\phi=0}^{2\pi} N_{cap}(\phi)}{\int_0^{2\pi}\int_{cap} |\cos(\theta)|\sin(\theta)d\theta d\phi}\]

The Implicit Cosine Correction in Photon Tallies

It may seem counter-intuitive that the final calculation for planar irradiance (\(E_d\) and \(E_u\)) is a direct sum of the tallied photon counts \(N(\theta, \phi)\), while the physical definition of irradiance explicitly includes a cosine factor. This is not an oversight, but a fundamental and elegant property of how photon transport is modeled across horizontal layers in a Monte Carlo simulation. The key is to recognize that the probability of a photon successfully crossing a horizontal detection plane is naturally proportional to the cosine of its zenith angle. An obliquely-traveling photon has a longer path to travel to cross a given vertical distance compared to a photon traveling straight down. Because it travels farther through the water, it has a higher chance of being absorbed or scattered away before it can be tallied.

This “survival probability” acts as an inherent, or implicit, cosine weighting on the photon population that reaches the detector. The act of simply summing the weights of the photons that successfully cross the plane has already accounted for the \(\cos(\theta)\) factor required for planar irradiance. This is also why the scalar irradiance calculation, which sums energy from all directions without a cosine weighting, must computationally remove this implicit factor by dividing the tallied counts \(N(\theta, \phi)\) by \(|\cos(\theta)|\).

Why AOMC Uses a \(1/cos(\theta)\) Weighting for Incident Photons

When photons are launched toward the water surface in the model, each photon represents a small parcel of radiant energy. However, photons arriving at different zenith angles do not bring the same amount of energy per unit horizontal area. A nadir photon (\(\theta=0^\circ\)) delivers its energy directly onto the surface, while an oblique photon (\(\theta=80^\circ\)) spreads the same energy across an area that is \(1/cos(\theta)\) times larger. Physically, this means that steep‑angle photons contribute far less downward irradiance than near‑vertical ones.

Below the water surface this geometric effect is automatically handled by the model: only photons whose paths intersect a horizontal plane are counted, and the probability of crossing that plane already carries the appropriate \(cos(\theta)\) dependence. At the surface, however, photons are forced to arrive at some (\(x,y\)) location regardless of angle, so this natural geometric filtering does not occur. Consequently, the model must explicitly impose the projection factor that relates radiance to irradiance.

The model does this by accumulating a normalization factor

\[ cositer(\lambda) = \sum{1/cos(\theta)} \]

taken over all photons incident on the surface. Although this looks like it “upweights” low incidence photons, it actually appears in the denominator of the final planar and scalar irradiance calculations. Thus, when many photons arrive at steep angles (large \(1/cos(\theta)\)), the normalization constant increases and the computed irradiances decrease correctly reflecting the fact that oblique light delivers less energy into the water column.

In short, the \(1/cos(\theta)\) weighting ensures that the model faithfully reconstructs the physically correct irradiance at the water surface. Without this correction, every incoming photon would carry identical weight, and the model would be unable to distinguish between strong near‑vertical illumination and weak grazing illumination leading to significant errors in computed underwater irradiances.