This page describes the variational cloud / hydrometeor analysis and short-range nowcast running in real time in the LAPS 500 m Front Range domain, and sets out what it adds beyond a simple two-dimensional advection nowcast, together with its current limitations. Every number quoted below is measured on case studies, and each is named so the claims can be re-checked.
The analysis is a 4-D variational (4D-Var) fit of a cloud and hydrometeor field to satellite and radar observations over a short assimilation window. The window is 5 minutes and the analysis time is the window end, so the product is valid at the most recent observation rather than 5 minutes behind it.
The control variable is the hydrometeor state at the window start, carried per species on the full 3-D grid: cloud liquid water, cloud ice, rain and snow. Graupel is carried through the model as a passive tracer and contributes to the reflectivity operator, while remaining outside the control vector.
Writing x0 for that initial state and xb for the background, the cost function assembles four kinds of term:
J(x0) = J_b(x0) background departure
+ SUM over obs times t of
w(t) * J_ir(t) IR brightness temperature
+ w_vis(t) * J_vis(t) visible reflectance factor
+ J_z(t) radar reflectivity, in Z-space
+ J_rh(x) humidity consistency
+ J_ph(x0) precipitation phase
Two distinct maps appear in the cost function, and it helps to keep them apart. The forward model propagates the state forward in time, from the window start to each observation time. The observation operator of section 2a then maps that state into observation space, rendering condensate as brightness temperature, reflectance factor and reflectivity. The adjoint runs back through both. The distinction is worth keeping: a model prognosticates, an operator evaluates. (In data-assimilation notation these are M and H respectively. Usage varies by field: the observation operator is sometimes termed a forward operator, particularly for radar, and remote-sensing work often calls it a forward model. On this page forward model always means the time propagation.)
Between observation times the state is carried by semi-Lagrangian advection of each species on the effective wind field described in section 3a, sub-stepped to respect the Courant limit. The trajectory weights are computed once and shared across species, so cloud and precipitation in the same air move together. The default forecast beyond the analysis time uses this same advection with no source or sink terms; section 2 sets out the prognostic processes that exist, and their status.
The forward model, the radiative-transfer operators and the reflectivity operator are all written in MLX and differentiated by automatic differentiation on the Metal GPU. The adjoint is therefore exact with respect to the code that runs, so no separately maintained adjoint can drift out of step. The practical cost of this appears in section 2b: every process has to stay differentiable, ruling out hard threshold gates.
Minimization is by nonlinear conjugate gradient over roughly ten iterations, arranged in two blocks: a cloud block over liquid water and ice, and a precipitation block over rain and snow. Each block runs its own line search, because their gradient magnitudes differ by many orders of magnitude and a shared step length serves neither. Step lengths are warm-started from the previous cycle, cutting the first iteration substantially, and a failed line search puts its block into a brief cooldown rather than stopping the minimization.
The useful distinction here is diagnostic versus prognostic. Diagnostic physics takes the state as given and determines how it appears: the microphysical assumptions inside the observation operator, turning condensate into radiances and reflectivity. Prognostic physics evolves the state, creating, converting, moving and removing condensate over time. Nearly all of the physics active every cycle is diagnostic; nearly all of the prognostic physics either sits behind a switch or awaits implementation. That split explains a good deal of the system's behaviour, and the three subsections below follow it.
The phase reference is worth setting out in detail, since it governs both the humidity term of section 1a and the saturation quantities the operators depend on. Referenced to liquid at every temperature, a flat 70% relative-humidity threshold grows progressively more severe with height:
| Level | Temperature | 70% over liquid equals |
|---|---|---|
| 100 hPa | −71.7 °C | 133% over ice |
| 150 hPa | −63.3 °C | 126% over ice |
| 200 hPa | −50.1 °C | 113% over ice |
| 250 hPa | −38.5 °C | 102% over ice |
| 300 hPa | −29.5 °C | 93% over ice |
A liquid-referenced penalty therefore demands ice supersaturation at every cirrus level, and does so most strongly where it is least defensible, at the coldest levels where the weighting factor also grows. A second error compounds it: saturation specific humidity computed over liquid runs too large below freezing, so the condensate-to-saturation ratio that scales the penalty runs too small. Such a term is at once too strict in its threshold and too weak in its magnitude, leaving the weight parameter able only to trade one error against the other.
The correction blends saturation vapour pressure between the liquid and ice formulations using the same ice water fraction above, with the ice branch taken from the LAPS Goff–Gratch esice routine so that the analysis and the underlying cloud analysis agree by construction. Blending across a ramp rather than switching at 0 °C keeps the cost smooth, and that matters because a kink in this same cost had previously stalled the line search. The effective liquid-relative threshold then varies with temperature: about 70% near freezing, 59% at −40 °C, and 37% at −72 °C. Section 6 reports its measured effect.
All of the following are implemented as elementwise operations so the adjoint flows through them, and each has its own switch. All of them stay off in the deployed configuration.
| Process | Form | Note |
|---|---|---|
| Sedimentation | 1-D upwind fall in z, sub-stepped for stability; rain 5 m/s, snow 1 m/s; mass past the lowest level leaves as surface precipitation | The module's own documentation calls this required for radar: rain at 5 m/s falls about six model levels in five minutes, so omitting it is a first-order error over the window |
| Vertical advection | Hydrometeors carried by the diagnosed vertical velocity | Two diagnosed w fields exist, orographic (terrain-following, decay scale 2500 m) and thermal/buoyant, while the transport itself stays off |
| Cloud condensation / evaporation | Relaxation toward a saturation- and updraft-implied equilibrium, with an additional one-sided evaporation sink in subsaturated air | This is the sink that would erode a spreading anvil edge |
| Adiabatic lifting source | dqc/dt = −w ∂qsat/∂z for saturated ascent, applied as a direct source so it scales linearly with w | Replaces an earlier gain-and-cap form that saturated for any w above 0.0025 m/s, roughly 70 times below real terrain-scale ascent, so every rising point saw the same equilibrium |
| Subsidence evaporation | Decay rate on existing condensate proportional to descent speed | The lee half of an orographic couplet. Without it the scheme adds cloud on a windward slope while leaving it downstream, so terrain cloud smears instead of staying anchored |
| Autoconversion | Cloud → precipitation above a content threshold, separate thresholds for liquid and ice | — |
| Accretion | Collection of cloud water by falling precipitation, Kessler-scaled | The first-cut rate was 40× too large, so the per-step removal fraction saturated and every cell containing rain lost all its cloud water in one step. At a defensible rate it is statistically indistinguishable from being off |
| Precipitation evaporation / sublimation | Sink in subsaturated air, returning mass to vapour | The other missing leading-edge sink, on the precipitation side |
The intended path is that the nowcasting machinery supersedes itself: the steering blend reduces to the real winds when its weight is set to zero, it is applied upstream of the advection routine so that improved dynamics modify the same wind field, and the steering vector is a placeholder for a quantity the dynamics should eventually produce rather than be told.
A single-layer 2-D advection nowcast moves everything at one speed. The atmosphere behaves otherwise. Measured on the case of 25 July 2026, 0050 UTC, the precipitating system and the cloud field move at speeds that differ by a factor of three:
| Field | Best-fit motion | Speed |
|---|---|---|
| Echo > 20 dBZ (15-min baseline) | +3.0 km E, −1.5 km N | 3.7 m/s |
| Convective cores > 35 dBZ | +2.5 km E, −1.5 km N | 3.2 m/s |
| IR cloud (anvil / cloud top) | +3.0 km E, −1.0 km N | 10.5 m/s |
Both are consistent with the wind profile: 250 hPa u = 15.2 m/s, 500–550 hPa 4.4–4.7 m/s, 700 hPa 0.7 m/s. The IR sees the anvil drifting with the upper winds while the precipitating system is steered by the mid levels.
Purely per-level 3-D advection fails for the opposite reason: under 19.6 m/s of shear it separates the echo levels by about 20 km in 15 minutes, tearing a storm apart vertically because nothing in a kinematic model holds it together. So the effective wind is a blend,
u_eff = u_3d + w(precip content) * (u_steer - u_3d)
where u_steer is a single deep-layer steering vector fitted from consecutive radar mosaics, and w ramps smoothly with local precipitation content. The result is 2-D coherent motion where a deep precipitating core should move as a unit, and the real 3-D winds in the anvil and in non-precipitating cloud, with a smooth transition so no artificial convergence is created at the storm edge. The blend engages on only about 2% of grid points, a deliberately narrow intervention.
Measured effect, 3-D reflectivity Equitable Threat Score at 20 dBZ, steered versus unsteered on the case of 25 July 2026, 0050 UTC:
| Lead | Unsteered | Steered | Change |
|---|---|---|---|
| +0 min (window end) | 0.6160 | 0.6160 | identical — clean control |
| +5 min | 0.4897 | 0.5289 | +8.0% |
| +10 min | 0.4010 | 0.4474 | +11.6% |
| +15 min | 0.3412 | 0.3917 | +14.8% |
The gain grows with lead, exactly as an accumulated displacement error should, and core-level (300–700 hPa) displacement spread tightens from 10.8 km to 3.8 km. On 28 July 2026 the vector was fitted from 22,299 echo points and reduced the mosaic-to-mosaic rms from 19.67 to 15.31 dBZ, a 22.2% improvement, before any forecast was run.
A 2-D advection nowcast can only move what it is given. The variational analysis can add what is missing. On 28 July 2026, 1800 UTC, the background carried far too little light echo; the increment supplied it, and the forecast then tracked the observations almost exactly. Column-maximum reflectivity, pixel counts:
| Valid | Observed ≥5 dBZ | 4D-Var | Model only (no increment) | Increment |
|---|---|---|---|---|
| 1755 (window start) | 56,419 | 51,208 | 33,930 | +17,278 |
| 1800 (analysis) | 60,459 | 38,825 | 29,755 | +9,070 |
| 1805 | 65,540 | 49,608 | 36,842 | +12,766 |
| 1810 | 69,011 | 61,187 | 45,789 | +15,398 |
| 1815 | 71,533 | 72,123 | 55,537 | +16,586 |
The "model only" column is the identical configuration run with the minimizer disabled, so the difference is the analysis increment and nothing else. It closes a genuine 22,000-pixel deficit against the observations, and the +20 minute forecast lands at 72,123 against 71,533 observed. That is information the observations carry, and assimilation is what delivers it.
The strictest control is to run the same forecast model from the LAPS cloud analysis itself, so that only the initial condition differs. On 25 July 2026, 1800 UTC, using a combined figure of merit (visible RMSE × 100 + IR RMSE):
| Lead | CA-initialized | 4D-Var | Change |
|---|---|---|---|
| +0 (analysis) | 7.72 | 6.74 | −13% — 4D-Var better |
| +5 min | 11.50 | 11.30 | −1.7% — 4D-Var better |
| +10 min | 14.67 | 14.90 | +1.6% — CA better |
| +15 min | 17.38 | 17.79 | +2.4% — CA better |
Section 2b lists the prognostic processes built into the model and held off. The reasons differ by process, and only one has so far been shown to degrade the forecast.
An important limitation applies to most of these tests: the analysis-time fit is a blunt instrument for them. A physically sound source or sink is absorbed by the minimizer: the analysis-time IR fit stayed flat at 6.84 K across a twenty-fold range of subsidence-evaporation rate, up to 97% decay per step, because the free condensate control compensated for the change the physics introduced. Only a pathological setting shows up at analysis time, and then only by breaking the fit outright. So a null result at the analysis time carries little diagnostic weight for these processes. They must be judged on the forecast, where the minimizer has finished and can absorb nothing. Several of the flat readings above were collected with the absorbing instrument, and the same redundancy is a large part of why the analysis and the forecast keep disagreeing about what helps.
The status of section 2b is therefore best described as awaiting full evaluation rather than as unpromising, and for sedimentation and the lifting source the existing evidence favours enabling them. The one caution is ordering: as section 5 argues, a non-conservative advection scheme is currently manufacturing condensate, so sinks introduced before that is fixed would be tuned against a numerical error rather than against the atmosphere.
One structural constraint shapes all of these: every process must be differentiable, because the adjoint is taken through the code that runs. Hard threshold gates have zero gradient below the knee, silently killing the adjoint wherever they bite, so the saturation and humidity gates are smooth ramps rather than clips, with the smoothing width set to zero recovering the original hard behaviour bit-for-bit. This is a real design cost of automatic differentiation, and it is why the rate constants above need re-derivation rather than direct adoption from a conventional microphysics scheme.
The forecast visibly grows a spreading anvil shield, and it grows it too fast: on 28 July 2026 the observed ≥5 dBZ area grew 27% in 20 minutes while the model grew 41%. The natural physical explanation is that an anvil spreading with no evaporation of cloud and no sedimentation of precipitation retains all its condensate at the leading edge. That explanation is correct as far as it goes, and those sinks really are absent from the free forecast.
But a stronger and more fixable effect is also present, and it would be masked rather than cured by adding sinks. Running with the minimizer disabled, so the forecast is pure advection of a conserved scalar, total condensate still grows monotonically:
| Valid | Total condensate | Change per 5 min |
|---|---|---|
| 1755 | 249.22 | — |
| 1800 | 254.54 | +2.13% |
| 1805 | 261.46 | +2.72% |
| 1810 | 271.53 | +3.85% |
| 1815 | 282.89 | +4.18% |
Advection of a conserved scalar conserves its total mass, so the growth is numerical in origin: the semi-Lagrangian scheme uses a non-monotone cubic (Catmull-Rom) interpolation that overshoots at sharp gradients. Two diagnostics localize it:
So the overshoot has both a missing-sink component and a scheme component, and the order of work matters: a monotone or flux-conserving advection scheme should come before, or at least alongside, evaporation and sedimentation. Adding sinks first would tune them to cancel a numerical error, and would then be mis-tuned as soon as the scheme was fixed. A flux-form conservative scheme conserves exactly, though it costs roughly ten times as much, so the practical question is whether a monotonicity limiter on the existing interpolation recovers most of the benefit.
Section 2a sets out the correction itself. Its measurement is reported here in full, because the outcome differed from the expectation twice over: in the channel it acted through, and in its size. It also accounts for an earlier humidity-weight sweep where moderate values did nothing and a large value collapsed the upper-level cloud while degrading the IR fit by 7%, behaviour consistent with a mis-specified penalty applied hard.
The expectation was that this would work through the humidity penalty. The larger effect arrived instead through the observation operator, because the same relative humidity feeds the convective cloud-fraction diagnosis. Ice-aware humidity reads higher in cold air, so diagnosed cloud fraction rose:
| Domain | Mean diagnosed cloud fraction | Cells below 0.8 |
|---|---|---|
| Full domain | 0.82 → 0.94 | 44% → 6% |
| Cloudiest 200×200 subset | 0.60 → 0.87 | 100% → 20% |
This is a change in how a state renders rather than in the state itself, and the fields confirm it. Inside the analysed region the two configurations agree on condensate to within 0.3% (7.254 against 7.234, with 125,629 against 125,430 cloudy points, difference RMS 3.1×10−8) while their visible fits differ by 24%.
Using the combined figure of merit, ice-aware against liquid-only:
| Domain | Background | Analysis | Window start |
|---|---|---|---|
| Cloudiest 200×200 subset | −31.3% | −16.5% | −16.5% |
| Full domain | −7.6% | +0.9% | −3.8% |
The subset promised a 16% improvement and the full domain returned about 1% the other way. Two effects account for the reversal, and both generalise beyond this correction:
Across the forecast leads on the full domain the ice-aware configuration runs about 1% behind at +5, +10, +15 and +20 minutes, with IR marginally ahead and visible marginally behind at every lead. On a separate daytime convective case the same isolated comparison came out roughly 0.7% ahead at analysis time. Two cases, comparable magnitude, opposite signs: the honest reading is that the combined score is indifferent to this change.
It stays in operational use for two reasons that survive that verdict. It is the physically correct formulation, and it removes a penalty that charged legitimate cirrus for existing. And it clearly improves the window-start fit, where the IR error moves from a 4% degradation to a 7% improvement, the end of the window that matters for getting the cloud evolution right.
One prediction from this work remains open. The corrected threshold should retain condensate near 125–175 hPa, where the air sits close to ice saturation, while thinning the genuinely dry levels at 100 hPa and 225–300 hPa. The stratiform test case produced almost no upper-level veil to begin with, 4 points against 1 above 200 hPa, so that test awaits a case carrying real cirrus. The 28 July 2026 case has 4,828 such points and is the intended vehicle.
Organised by component, following sections 1 to 3, so that each entry names the part of the algorithm it belongs to. Only the first is a strict prerequisite for the others.
This comes first, since section 5 shows the semi-Lagrangian scheme manufacturing condensate: a sink enabled before the scheme improves would be tuned against a numerical error, and mis-tuned as soon as the scheme was fixed. The cheap candidate is a monotonicity limiter on the existing cubic interpolation, clipping it to the local range of the donor stencil; a flux-form conservative scheme conserves exactly at roughly ten times the cost. The readout is simple: total condensate under a pure-advection run should stay flat.
Once the numerics allow it, the processes of section 2b are worth enabling in order of the evidence behind them:
Beyond those, graupel generation and fall from section 2c are the leading candidate for the complete absence of core intensification described in section 4, since riming is what builds a convective core's reflectivity. Graupel currently advects as a passive tracer and appears in the reflectivity operator, with creation, removal and fall all still to be added. This is the single largest identified gap between the forecast and the observations.
Every change measured so far trades IR against visible. The visible residual is dominated by placement rather than amplitude, so height-aware cloud placement is the structural candidate for improving IR while keeping the visible gains. The warm IR bias drift with lead, present in both case studies, is the symptom to track. The phase reference corrected in section 2a also carries an untested prediction here: it should retain cirrus near 125–175 hPa while thinning the genuinely dry levels, and the 28 July 2026 case carries 4,828 suitable points to test it on.
The precipitation block intermittently accepts no step, leaving the radar term inert for that cycle. It is now instrumented with a gradient-alignment diagnostic, and the reported descent direction stays sound while the line search still fails, so the cost and its gradient disagree somewhere in the precipitation path. Diagnosis awaits a case carrying active precipitation, since a precipitation-free case leaves the block nothing to fit and the same symptom means nothing there.
Two measurement rules follow from what section 4a establishes. Physics is judged on the forecast rather than the analysis fit, since the minimizer absorbs a physically sound source or sink and an analysis-time null says little about it. And the unassimilated-cloud-analysis control of section 3c is re-run after each substantive change, as the cleanest single measure of whether the assimilation earns its cost.
This remains a nowcast throughout. The target range stays at zero to a few hours, where observed structure still carries most of the skill, and adding dynamics leaves that unchanged. What changes is how much of the forecast advection has to carry on its own. Today it carries almost all of it, with a fitted steering vector standing in for storm motion.
The steering blend of section 3a is deliberately built to give way. Its weight reduces to zero to recover the real 3-D winds exactly, it sits upstream of the advection routine so improved dynamics modify the same wind field, and the fitted vector stands in for a quantity the dynamics should eventually produce on their own. As vertical transport, latent heating feedback and a mass-consistent wind field arrive, advection becomes one term among several rather than the whole forecast, and the prescribed storm motion can relax toward zero weight as the dynamics take the work over.
All figures are single-case measurements unless stated otherwise. The combined figure of merit weights visible error roughly 2:1 against IR in practice, and a different weighting would move the crossover lead in section 3c.