MADIS · Blizzard Jonas · 2016-01-23 1200Z

The Dense Mesonet as a
CAD-Aware Atmospheric MRI

CPA Weather Lab · Camp Hill PA · May 2026 · 1,676 PA/Mid-Atlantic Stations

An MRI machine does not sense the body at the slice it images — it reconstructs internal structure from a dense field of sensors at the surface. A sufficiently dense mesonet does the same thing for the atmosphere. Every station is a transducer. The network as a whole is a tomographic instrument. The equations below are the reconstruction kernel.

📡 Network Snapshot

Global records152,729

PA/Mid-Atl records14,722

Unique stations1,676

Mean spacing~6 km

Providers17 types

🌡 Observed Ranges

Temperature−17.8 → +35.6°C

θ (pot. temp)255 → 311 K

Wind speed0 → 25 m/s

Altimeter871 → 1047 hPa

PrecipAccum max71.1 mm

📊 Variable Availability

Temperature91.2%

Wind speed87.2%

Altimeter75.4%

Dew point58.4%

Solar radiation13.9%

255–311 Kθ range PA domain

−17.8°CT_w coldest Jonas

0–401 J/m³KE density range

6 kmMean station spacing

28Derivable quantities

5-layerMRI scan protocol

FILTER:

TIER 1

Single-Station Thermodynamic Derivations

These require only the local observation vector. Every station can compute them independently without knowledge of any neighbor.

2.1 Saturation Vapor Pressure e_s Thermodynamic ▼

The August-Roche-Magnus approximation — accurate to ±0.1% from −40°C to +40°C. This is the absolute ceiling on moisture content at a given temperature: the atmosphere cannot hold more than e_s without condensation.

e_s = 6.112 × exp( 17.67 × T_C / (T_C + 243.5) ) [hPa] # Python: e_s = 6.112 * np.exp(17.67 * T_C / (T_C + 243.5))

Jonas observed (PA): T ranged −17.8°C to +35.6°C → e_s = 1.52 to 58.1 hPa (39× range). This gradient controls the rain/snow boundary and the rate of cold dome erosion via diabatic warming.

T (air temp) — 91.2%

2.2 Actual Vapor Pressure e Thermodynamic ▼

The actual partial pressure of water vapor in the air. Can be derived from dew point (preferred) or from RH + e_s.

# From dew point (preferred): e = 6.112 × exp( 17.67 × T_d,C / (T_d,C + 243.5) ) [hPa] # From relative humidity: e = e_s × (RH / 100) [hPa]

Jonas range (PA): 0.07 – 10.45 hPa. Dew point available at 44% of PA subset stations; RH at 86% — prefer T_d when available for accuracy.

T_d (dew point) — 58.4% RH — 84.7%

2.3 Mixing Ratio w — The Conserved Moisture Tracer Thermodynamic ▼

Mass of water vapor per unit mass of dry air. Conserved in adiabatic processes — the fingerprint of an air parcel's origin. You can trace a low-level jet's Gulf moisture signature northward even when the air has cooled 20°C.

w = 0.622 × e / (P − e) [kg/kg] # Typical winter PA values: 0.001–0.006 kg/kg (1–6 g/kg) # Jonas PA range: 0.0000 – 0.0065 kg/kg

Physical meaning: The 1000:1 contrast between mountain dry air and coastal moisture sources drives the moisture flux divergence that determines snowfall banding. w is the tracer that reveals airmass origin — a cold dome injected from Canada will have w < 0.002 kg/kg; marine air will carry w > 0.005 kg/kg.

T_d (dew point) — 58.4% P (altimeter) — 75.4%

2.4 Vapor Pressure Deficit VPD Thermodynamic Cold Dome ▼

The "dryness" of the air — how far the atmosphere is from saturation. VPD drives evaporation, controls cloud formation, and diagnoses the cold dome interior.

