Ash Plume Light Extinction and Scattering in GLSL

This article details the mathematical formulation and GLSL implementation of light extinction and scattering through dense volcanic ash plumes. Volcanic ash presents unique optical properties compared to standard water-droplet clouds, characterized by high optical depth, strong absorption from pulverized rock and glass, and heavy forward-to-backward anisotropic scattering. Below, we break down the Radiative Transfer Equation (RTE) adapted for dense particulate media, discuss the Beer-Lambert law and empirical multiple scattering models, and provide an optimized GLSL raymarching implementation.

Radiative Transfer and Extinction Formulation

Light traveling through a participating medium like a volcanic ash plume loses energy via extinction (absorption and out-scattering) and gains energy via in-scattering. The change in radiance \(L\) at position \(\mathbf{x}\) along direction \(\mathbf{\omega}\) is defined by the Radiative Transfer Equation:

\[\frac{dL(\mathbf{x}, \mathbf{\omega})}{ds} = -\mu_e(\mathbf{x}) L(\mathbf{x}, \mathbf{\omega}) + \mu_s(\mathbf{x}) L_i(\mathbf{x}, \mathbf{\omega})\]

Where:

  • \(\mu_e = \mu_a + \mu_s\): The extinction coefficient, composed of the absorption coefficient \(\mu_a\) and the scattering coefficient \(\mu_s\).
  • Albedo (\(\alpha = \frac{\mu_s}{\mu_e}\)): Volcanic ash has a significantly lower single-scattering albedo (\(\alpha \approx 0.1 - 0.3\)) than water clouds (\(\alpha \approx 0.99\)), causing rapid darkening and light attenuation.

Transmittance (Beer-Lambert Law)

The fractional light transmitted between two points over distance \(d\) through a homogeneous segment is calculated using the Beer-Lambert Law:

\[T = e^{-\int_0^d \mu_e(s) \, ds} \approx e^{-\mu_e \cdot \rho \cdot d}\]

Where \(\rho\) is the sampled particulate density at step \(s\).

Phase Function (Henyey-Greenstein)

The directional distribution of scattered light is driven by the phase function \(P(\theta)\). For volcanic ash, a two-lobe Henyey-Greenstein phase function accounts for forward diffraction around ash particles and backward reflection off jagged mineral surfaces:

\[P(\theta, g) = \frac{1 - g^2}{4\pi (1 + g^2 - 2g \cos\theta)^{3/2}}\]

\[P_{\text{ash}}(\theta) = k \cdot P(\theta, g_{\text{forward}}) + (1 - k) \cdot P(\theta, g_{\text{backward}})\]

Where \(g \in (-1, 1)\) represents the asymmetry parameter, with \(g_{\text{forward}} \approx 0.6\) to \(0.8\), \(g_{\text{backward}} \approx -0.2\), and weight \(k \approx 0.7\).

Dense Particulate Multiple Scattering Approximation

In dense ash, multiple scattering normally requires expensive path tracing. In real-time GLSL rendering, this is approximated using a combination of:

  1. The "Powder Effect": Darkens edges when looking through dense layers away from light: \[S_{\text{powder}} = 1.0 - e^{-2 \cdot \rho \cdot d}\]
  2. Multi-Octave Scattering: Blurs out-scattering into lower extinction tiers to mimic forward multiple-scattered light penetration.

GLSL Implementation

The following GLSL fragment shader functions demonstrate how to calculate light extinction, in-scattering, and direct transmittance along a view ray march through a density field.

#version 450 core

// Material Constants for Volcanic Ash
const float ABSORPTION_COEFF = 0.85;  // Strong absorption (dark silicates/basalt)
const float SCATTERING_COEFF = 0.35;  // Single scattering
const float EXTINCTION_COEFF = ABSORPTION_COEFF + SCATTERING_COEFF;
const float ALBEDO           = SCATTERING_COEFF / EXTINCTION_COEFF;

// Henyey-Greenstein Phase Function
float PhaseHenyeyGreenstein(float cosTheta, float g) {
    float g2 = g * g;
    return (1.0 - g2) / (4.0 * 3.14159265 * pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5));
}

// Dual-Lobe Phase Function for Ash Particulates
float AshPhase(float cosTheta) {
    float forward  = PhaseHenyeyGreenstein(cosTheta,  0.65);
    float backward = PhaseHenyeyGreenstein(cosTheta, -0.25);
    return mix(backward, forward, 0.75);
}

// Powder effect to approximate dense particulate back-scattering
float PowderEffect(float density, float stepSize) {
    return 1.0 - exp(-density * stepSize * 2.0);
}

// Sample density field (e.g., 3D procedural noise or texture)
float SampleAshDensity(vec3 pos);

// Shadow Ray: Computes optical depth toward the light source
float LightTransmittance(vec3 samplePos, vec3 lightDir, float shadowStepSize, int shadowSteps) {
    float opticalDepth = 0.0;
    
    for (int i = 0; i < shadowSteps; ++i) {
        samplePos += lightDir * shadowStepSize;
        float d = SampleAshDensity(samplePos);
        opticalDepth += d * EXTINCTION_COEFF * shadowStepSize;
    }
    
    return exp(-opticalDepth);
}

// View Raymarching Function
vec4 RenderAshPlume(vec3 rayOrigin, vec3 rayDir, vec3 lightDir, vec3 lightColor, float maxDist) {
    const int STEPS = 64;
    const int SHADOW_STEPS = 6;
    float stepSize = maxDist / float(STEPS);
    float shadowStepSize = stepSize * 1.5;
    
    vec3 currentPos = rayOrigin;
    float transmittance = 1.0;
    vec3 totalScattering = vec3(0.0);
    
    float cosTheta = dot(rayDir, lightDir);
    float phase = AshPhase(cosTheta);
    
    for (int i = 0; i < STEPS; ++i) {
        float density = SampleAshDensity(currentPos);
        
        if (density > 0.001) {
            // Local extinction step
            float stepExtinction = density * EXTINCTION_COEFF * stepSize;
            float stepTransmittance = exp(-stepExtinction);
            
            // Light extinction from the primary light source
            float lightAtten = LightTransmittance(currentPos, lightDir, shadowStepSize, SHADOW_STEPS);
            
            // Direct in-scattering with powder effect for high density
            float powder = PowderEffect(density, stepSize);
            vec3 inScatteredLight = lightColor * lightAtten * phase * powder * ALBEDO;
            
            // Energy-conserving integration over current step
            vec3 stepScattering = inScatteredLight * (1.0 - stepTransmittance);
            totalScattering += transmittance * stepScattering;
            
            // Accumulate transmittance loss along view ray
            transmittance *= stepTransmittance;
            
            // Early exit if the medium becomes completely opaque
            if (transmittance < 0.005) {
                break;
            }
        }
        
        currentPos += rayDir * stepSize;
    }
    
    return vec4(totalScattering, 1.0 - transmittance);
}

Integration Details

  • Step Sizing: Due to steep extinction gradients in volcanic plumes, uniform step sizes often cause visual banding. Applying a small randomized offset (dither) to the ray starting point rayOrigin += rayDir * stepSize * InterleavedGradientNoise() mitigates stepping artifacts.
  • Ambient Term: Ash clouds block nearly all indirect sky light. Multiply ambient contributions by the integrated transmittance to prevent flat ambient colors deep within the inner core of the plume.