← Back to Part 5: Clouds & Precipitation

5.2 Precipitation

Precipitation formation bridges the gap between cloud microphysics and the global hydrological cycle. Cloud droplets (~10 μm) must grow by a factor of 100 in radius to become raindrops (~1 mm), a volume increase of 10&sup6;. This enormous growth requires mechanisms beyond simple vapor diffusion: collision-coalescence in warm clouds and the Bergeron-Findeisen ice crystal process in cold clouds. Understanding these pathways is essential for weather prediction, remote sensing, and climate modeling.

Warm Rain vs Cold Rain Processes

Warm Rain Process

Occurs entirely at temperatures above 0°C. Precipitation forms exclusively through collision-coalescence. Dominant in tropical maritime clouds where cloud bases are warm, CCN concentrations are low (producing fewer but larger droplets), and cloud depths are sufficient for collector drops to grow to precipitation size.

Timescale: 20-40 minutes from nucleation to rainfall. Responsible for most tropical oceanic rain. Drop spectra tend to be broad with large median diameters.

Cold Rain (Ice) Process

Involves ice-phase microphysics via the Bergeron-Findeisen mechanism. Ice crystals grow rapidly by vapor deposition in mixed-phase clouds, then further by aggregation and riming. Snow melts to rain below the 0°C level (bright band on radar). Dominant at mid and high latitudes where cloud tops extend well above freezing.

Timescale: 30-60 minutes. Accounts for most precipitation globally (~70% by some estimates). Produces stratiform rain with smaller, more uniform drop sizes.

Collision-Coalescence: Warm Rain Process

In warm clouds (entirely above 0°C), precipitation forms exclusively through collision-coalescence. Larger drops fall faster than smaller ones, sweeping them up in their path. This process requires a broad initial drop size distribution and becomes efficient when some drops exceed ~20 μm in radius (collector drops).

Terminal Velocity

The terminal velocity of a falling drop depends on its size. Different regimes apply:

Stokes Regime (r < 30 μm)

$$v_t = \frac{2}{9}\frac{\rho_w g r^2}{\mu}$$

Viscous drag dominates. Terminal velocity proportional to r². Cloud droplets: v_t ~ 0.01-1 cm/s.

Intermediate (30-500 μm)

$$v_t \approx 8000 \, r \quad (\text{m/s, r in m})$$

Transition regime. Drizzle drops: v_t ~ 0.5-2 m/s. Empirical fits needed.

Large Drops (D > 1 mm)

$$v_t \approx 9.65 - 10.3\,e^{-0.6D}$$

Gunn-Kinzer empirical fit (D in mm, v in m/s). Raindrops: v_t ~ 4-9 m/s. Max ~9.2 m/s for D ~ 5 mm.

Collection Kernel

The rate at which a collector drop of radius R sweeps up smaller drops of radius r is described by the gravitational collection kernel:

Gravitational Collection Kernel

$$K(R, r) = \pi(R + r)^2 \, |V(R) - V(r)| \, E(R, r)$$

where V(R), V(r) are terminal velocities and E(R,r) is the collection efficiency (product of collision efficiency and coalescence efficiency).

$$E(R, r) = E_c(R, r) \times E_{\text{coal}}(R, r)$$

E_c depends on the size ratio p = r/R. For p < 0.05, collision efficiency is very low (~0.01) because small drops follow streamlines around the collector. E_c increases rapidly for p > 0.1 and approaches 1 for p > 0.6.

Stochastic Collection Equation (Smoluchowski)

The Smoluchowski coagulation equation, adapted for cloud physics as the stochastic collection equation (SCE), describes the time evolution of the drop mass distribution function n(x,t) where x is drop mass:

Stochastic Collection Equation (SCE / Smoluchowski)

$$\frac{\partial n(x,t)}{\partial t} = \frac{1}{2}\int_0^x K(x-y, y)\,n(x-y,t)\,n(y,t)\,dy - n(x,t)\int_0^\infty K(x,y)\,n(y,t)\,dy$$

The first term represents gain (two smaller drops merging to form a drop of mass x). The second term represents loss (drop of mass x collecting another drop and leaving the size bin).

Continuous Collection Model

For a single large collector drop of radius R falling through a population of small cloud droplets with liquid water content LWC, the continuous collection equation gives:

Continuous Collection Growth

$$\frac{dR}{dt} = \frac{\text{LWC} \cdot E(R) \cdot V(R)}{4 \, \rho_w}$$

Growth rate accelerates as R increases (because V(R) increases with R), leading to a positive feedback. A collector drop can grow from 20 μm to 1 mm in approximately 20-30 minutes with LWC ~ 1 g/m³.