VPD = e_s − e [hPa] # VPD < 0: supersaturation → fog, ice crystal formation, heavy snow # VPD ≈ 0: saturated cold dome core # VPD > 3: dry slot, storm weakening

Jonas range (PA): −0.75 to +7.4 hPa. Negative VPD at the core of the cold dome marks active ice crystal nucleation. High VPD on the warm sector side marks the dry slot intruding from the southwest.

T_d — 58.4% T — 91.2%

2.5 Virtual Temperature T_v Thermodynamic ▼

The temperature dry air would need to have the same density as moist air. Thermodynamically correct for all buoyancy and pressure calculations. Must be used when computing density-derived quantities.

T_v = T × (1 + 0.61 × w) [K] # Correction is only 0.1–0.4 K in winter, but critical for # divergence and density calculations

Jonas range (PA): 255.5 – 282.0 K.

T — 91.2% w (needs T_d) — 58.4%

2.6 Potential Temperature θ — The Airmass Label Thermodynamic Cold Dome ▼

The temperature a parcel would have if brought adiabatically to 1000 hPa. Conserved in dry adiabatic processes. Plotting θ instead of T removes topographic contamination — this is the fundamental "tissue type" label in the atmospheric MRI.

θ = T × (1000 / P)^0.286 [K] # Exponent = R_d / c_p = 287 / 1004 = 0.286 # Jonas PA cold dome wall: ~256 K (cold) vs ~268 K (warm sector) # 12 K θ contrast = the cold dome boundary

Jonas range (PA): 255.3 – 311.0 K. The cold dome appears as a sharp θ gradient parallel to the Appalachian ridges — invisible in raw T due to elevation contamination.

T — 91.2% P (altimeter) — 75.4%

2.7 Equivalent Potential Temperature θ_e — The MRI Contrast Agent Thermodynamic Precipitation ▼

Adds latent heat to θ — conserved in moist adiabatic processes and through precipitation. Fronts and convective instability are visible in θ_e even when T gradients are obscured by terrain. The θ_e boundary between the cold dome and warm sector defined the Jonas rain/snow line with far greater precision than any temperature threshold.

θ_e ≈ θ × exp( L_v × w / (c_p × T) ) [K] # L_v = 2.5×10⁶ J/kg (latent heat of vaporization) # c_p = 1004 J/kg/K (specific heat dry air) # Jonas range (PA): 257.7 – 295.5 K # Layer where dθ_e/dz < 0 → convective instability → thunder-snow

Jonas PA range: 257.7 – 295.5 K. The θ_e gradient identified the rain/snow line 2+ hours before T_w crossed 0°C at individual stations.

θ (needs T, P) w (needs T_d) — 58.4%

2.8 Air Density ρ — The Tissue Density Thermodynamic ▼

The mass of air per unit volume — the "tissue density" in the atmospheric MRI analogy. The cold dome is literally heavier air. Mapping ρ from station altimeter readings gives the mass distribution of the atmosphere.

ρ = P / (R_d × T_v) [kg/m³] # R_d = 287.05 J/kg/K # Jonas range (PA): 1.099 – 1.369 kg/m³ # High-elevation stations: ρ ~10–15% lower than valley stations

Jonas range (PA): 1.099 – 1.369 kg/m³. The cold, dense low-elevation air is ~25% denser than warm high-elevation air — the cold dome is literally a heavy fluid pooling in the valleys and against the mountains.

P (altimeter) — 75.4% T — 91.2%

2.9 Wind Vector Components u, v Dynamics ▼

Convert meteorological wind direction (from-direction, clockwise from north) to mathematical Cartesian components for vector operations.

u = −|V| × sin(φ_wind) [m/s, positive = eastward] v = −|V| × cos(φ_wind) [m/s, positive = northward] # φ_wind in degrees (meteorological convention: FROM) # Jonas PA mean: ū = −1.43 m/s, v̄ = −2.22 m/s (net WNW flow)

windSpeed — 87.2% windDir — 86.7%

2.10 Kinetic Energy Density KE Dynamics Cold Dome ▼

Mechanical energy in the boundary layer. Near-zero KE inside the cold dome core is the "dead zone" diagnostic for cold air damming — the atmosphere has stagnated.

KE = ½ × ρ × |V|² [J/m³] # Jonas range (PA): 0 – 401 J/m³ # Peak KE: warm sector coastal jet # Near-zero KE: cold dome interior → CAD diagnostic

Jonas range (PA): 0 – 401 J/m³. Mean = 12.3 J/m³. The KE spatial field directly maps the cold dome extent — zero KE in the stagnant cold air, high KE in the overlying warm jet.

ρ (needs P, T) windSpeed — 87.2%

2.11 Wet-Bulb Temperature T_w — The Precip Type Oracle Thermodynamic Precipitation ▼

The temperature reached by evaporating water into the air until saturation. T_w = 0°C is the rain/snow/sleet boundary — not T, not T_d. Mapping T_w across the network is the most direct precipitation-type diagnostic available from surface observations.

T_w ≈ T_C − VPD / 65 [°C] (simplified Normand) # More precise: iterative psychrometric solution # T_w = 0°C: rain/snow boundary # T_w = −2°C: dendritic snowflake growth maximum # Jonas range (PA): −17.8 to +7.8°C

⚠ T_w requires BOTH T and T_d. Only 44% of Jonas PA stations provide this. The T_w = 0°C isoline is the most important contour on the blizzard MRI scan.

T — 91.2% T_d — 58.4% (or RH)

TIER 2

Network-Derived Vector and Gradient Fields

Require spatial gradients across ≥3 non-collinear stations. The dense network becomes an MRI machine here. Mean 6 km PA spacing resolves mesoscale features 15–20 km and larger.

3.1 Spatial Gradient via Least-Squares Plane Fitting Tensor/Field ▼

For any scalar field φ sampled at irregular station locations, the gradient is estimated by fitting a local plane through k nearest neighbors. This is the master operator for all Tier-2 fields.

# Least-squares plane through k neighbors: [∂φ/∂x, ∂φ/∂y] = (A^T A)^{-1} A^T φ # A = design matrix [Δx_i, Δy_i, 1] for each neighbor i # Convert lat/lon to meters: dx = (lon2−lon1) × 111320 × cos(lat_mean) # eastward dy = (lat2−lat1) × 110540 # northward

PA network density check: 1,676 stations in ~55,000 km² → mean spacing ~6 km → Nyquist wavelength ~12 km → resolves features ≥15 km. The Piedmont is well-resolved; ridge tops at 15–30 km spacing need larger kernels.

3.4 Horizontal Wind Divergence δ — Vertical Motion Diagnostic Dynamics Cold Dome ▼

The spreading or contracting of wind at the surface. By mass continuity, surface convergence implies upward vertical motion (clouds, precipitation, cold dome deepening), and surface divergence implies downward motion (subsidence, dome erosion).

δ = ∇_h · V = ∂u/∂x + ∂v/∂y [s⁻¹] # Mass continuity: δ = −∂w_vert/∂z # δ < 0 (convergence) → upward motion → precipitation, deep cold dome # δ > 0 (divergence) → downward motion → dome erosion, clearing # Synoptic scale: ~10⁻⁵ s⁻¹ # Mesoscale bands: ~10⁻⁴ s⁻¹ (resolvable at 6 km spacing)

u, v (wind components) ≥3 nearby stations

3.5 Relative Vorticity ζ — Rotation Signature Dynamics ▼

The rotation of the wind field — cyclonic (CCW) or anticyclonic (CW). The leading edge of the Jonas surface low appeared as a positive vorticity maximum sweeping through PA 12–18 hours before peak winds.

ζ = ∂v/∂x − ∂u/∂y [s⁻¹] # ζ > 0 (cyclonic, NH): upward motion, precipitation # ζ < 0 (anticyclonic): subsidence, cold dome strengthening