Ice Nucleation and the Bergeron-Findeisen Process

In mixed-phase clouds, the Bergeron-Findeisen process exploits the saturation vapor pressure difference over ice versus liquid water at subfreezing temperatures.

Homogeneous Nucleation

Spontaneous freezing of pure water droplets occurs below -38 to -40°C. Requires no foreign particles. Dominates in the upper troposphere (cirrus formation).

Heterogeneous Nucleation

Ice nucleating particles (INPs) catalyze freezing at -10 to -35°C. Effective INPs include mineral dust (feldspar), biological particles, and volcanic ash. INP concentration is typically 0.01-10 per liter.

Bergeron-Findeisen: Vapor Pressure Difference

$$e_{s,\text{water}}(T) > e_{s,\text{ice}}(T) \quad \text{for } T < 0°\text{C}$$

The difference peaks near -12°C at about 0.27 hPa (~10% of e_s,ice). In a cloud saturated with respect to water, the environment is supersaturated with respect to ice by 10-15%, driving rapid ice crystal growth.

Ice Crystal Growth Rate

The diffusional growth rate of an ice crystal in a supersaturated environment follows:

$$\frac{dm}{dt} = 4\pi C \, D_v \, (S_i - 1) \, \rho_{s,i}(T) \cdot f_v$$

where C is the crystal capacitance (shape-dependent), D_v is vapor diffusivity, S_i is ice supersaturation, and f_v is the ventilation factor. The capacitance C equals the radius for spheres but differs for plates, columns, and dendrites.

Ice Crystal Habits

Temperature Range	Primary Habit	Growth Axis
0 to -4°C	Thin plates, dendrites	a-axis (basal)
-4 to -10°C	Columns, needles	c-axis (prism)
-10 to -22°C	Plates, sectors, dendrites	a-axis (basal)
-22 to -40°C	Columns, thick plates	c-axis (prism)
< -40°C	Polycrystals, bullet rosettes	Mixed

Raindrop Size Distribution and Radar Reflectivity

The size distribution of raindrops determines rainfall rate, radar reflectivity, and kinetic energy. The Marshall-Palmer (1948) exponential distribution provides a good fit for stratiform rain:

Marshall-Palmer Distribution

$$N(D) = N_0 \exp(-\Lambda D)$$

N_0 = 8000 m³ mm¹ (intercept) and Λ = 4.1 R^(-0.21) mm¹ (slope parameter, R in mm/hr). D is drop diameter in mm.

Gamma Distribution (Generalized Form)

The three-parameter gamma distribution generalizes the exponential form and better captures observed curvature in measured DSDs, especially for convective rain:

$$N(D) = N_0 \, D^{\mu} \exp(-\Lambda D)$$

The shape parameter μ controls curvature: μ = 0 recovers Marshall-Palmer, μ > 0 produces concave-up DSDs (deficit of small and large drops), μ < 0 produces concave-down DSDs. Typical values: μ = 2-6 for convective rain.

Rainfall Rate and Radar Reflectivity

Rainfall Rate

$$R = \frac{\pi}{6}\int_0^\infty D^3 \, v(D) \, N(D) \, dD$$

Radar Reflectivity Factor

$$Z = \int_0^\infty D^6 \, N(D) \, dD$$

Z-R Relationships

The power-law Z = aR<sup>b</sup> is used by operational radar for rainfall estimation:

$$Z = a \, R^b$$

Rain Type	a	b	Reference
Marshall-Palmer (stratiform)	200	1.6	WSR-88D default
Convective rain	300	1.4	Joss & Waldvogel
Tropical rain	250	1.2	Rosenfeld et al.
Snow (liquid equiv.)	75	2.0	Sekhon & Srivastava

Fortran: Stochastic Collection Equation Solver

This Fortran program solves the discretized stochastic collection equation using the Berry-Reinhardt method, tracking the time evolution of the drop spectrum from an initial cloud droplet distribution through to precipitation-sized drops.