u, v (wind components)

3.3 Moisture Flux Vector F_q — Water Vapor Transport Dynamics Precipitation ▼

The horizontal transport of water vapor — the quantity that determines where moisture is delivered to storm systems. Its divergence gives precipitation potential: convergence = moisture piling up = enhanced precipitation.

F_q = ρ × w × V [kg/m²/s] # Divergence = moisture budget: ∇ · F_q = ∂(ρwu)/∂x + ∂(ρwv)/∂y [kg/m³/s] # ∇·F_q < 0: convergence → enhanced precipitation # ∇·F_q > 0: divergence → drying, subsidence

Jonas interpretation: The low-level jet carried a northward moisture flux from the Atlantic. The convergence maximum at the Chesapeake/Susquehanna axis is where snowfall rates peaked (30–60 cm in 24 hours).

ρ (needs P, T) w (needs T_d) u, v (wind)

3.6 Cold Dome Gradient Tensor G_θ — The Diffusion Tensor Image Tensor/Field Cold Dome ▼

The potential temperature gradient tensor encodes the cold dome's strength and geometry. Its eigenvectors reveal the orientation of the cold dome wall and the direction of maximum gradient — analogous to diffusion tensor imaging (DTI) in brain MRI, which reveals white matter fiber tract orientations.

Physical meaning: Just as DTI reveals the brain's fiber tracts, G_θ reveals the cold dome's internal structure, resistance to erosion, and connectivity to Appalachian terrain channeling. The eigenvalue ratio (λ₁/λ₂) measures how "fibrous" the cold dome boundary is vs how isotropic.

θ (needs T, P) ≥4 neighbor stations

3.7 Frontogenesis Function F — Cold Dome Longevity Diagnostic Dynamics Cold Dome ▼

The rate of change of the temperature gradient magnitude — is the cold dome wall sharpening (frontogenesis) or weakening (frontolysis)? This is the single most important diagnostic for cold dome longevity.

D_s = ∂u/∂x − ∂v/∂y # stretching deformation D_h = ∂v/∂x + ∂u/∂y # shearing deformation D = sqrt(D_s² + D_h²) [s⁻¹] F = −½ |∇θ| [(D_s cos 2β + D_h sin 2β) − δ] # β = angle of isentrope relative to x-axis # F > 0: frontogenesis → dome wall intensifying # F < 0: frontolysis → dome wall eroding

Jonas application: Frontogenesis was active along the PA/MD border for 18+ hours, locking the rain/snow line in place and enabling historic snowfall accumulations along the I-95 corridor. F maps the "immune system" of the cold dome — how strongly it resists erosion.

θ gradient u, v gradient divergence δ

3.9 Sensible Heat Flux H — Cold Dome Erosion Rate Dynamics Cold Dome ▼

The turbulent transport of heat from surface to air (or vice versa). During the day when solar heating raises T_s above cold dome air, H > 0 erodes the dome base from below. The rate is calculable from the station network.

H = ρ × c_p × C_H × |V| × (T_s − T_a) [W/m²] # c_p = 1004 J/kg/K # C_H ≈ 0.001–0.003 (bulk transfer coefficient) # T_s from soil sensor or estimated; T_a = air temperature # Cold dome warm rate: ∂θ/∂t ≈ −H / (ρ × c_p × Δz) [K/s] # H=50 W/m², Δz=1000 m → ~4 K warming over 24 hours

T, windSpeed T_soil — 1.4%

3.10 Temperature Advection ADV_T Dynamics Cold Dome ▼

The rate of temperature change due to wind transport. Warm advection (WAA) erodes the cold dome from the south; cold advection (CAA) reinforces it. The classic CAD signature is WAA at the surface being overwhelmed by CAA aloft.

ADV_T = −V · ∇T = −u∂T/∂x − v∂T/∂y [K/s] # WAA (warm advection): ADV_T > 0 → cold dome erosion # CAA (cold advection): ADV_T < 0 → cold dome building

Jonas maximum WAA: Delaware and MD Eastern Shore — first to transition from snow to rain as warm oceanic air undercut the cold dome from the southeast. The ADV_T field revealed this 3–4 hours before T rose above freezing.

T gradient u, v

TIER 3

Spatiotemporal Tensors and Composite Fields

4.1 Atmospheric State Tensor — The MRI Voxel Tensor/Field ▼

After Barnes interpolation to a regular grid, each grid cell holds a state tensor encoding both the field value AND its spatial gradient — exactly what an MRI diffusion tensor voxel encodes. This is the atomic unit of the atmospheric MRI image.

4.2 Barnes Objective Analysis — The Spatial Filter Tensor/Field ▼

The standard method for converting irregular station observations to a regular grid. Two-pass scheme: first pass creates a background field, second pass corrects it. The spatial filter analogous to the MRI slice selection pulse.

# Pass 1 (background): φ^(1)(x,y) = Σ φ_i × w_i / Σ w_i w_i = exp(−r_i² / κ) # Pass 2 (correction): φ^(2) = φ^(1) + Σ(φ_i − φ^(1)_i) × w'_i / Σw'_i w'_i = exp(−r_i² / (γκ)) # r_i = station-to-grid distance # κ = smoothing parameter # γ ≈ 0.3 convergence factor # Jonas PA: κ = (1.5 × 6km)² = 81 km² → resolves features ≥20 km

MetPy implementation: metpy.interpolate.interpolate_to_grid(lon, lat, data, interp_type='barnes', search_radius=40000)

4.3 Wind Shear Tensor J_V — Complete Wind Gradient Tensor/Field Dynamics ▼

The complete 2×2 Jacobian of the horizontal wind. Decomposed into symmetric (strain) and antisymmetric (rotation) parts, it contains ALL dynamically relevant information: divergence, vorticity, and deformation in a single tensor.

J_V = | ∂u/∂x ∂u/∂y | | ∂v/∂x ∂v/∂y | # Trace(J_V) = divergence δ = ∂u/∂x + ∂v/∂y # Antisym component = vorticity ζ/2 # Sym traceless = deformation (D_s, D_h) # Decomposition: J_V = E (strain) + R (rotation) E = ½(J_V + J_V^T) # symmetric R = ½(J_V − J_V^T) # antisymmetric

Power: One tensor, computed at each grid point, encodes divergence + vorticity + deformation simultaneously. Requires 4 spatial derivatives — achievable with ~6–8 stations within a 10 km radius.

4.5 Moisture Balance Budget — The Storm's Water Account Dynamics Precipitation ▼

The complete budget for atmospheric water at the surface level. All three terms are partially observable from surface stations — without a sounding. The residual tells whether the column is gaining or losing precipitable water.

∂W/∂t = −∇ · F_q + E − P # W = precipitable water [mm] # F_q = moisture flux = ρwV # E = evaporation (Penman-Monteith, where solar sensors exist) # P = precipitation rate (precipRate, directly measured) # ∂W/∂t > 0: column moistening → intensifying storm # ∂W/∂t < 0: column drying → storm ending

NEW

New Equations from MADIS Variables Not Yet in Your Compendium

Derived from MADIS variables available in the Jonas dataset that go beyond the original equations.md compendium.

5.1 Pressure Tendency → Surface Development Rate New Dynamics ▼

MADIS provides pressChange3Hour [Pa] at ASOS/AWOS stations. Converting to a rate gives the real-time storm intensification diagnostic.

∂P/∂t = pressChange3Hour / 10800 [Pa/s] # < −4 hPa/3h: explosive cyclogenesis # > +2 hPa/3h: rapid cold dome reinforcement # Map spatially: reveals storm center, intensification rate

5.3 Snowfall Rate from T_w × Q_rate — The Crystal Efficiency Map New Precipitation ▼