program stochastic_collection
  ! ============================================================
  ! Solve the discretized Stochastic Collection Equation (SCE)
  ! for the time evolution of a cloud drop spectrum
  ! Uses logarithmic mass grid and Golovin kernel for validation
  ! ============================================================
  implicit none

  integer, parameter :: nbins = 80
  real(8), parameter :: pi = 3.14159265358979d0
  real(8), parameter :: rho_w = 1000.0d0
  real(8), parameter :: dt = 1.0d0        ! Time step (s)
  integer, parameter :: nsteps = 1800     ! 30 minutes

  real(8) :: mass(nbins)       ! Bin center mass (kg)
  real(8) :: radius(nbins)     ! Bin center radius (m)
  real(8) :: n(nbins)          ! Number concentration (#/m^3/bin)
  real(8) :: n_new(nbins)      ! Updated spectrum
  real(8) :: K_ij              ! Collection kernel
  real(8) :: dm, m_sum, LWC, Ntot, Dmean
  real(8) :: m_min, m_ratio
  integer :: i, j, k, istep

  ! --- Set up logarithmic mass grid ---
  m_min = (4.0d0/3.0d0) * pi * rho_w * (5.0d-6)**3  ! 5 um drop
  m_ratio = 2.0d0**(1.0d0/3.0d0)   ! Mass doubling every 3 bins

  do i = 1, nbins
    mass(i) = m_min * m_ratio**(i - 1)
    radius(i) = ((3.0d0 * mass(i)) / (4.0d0 * pi * rho_w))**(1.0d0/3.0d0)
  end do

  ! --- Initial cloud droplet distribution (lognormal-like) ---
  do i = 1, nbins
    n(i) = 300.0d6 * exp(-0.5d0 * ((log(radius(i)/10.0d-6))/0.3d0)**2)
    n(i) = n(i) * (mass(i) * log(m_ratio))  ! Convert to per-bin
  end do

  ! --- Time integration loop ---
  open(unit=10, file='sce_output.dat', status='replace')
  write(10,'(A)') '# time(s)  LWC(g/m3)  Ntot(#/cm3)  Dmean(um)'

  do istep = 1, nsteps
    n_new = n

    ! --- Collection: gain and loss terms ---
    do i = 1, nbins
      do j = 1, i
        ! Gravitational collection kernel (simplified)
        K_ij = pi * (radius(i) + radius(j))**2  &
             * abs(terminal_vel(radius(i)) - terminal_vel(radius(j))) &
             * 0.7d0   ! Mean collection efficiency

        ! Find target bin for merged drop
        m_sum = mass(i) + mass(j)
        do k = 1, nbins - 1
          if (mass(k) <= m_sum .and. m_sum < mass(k+1)) exit
        end do
        if (k >= nbins) k = nbins

        ! Loss from bins i and j
        dm = K_ij * n(i) * n(j) * dt
        if (i == j) dm = 0.5d0 * dm  ! Self-collection factor

        n_new(i) = n_new(i) - dm
        if (i /= j) n_new(j) = n_new(j) - dm

        ! Gain in target bin k
        n_new(k) = n_new(k) + dm
      end do
    end do

    ! Enforce non-negativity
    do i = 1, nbins
      if (n_new(i) < 0.0d0) n_new(i) = 0.0d0
    end do
    n = n_new

    ! --- Diagnostics every 60 seconds ---
    if (mod(istep, 60) == 0) then
      LWC = 0.0d0; Ntot = 0.0d0; Dmean = 0.0d0
      do i = 1, nbins
        LWC  = LWC  + n(i) * mass(i)
        Ntot = Ntot + n(i)
        Dmean = Dmean + n(i) * 2.0d0 * radius(i)
      end do
      if (Ntot > 0.0d0) Dmean = Dmean / Ntot
      write(10,'(F8.0, 3E14.5)') dble(istep), LWC*1.0d3, Ntot*1.0d-6, Dmean*1.0d6
      write(*,'(A,F6.0,A,F8.4,A,F10.1,A,F8.1)') &
        't=', dble(istep), 's  LWC=', LWC*1d3, ' g/m3  N=', Ntot*1d-6, '/cm3  Dm=', Dmean*1d6
    end if
  end do

  close(10)
  write(*,'(A)') 'SCE integration complete -> sce_output.dat'

contains

  function terminal_vel(r) result(vt)
    real(8), intent(in) :: r
    real(8) :: vt, D_mm
    D_mm = 2.0d0 * r * 1.0d3
    if (D_mm < 0.06d0) then
      vt = 2.0d0 * rho_w * 9.81d0 * r**2 / (9.0d0 * 1.8d-5)
    else if (D_mm < 1.0d0) then
      vt = 4.5d0 * sqrt(D_mm)
    else
      vt = 9.65d0 - 10.3d0 * exp(-0.6d0 * D_mm)
    end if
  end function terminal_vel

end program stochastic_collection

Interactive Simulation: Raindrop Size Distribution

Python

Plot Marshall-Palmer exponential drop size distributions for different rainfall rates and compute radar reflectivity Z. Shows the Z-R relationship used by operational weather radars.

raindrop_size_distribution.py88 lines

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma as gamma_func

# Marshall-Palmer DSD: N(D) = N0 * exp(-Lambda * D)
N0 = 8000.0  # m^-3 mm^-1 (intercept parameter)

def marshall_palmer(D, R):
    """N(D) for rainfall rate R (mm/hr). D in mm."""
    Lambda = 4.1 * R**(-0.21)  # mm^-1
    return N0 * np.exp(-Lambda * D)

def terminal_velocity(D_mm):
    """Terminal velocity (m/s) for diameter D in mm."""
    return np.where(D_mm < 1.0, 4.5 * np.sqrt(np.maximum(D_mm, 0.01)),
                    9.65 - 10.3 * np.exp(-0.6 * D_mm))

# Diameter array
D = np.linspace(0.1, 7.0, 500)
dD = D[1] - D[0]

# Rainfall rates to plot
rain_rates = [1, 5, 10, 50]
colors = ['#38bdf8', '#22d3ee', '#2dd4bf', '#a78bfa']

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)