Converts the network's precipRate [m/s water equivalent] to snowfall rate using the spatially resolved T_w field. Simultaneously maps crystal type across the storm domain.

z_snow_rate = Q_rate × SWE_ratio(T_w) [m_snow/s]

T_w (°C)	Crystal Type	SWE Ratio	Snowfall Efficiency
−1 to 0	Wet snow / graupel	5–8:1	Low
−3 to −1	Dendritic plates	10–13:1	Moderate
−8 to −3	Plates / columns	8–12:1	Moderate
< −8	Needles / columns	15–25:1	High

Jonas at 12Z: Q_rate reached 2.8×10⁻⁵ m/s at peak (>1 inch/hr liquid), T_w = −2 to −3°C → snowfall rate ≈ 2.5–3.6 cm/hr. Areas with T_w = −5°C saw the same QPF produce 60% more snow depth.

5.4 GPS Precipitable Water — Direct Column Measurement New Thermodynamic ▼

MADIS includes totalColumnPWV [cm] and wetSignalDelay [m] from GPS-MET receivers — a direct measurement of total column water vapor with ~1 mm accuracy. Gold standard calibration for surface proxy PW.

PWV = 10⁴ / (6.3 + T_surface/16.4) × wetSignalDelay [mm] # 14 GPS-MET stations in Jonas full dataset # Accuracy: ~1 mm PWV # Use to calibrate surface proxy: PW_proxy = ρ₀ × w₀ × H_w

GPS-MET stations — sparse

5.5 Bulk Richardson Number Ri_b — Boundary Layer Stability New Cold Dome ▼

Computed from two stations at different elevations — diagnoses whether the cold dome boundary layer is stable, neutral, or convectively unstable.

Ri_b = g × ΔT × Δz / (T̄ × ΔV²) # g = 9.81 m/s², ΔT = T_high − T_low, ΔV = V_high − V_low # Ri_b > 0.25: stable BL → cold dome intact # Ri_b ≈ 0: neutral → vigorous mixing, dome eroding # Ri_b < 0: unstable → rapid cold dome destruction

T, V at two elevation levels

5.7 Penman-Monteith ET — Drought / Water Budget New Dynamics ▼

The physically rigorous evapotranspiration equation. During summer drought periods, stations with solar radiation sensors form a sub-network computing the atmospheric drying rate across the domain.

ET = [Δ(R_n−G) + ρ·c_p·VPD/r_a] / [λ(Δ+γ)] [mm/day] # Δ = de_s/dT [hPa/K] # R_n = net radiation [W/m²] (from solarRadiation) # G = ground heat flux (~10% of R_n, or from soil T sensor) # r_a = aerodynamic resistance ≈ ln(z/z_0)² / (k² × windSpeed) # λ = 2.45 MJ/kg (latent heat evaporation) # γ = 0.067 kPa/°C (psychrometric constant)

solarRadiation — 13.9% T, VPD, windSpeed soilTemperature — 1.4%

5.2 Precipitation Peakedness — Stratiform vs Convective Texture New Precipitation ▼

Uses MADIS multi-period accumulation fields (precip1min, precip3hr, etc.) to diagnose whether precipitation is steady/stratiform or convective/banded — without a radar.

Peakedness = Ī_1min / Ī_3hr_mean # Peakedness ≈ 1: steady stratiform (lake-effect, frontal) # Peakedness > 3: convective banding (thunder-snow, heavy convection) # Note: precip1min sparse in Jonas PA subset

PROTOCOL

The 5-Layer Atmospheric MRI Scan

Applied to Blizzard Jonas and all major storm types. Each layer is a separate diagnostic panel — read in sequence like an MRI stack.

Layer 1 · Mass Field
Cold Dome Geometry
The structural anatomy of the cold air mass — its boundaries, depth, and density distribution. The "T1-weighted image" of the storm.

          θ field
          ρ field
          Barnes P
          G_θ tensor
        
Layer 2 · Moisture Field
Snow/Rain Boundary
The T_w = 0°C isoline is the rain/snow line. Maps crystal type and SWE ratio across the domain. The "contrast-enhanced" panel.

          T_w isoline
          w gradient
          VPD field
          SWE ratio
        
Layer 3 · Wind Field
Transport & Deformation
The circulatory system — where air is moving, rotating, and deforming. Equivalent to blood flow imaging in the body.

          (u,v) vectors
          ∇·V divergence
          ζ vorticity
          J_V tensor
        
Layer 4 · Flux Field
Convergence Zones
Where moisture is accumulating (precipitation enhancement) and where temperature is changing (dome erosion or reinforcement).

          ∇·F_q
          ADV_T
          ADV_q
          H flux
        
Layer 5 · Instability Field
Banding Potential
Thunder-snow and convective snowfall bands — the "hot spots" in the scan. Detectable only with the dense network resolving 15–20 km features.

          θ_e gradient
          CAPE proxy
          Frontogenesis F
          KE field
        

Cold Dome Erosion Mechanisms

1. Bottom Erosion (H > 0)

Solar heating raises T_s above cold dome air → H > 0 → warm rate ≈ 4 K/24h. Visible as warm spots in valley stations while ridge stations remain cold. Detectable at 1–2 hr resolution.

2. Side Erosion (WAA)

Warm advection from SW → gradual westward retreat of θ contours. ADV_T field shows the erosion front in real time. Most common path during spring and fall CAD events.

3. Top Erosion (Subsidence)

Diagnosed at surface by pressure rise + T rise WITHOUT warm advection. The cold dome is being compressed from above by descending air. Invisible to a sparse network — needs 6 km density.

4. Precipitation Erosion

Evaporative cooling maintains dome depth even as other mechanisms warm it. The precipRate × VPD term — measured directly by the network — counteracts erosion and can lock the cold dome in place for 24–48 hours.

IMPL

Variable Priority Tree & Pipeline

Given partial data availability, compute derived fields in this priority order. Each level unlocks the next.

Level 0

T, z, φ, λ — always available

Temperature · Elevation · Lat/Lon

91.2%

Level 1

θ = T(1000/P)^0.286 · e_s · ρ

Needs P_alt (altimeter)

75.4%

Level 2

u, v wind components

Needs windSpeed + windDir

87.2%

Level 3

KE = ½ρ|V|² · PGF

Needs Levels 1+2

85%

Level 4

e · w · T_v · T_w · VPD · θ_e

Needs T_d (dew point) or RH

58.4%

Level 5

∇T · ∇·V · ζ · ADV_T · ∇P

Needs ≥3 neighboring stations at each level

network

Level 6

F_q · ∇·F_q · ADV_q · Frontogenesis F

Needs Level 4 moisture + Level 5 gradients

~40%

Level 7

G_θ tensor · J_V tensor · State tensor S(x,y,t)

Full network coverage + Barnes analysis

gridded

# MADIS Atmospheric MRI Pipeline — Python skeleton
import numpy as np
import pandas as pd
from metpy.interpolate import interpolate_to_grid

# 1. Load MADIS subset
df = pd.read_parquet('madis_20160123_1200_pa_midatlantic_subset.parquet')

# 2. Tier-1: single-station derived fields
T_C = df['temperature'] - 273.15
P   = df['altimeter'] / 100  # Pa → hPa

df['e_s']   = 6.112 * np.exp(17.67 * T_C / (T_C + 243.5))
df['theta'] = (df['temperature']) * (1000 / P)**0.286
df['rho']   = (P * 100) / (287.05 * df['temperature'])
df['u']     = -df['windSpeed'] * np.sin(np.radians(df['windDir']))
df['v']     = -df['windSpeed'] * np.cos(np.radians(df['windDir']))
df['KE']    = 0.5 * df['rho'] * df['windSpeed']**2