# Panel 1: Marshall-Palmer DSD for different rain rates
print("Raindrop Size Distribution Analysis")
print("=" * 55)
print(f"{'R (mm/hr)':>10} {'Lambda':>10} {'Z (dBZ)':>10} {'Z-R check':>12}")
print("-" * 55)

Z_values = []
for R, c in zip(rain_rates, colors):
    ND = marshall_palmer(D, R)
    # Compute radar reflectivity Z = integral(D^6 * N(D) dD)
    Z = np.sum(ND * D**6 * dD)  # mm^6 m^-3
    Z_dBZ = 10 * np.log10(Z)
    Z_values.append(Z_dBZ)
    Lambda = 4.1 * R**(-0.21)
    # Z-R relationship check: Z = 200 * R^1.6
    Z_check = 200 * R**1.6
    print(f"{R:>10} {Lambda:>10.2f} {Z_dBZ:>10.1f} {10*np.log10(Z_check):>12.1f}")
    axes[0].semilogy(D, ND, color=c, lw=2.5,
                     label=f'R={R} mm/hr (Z={Z_dBZ:.1f} dBZ)')

axes[0].set_xlabel('Diameter D (mm)', color='#cbd5e1', fontsize=12)
axes[0].set_ylabel('N(D) (m$^{-3}$ mm$^{-1}$)', color='#cbd5e1', fontsize=12)
axes[0].set_title('Marshall-Palmer DSD: N(D) = N\u2080 exp(-\u039bD)', color='#cbd5e1', fontsize=13)
axes[0].set_ylim(1e-1, 1e5)
axes[0].legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

# Panel 2: Z-R relationship
R_arr = np.linspace(0.5, 80, 300)
Z_mp = 200 * R_arr**1.6       # Marshall-Palmer (stratiform)
Z_conv = 300 * R_arr**1.4     # Convective
Z_trop = 250 * R_arr**1.2     # Tropical

axes[1].semilogy(R_arr, Z_mp, color='#38bdf8', lw=2.5, label='Z=200R\u00b9\u22c5\u2076 (stratiform)')
axes[1].semilogy(R_arr, Z_conv, color='#f87171', lw=2.5, ls='--', label='Z=300R\u00b9\u22c5\u2074 (convective)')
axes[1].semilogy(R_arr, Z_trop, color='#a78bfa', lw=2.5, ls=':', label='Z=250R\u00b9\u22c5\u00b2 (tropical)')

# Mark typical thresholds
for z_val, lbl in [(200, 'Light'), (2000, 'Moderate'), (20000, 'Heavy')]:
    axes[1].axhline(y=z_val, color='#475569', ls='--', alpha=0.4)
    axes[1].text(75, z_val * 1.3, lbl, color='#94a3b8', fontsize=8)

axes[1].set_xlabel('Rainfall Rate R (mm/hr)', color='#cbd5e1', fontsize=12)
axes[1].set_ylabel('Reflectivity Z (mm\u2076/m\u00b3)', color='#cbd5e1', fontsize=12)
axes[1].set_title('Z-R Relationships', color='#cbd5e1', fontsize=13)
axes[1].legend(fontsize=9, facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor=fig.get_facecolor())
print()
print("Z-R relationships: Z = aR^b")
print("  Stratiform:  a=200, b=1.6 (WSR-88D default)")
print("  Convective:  a=300, b=1.4")
print("  Tropical:    a=250, b=1.2")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← 5.1 Cloud Physics 5.3 Severe Weather →