# 3. Moisture fields (where T_d available)
mask = df['dewpoint'].notna()
Td_C = df.loc[mask, 'dewpoint'] - 273.15
df.loc[mask, 'e']    = 6.112 * np.exp(17.67 * Td_C / (Td_C + 243.5))
df.loc[mask, 'w']    = 0.622 * df.loc[mask,'e'] / (P[mask] - df.loc[mask,'e'])
df.loc[mask, 'T_v']  = df.loc[mask,'temperature'] * (1 + 0.61*df.loc[mask,'w'])
df.loc[mask, 'VPD']  = df.loc[mask,'e_s'] - df.loc[mask,'e']
df.loc[mask, 'T_w']  = T_C[mask] - df.loc[mask,'VPD'] / 65
df.loc[mask, 'theta_e'] = df.loc[mask,'theta'] * np.exp(
    2.5e6 * df.loc[mask,'w'] / (1004 * df.loc[mask,'temperature']))

# 4. Barnes interpolation to 5km grid
lon, lat = df['longitude'].values, df['latitude'].values
grid_theta = interpolate_to_grid(lon, lat, df['theta'].fillna(method='nearest'),
    interp_type='barnes', hres=5000, search_radius=40000)

REF

Quick Reference — All Equations

Equation	Formula	Units	Requires	Availability
e_s	6.112 exp(17.67 T_C / (T_C+243.5))	hPa	T	91%
e	6.112 exp(17.67 T_d / (T_d+243.5))	hPa	T_d	58%
w	0.622 e / (P−e)	kg/kg	T_d, P	58%
T_v	T (1 + 0.61w)	K	T, w	58%
θ	T (1000/P)^0.286	K	T, P	75%
θ_e	θ exp(Lw/c_pT)	K	T, P, w	58%
ρ	P / (R_d T_v)	kg/m³	P, T_v	75%
u	−\|V\| sin(φ_wind)	m/s	speed, dir	87%
v	−\|V\| cos(φ_wind)	m/s	speed, dir	87%
KE	½ρ\|V\|²	J/m³	ρ, V	85%
T_w	T_C − VPD/65	°C	T, T_d	58%
VPD	e_s − e	hPa	T, T_d	58%
F_q	ρ w V	kg/m²/s	ρ, w, V	~40%
δ (divergence)	∂u/∂x + ∂v/∂y	s⁻¹	u, v, spacing	network
ζ (vorticity)	∂v/∂x − ∂u/∂y	s⁻¹	u, v, spacing	network
ADV_T	−u ∂T/∂x − v ∂T/∂y	K/s	T, V, spacing	network
H (heat flux)	ρ c_p C_H \|V\| (T_s−T_a)	W/m²	T, V, T_soil	1.4%
Frontogenesis F	−½\|∇θ\|(D_s cos2β + D_h sin2β − δ)	K/m/s	θ, V, spacing	network
Snow rate	Q_rate × SWE(T_w)	m/s	Q_rate, T_w	~50%
GPS PWV	10⁴/(6.3+T/16.4) × wetSignalDelay	mm	GPS-MET	sparse
Ri_b	g ΔT Δz / (T̄ ΔV²)	—	T, V at 2 levels	pairs
ET (P-M)	[Δ(R_n−G) + ρc_p VPD/r_a] / [λ(Δ+γ)]	mm/day	T, RH, I_solar, V	14%
∂P/∂t	pressChange3hr / 10800	Pa/s	ASOS press	ASOS only
G_θ tensor	[[∂θ/∂x, ∂θ/∂y],[∂θ/∂y, −∂θ/∂x]]	K/m	θ, ≥4 stations	network
J_V tensor	[[∂u/∂x, ∂u/∂y],[∂v/∂x, ∂v/∂y]]	s⁻¹	u, v, ≥6 stations	network

Availability color key: ■ ≥80% robust ■ 60–80% good ■ 40–60% moderate ■ <15% sparse ■ network requires dense coverage

The Dense Mesonet as aCAD-Aware Atmospheric MRI