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.
Figure 1.1. The assimilation window: observations at each time are compared against the model state advected from the single set of control variables at the window start, and the adjoint returns the gradient of the total cost with respect to those controls.
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 |
Table 2.1. Relative humidity over liquid that corresponds to ice saturation, by level. The 70% liquid threshold the older scheme used reaches 133% over ice at 100 hPa, which is why a single liquid-referenced threshold cannot serve the cold levels.
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.
Height inversion from observed cloud drift. Deployed 1 September 2026. Infrared radiances constrain cloud-top temperature, and therefore height, far better than they constrain cloud amount, but the two trade off against each other: a cloud can be made to match an observed brightness temperature either by moving it or by thickening it, and the minimization has no strong preference. Observed cloud motion breaks that degeneracy from outside the radiance, because the speed and direction at which a cloud drifts identify the wind level it sits at.
The analysis therefore measures cloud drift by cross-correlation between successive visible images, compares it against the drift implied by the analysed cloud height, and where the two disagree redistributes condensate toward the level whose wind matches the observation. On the case it was developed for, this reduced the disagreement between observed and analysed cloud motion from 12.2 m/s to 1.5 m/s.
The redistribution is applied only where a stack of gates agrees that it is warranted, each gate stating a physical condition rather than a tuned threshold: the cloud must be optically thin, actually advecting, measurably misplaced, and the required motion must be achievable at some wind level in the column. All are smooth ramps rather than steps, because hard gates left visible discontinuities at their boundaries.
A terrain ceiling carries most of the benefit, and it is the one gate whose case is physical rather than empirical: cloud over high terrain is anchored orographic cloud, not advecting cirrus, and lifting it to a fast-wind level makes it advect away when it should stay put. Excluding columns above 2800 m improved the analysis figure of merit, preserved the motion gain, and halved the tail of forecast damage. Three earlier attempts to gate on cloud properties instead — observed brightness temperature, and target-level ice phase under two configurations — all failed the same way: with the other gates already keeping under 8% of columns, any further multiplicative ramp acts as a volume control rather than a discriminator. That is recorded here because it is a negative result worth not repeating.
The caveat is that terrain and regime are confounded in the single case this was established on, since the high terrain is also where that case's orographic cloud sat. A case with genuine mountain cirrus would separate "do not lift anchored orographic cloud" from "do not lift over mountains"; the second reading would suppress real cirrus over the Rockies. No such case is archived yet, and acquiring one is the highest-value outstanding item.
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 | Switch | Form | Note |
|---|---|---|---|
| Sedimentation | --precip-sink (on by default in the physics module; autoconversion is what feeds it) | 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 | --vertical-advect needs --fcst-conserve for the forecast path | 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 | --window-condense (in the window) --fcst-condense (in the forecast) | 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 | --lift-source, --k-lift | 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 | --oro-couplet, --oro-k-subs | 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 | --precip-sink, --k-auto | Cloud → precipitation above a content threshold, separate thresholds for liquid and ice | — |
| Accretion | --accrete, --k-accr | 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 | no switch module parameter do_pevap, reachable only in code | Sink in subsaturated air, returning mass to vapour | The other missing leading-edge sink, on the precipitation side |
Table 2.2. Prognostic processes that are implemented and differentiable but held off by default, with the switch that enables each and the form it takes. They are held off pending the conservative advection scheme of section 7a, not because they are unfinished.
Two cautions that the switch names do not convey. Condensation has two switches covering different intervals — one inside the assimilation window, one in the forecast that follows it — and enabling the wrong one leaves the interval of interest as pure advection. And the humidity-control prototype runs a reduced cost function: it omits the cloud-presence and motion terms, so switches feeding those are inert whenever it is active, however well they work on their own. A switch is only live if the term it feeds is in the cost function actually being minimized.
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 |
Table 3.1. Best-fit motion of the reflectivity field against the cloud field over a 15-minute baseline. The two disagree in both speed and direction, which is what the hybrid advection of this section exists to represent.
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% |
Table 3.2. Forecast skill with and without the steering-wind blend, by lead. The window end is identical by construction and serves as the control that the comparison is clean.
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 |
Table 3.3. Echo coverage above 5 dBZ at each window time: observed, analysed, and the model run without the increment. The difference between the last two is the information the assimilation supplies.
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 control is to run the same forecast model from the unassimilated cloud field, so that only the initial condition differs. That comparison run is called bgfcst, the background forecast. It is seeded from the LAPS cloud analysis at the analysis time and then shares everything else with the 4D-Var forecast: the same steering-blended winds, the same five-minute stepping, and the same physics path. Measured on 25 July 2026, 1800 UTC, over a domain-scale box, using a combined figure of merit (visible RMSE × 100 + IR RMSE), with both columns scored against the same observations at the same valid times:
| Lead | Valid | bgfcst | 4D-Var | Change |
|---|---|---|---|---|
| +0 (analysis) | 1800 | 7.55 | 6.56 | −13.1% — 4D-Var better |
| +5 min | 1805 | 12.93 | 9.65 | −25.4% — 4D-Var better |
| +10 min | 1810 | 15.17 | 12.93 | −14.8% — 4D-Var better |
| +15 min | 1815 | 17.35 | 15.59 | −10.1% — 4D-Var better |
Table 3.4. Forecast error against a run started from the unassimilated first guess, by lead. Negative change favours the assimilated start.
An earlier version of Table 3.4 showed the advantage crossing over near +10 minutes, with the comparison run ahead at +10 and +15. That version was not a controlled comparison. The two forecasts ran through different code paths and differed in four ways rather than one: the comparison run used raw LGA winds against the 4D-Var run's steering-blended winds, took a single call of the full lead against five-minute frames, ran no condensation physics, and was seeded from the back-advected window-start background while carrying a reftime stamp for the analysis time — so each of its leads was labelled about one window later than it was valid. All four differences favoured 4D-Var, so the crossover was reported despite them rather than because of them. They were removed on 2 September 2026.
Table 3.4 is therefore a new measurement rather than a corrected version of the old one: the configuration also changed over the intervening month, and the scoring box is not guaranteed to match. The analysis row is the one quantity a forecast-side fix cannot touch, and it reproduces — −13.1% against the earlier −13% — which is what makes the two tables comparable enough to set side by side at all. The seeding switch that makes the reftime stamp truthful is off by default and has to be set explicitly; any future run of this comparison must set it.
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% |
Table 5.1. Total condensate through the forecast, showing the semi-Lagrangian scheme manufacturing mass. The per-interval change is positive under pure advection, where it should be zero.
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% |
Table 6.1. Diagnosed cloud fraction before and after the ice-phase saturation correction, with the count of cells below 0.8. The correction acts through the cloud-fraction channel rather than through the condensate directly.
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% |
Table 6.2. Size of the ice-phase correction on the full domain and on the cloudiest subset. The effect is domain-dependent, so the subset figure is not a general magnitude.
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 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 ice-subsaturated 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 ice-subsaturated 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.
A case observed on 2 September 2026 at 1400 UTC turns the argument of 7f into a measurement. Mid-level cloud lay across the mountains, and the structure has two parts moving at different speeds: the wave crests are locked to the terrain while the cloud envelope advects through them. The forecast described in section 3a carries everything at one velocity, so it reproduces the envelope and carries the crests away with it.
Writing it as an equation makes the missing term explicit. What the observation requires is
dq/dt + U . grad(q) = S(x)
with S a source and sink fixed in space by the terrain, and its amplitude set by the
moisture the envelope delivers. Our forecast is the same expression with S = 0. The
two degenerate alternatives both fail: pure advection moves the crests, and holding the
field still would fix the crests but destroy the envelope, which we currently get right.
Scoring the forecast against persistence, on observed-cloudy pixels, separated into terrain relief terciles. Positive means the 4D-Var forecast is ahead; the table gives infrared RMS margin in kelvin against forecast lead:
| Case | Zone | analysis | +5 | +10 | +15 |
|---|---|---|---|---|---|
| 2 Sep, wave | high relief | +0.17 | −0.76 | −1.00 | −1.44 |
| flat | +2.44 | +2.49 | +3.38 | +3.19 | |
| 28 Jul, convective | high relief | +3.45 | — | — | +7.30 |
| 26 Jul, stratiform | high relief | +0.61 | — | — | +0.31 |
Table 7.1. Translation gain by terrain zone and lead for the two cases. Gain is the correlation at the best along-wind shift minus the correlation at zero shift; near-zero means the feature stands still.
Three things follow. On the wave case the forecast loses to persistence over high terrain and the gap widens with lead, while over flat ground in the same domain, at the same leads, it wins and its margin grows. The flat zone is the control: the model, the winds and the scheme are identical there, so this is not a general forecast failure but one tied to terrain. And the effect is not generic to mountains — on an ordinary convective case the forecast beats persistence over high relief by a growing margin, the opposite sign. The stratiform case sits between them. What distinguishes the wave case is that the terrain-locked part of the pattern carries enough of the variance to dominate the error.
A second signature points at the missing sink specifically. Splitting a ring around a Sawatch summit into upwind and lee halves by the observed flow direction, the forecast at +15 retains a 1.87 K advantage over persistence upwind and is level with it in the lee. Losing skill preferentially downstream is what a missing evaporation sink predicts, and no wind error produces that asymmetry. The samples are small — a few hundred pixels per half — so this is supporting evidence rather than a result on its own.
Cloud motion measured by cross-correlation returns 0.7 m/s in the visible and 2.2 m/s in the infrared where the analysed field moves at 8.8 m/s. That is not a failure of the measurement: on a stationary pattern multiplied by a moving envelope, a single cross-correlation returns a variance-weighted blend of the two. The consequence is that cloud-drift speed is not a velocity in this regime, and any inference that treats it as one — including the height inversion of section 2a — is inferring from a mixture. The speed gate rejects the slowest columns, but 22% of them pass it.
The source term needs three things.
The wave structure is present. The forcing is a linear mountain-wave solution whose phase tilts with height, driven by a steering layer taken from 450–675 hPa, with a two-layer Scorer resonance and a closed-form trapped mode. It is not the boundary-layer upslope diagnostic that an in-phase terrain forcing would give. It reports a resonant wavelength of 10.1 km against 7–8 km observed. Two defects remain and both are measured. The amplitude is clipped: in the output field the vertical velocity reaches its ceiling at every level from 375 to 850 hPa and its ninetieth percentile sits exactly on that ceiling through 525–600 hPa, so the crests are flat-topped rather than sinusoidal. And the vertical decay is too weak: wave energy peaks at 575 hPa inside the steering layer, falls by roughly a factor of three by 300 hPa and then flattens near 1.2 m s-1 instead of decaying to zero, so a non-trapped residual reaches the tropopause.
The evaporating half of the couplet is bypassed rather than merely disabled. With saturation adjustment enabled, which is the shipped configuration, the explicit lift and subsidence terms are computed and then discarded; descent instead removes cloud physically, through conservation of potential temperature (a sinking parcel arrives warmer, saturation specific humidity rises, the diagnosed condensate falls). That mechanism is correctly wired but is outcompeted: across the amplitude ladder the correlation between vertical motion and cloud change strengthens monotonically but only to −0.14, and the mean cloud change in descending air remains positive at every lead. Descent slows cloud growth by roughly a factor of five; it is not yet a net sink.
The third requirement — the moisture field for the source to act on — is the one that does not hold up. A parcel crossing a wave conserves its total water while the partition between vapour and condensate oscillates, so total water should be smooth across a wave cloud. Measured on the 2 September case, it is not: the small-scale total-water perturbation correlates with the condensate at +0.32 to +0.39 through the wave-cloud layer (475–550 hPa), where a conserved formulation would give approximately zero, and its decorrelation length is only about three times the cloud's own, against an observed wavelength of 7–8 km.
The origin is the condensate control rather than the humidity control, which can be shown by comparing configurations. Where the humidity control is effectively inert the correlation is 0.75 to 0.92: specific humidity stays at its smooth background value, so the total-water perturbation reduces to the condensate itself and the correlation approaches one by construction. Where the humidity control is active it falls to 0.32 to 0.39. In other words condensate is placed without debiting the vapour that formed it, and the humidity analysis is left chasing that debt afterwards rather than the two being consistent by construction. The remedy is the change of variable below, which makes the inconsistency unrepresentable rather than something a second control has to correct. Box-mean total water over 400–700 hPa also falls by about 6% across the first five-minute interval, though that figure is not yet corrected for flux through the box boundary and should be treated as indicative.
The intended remedy is a change of variable, and it has well-established prior art. Cloud-resolving models of this class advance conserved quantities — liquid/ice water static energy together with total non-precipitating water and total precipitating water — and diagnose the vapour, liquid and ice split from them; this is the formulation used by SAM (Khairoutdinov and Randall, 2003). Our control vector does the opposite: four condensate species and specific humidity, each free to move independently, so nothing in the formulation prevents an analysis from placing condensate without the vapour needed to sustain it. Carrying total water as the analysed variable would make that inconsistency unrepresentable rather than merely penalised. Precipitating water would join the conserved set once sedimentation is treated, at which point the budget closes over the column rather than over the grid box. The constraint that shapes how much of this can be borrowed is that every operation must stay differentiable for the variational minimisation, so an iterative saturation adjustment is not directly usable and the smooth, analytically differentiable form used here is retained.
One approach was tried and does not work. Mapping cross-correlation speed against terrain to separate the crest and envelope components fails for a structural reason: a terrain-locked crest pattern is itself a stationary feature, so removing the static field to isolate moving cloud removes the crests along with it, and what remains is the envelope in both fields. Without that removal the correlation locks onto the land surface instead. The components do not separate by this route, and the terrain-relative error scoring above is what replaced it.
Update, September 2026: the separation does work, by scale rather than by speed against terrain, and the change of variable is now built. Band-passing the two window-time visible fields into octaves and measuring the translation gain in each — the gain being the correlation at the best along-wind shift minus the correlation at zero shift — separates the two components cleanly on this case. The 4–8 km band, which is the crest scale, gains nothing from being shifted: its gain is −0.023. Every other scale gains, and the gains grow with scale: +0.039 at 2–4 km, +0.039 at 8–16 km and +0.057 at 16–32 km. So the crest band stands and everything both larger and smaller advects through it, which is the structure asserted at the head of this section, now measured directly from the observations with no model involved. Two cautions belong with that number. The raw gain is biased upward, since taking the best of twenty-one trial shifts finds some correlation in noise; the figures above subtract the same statistic computed across the wind, where advection cannot contribute, and the bias so measured runs from 0.000 to 0.023. And the position of the peak carries no quantitative meaning: calibrated against synthetic mixtures of a standing and an advecting field, the peak jumps from one component to the other rather than moving smoothly between them, so it ranks the two regimes without measuring their proportions.
The same measurement on a second case, 23 August 2026 at 1525 UTC, returns no standing band at all, the crest scale least of any: +0.108 at 4–8 km against +0.056 and +0.038 at the two scales above it. That case is a growth episode rather than a steady wave, and the contrast is what makes the first case's single standing scale meaningful rather than an artefact of the method.
The change of variable is implemented. Total water and liquid water potential temperature are now the analysed pair, with temperature, vapour, liquid and ice diagnosed from them, and the forward model advects the pair and re-diagnoses the species at every step rather than only at the initial time. The spatially fixed source S(x) that this section argues for is also built, as a vertical motion field placed from the observed visible imagery: its shape is the along-wind derivative of the observed cloud, which makes it zero-mean by construction and so balances ascent against descent, and its amplitude comes from the saturation inversion. On the wave case it improves the unassimilated ten-minute lead, taking the 4–8 km visible correlation from 0.401 to 0.438 and cloud-top correlation over cloudy pixels from 0.755 to 0.777.
What the configuration is worth, against the deployed analysis and against the cloud analysis it starts from. All three are archived for the wave case, so all three can be scored against the same observations over the same box with the same forward operator. The cloud analysis is included by re-running with the minimisation switched off, which leaves the analysis at its background, and with the seeding switched off so that background is the cloud analysis itself rather than a primed version of it. Root-mean-square error against the observations, infrared in kelvin and visible in reflectance, over the whole box and over the observed-cloudy pixels alone:
| Lead | Run | IR RMSE (K) | VIS RMSE (reflectance) | ||
|---|---|---|---|---|---|
| all | cloudy | all | cloudy | ||
| 0000 (1355Z) | LAPS cloud analysis | 1.763 | 2.100 | 0.0381 | 0.0791 |
| operational 4D-Var | 1.926 | 2.609 | 0.0430 | 0.0854 | |
| test 4D-Var | 1.741 | 2.188 | 0.0332 | 0.0487 | |
| 0005 (1400Z) | LAPS cloud analysis | 1.816 | 2.390 | 0.0404 | 0.0843 |
| operational 4D-Var | 1.628 | 2.854 | 0.0381 | 0.0794 | |
| test 4D-Var | 1.600 | 2.561 | 0.0279 | 0.0556 | |
Table 7.2. Error against the observations for the cloud analysis, the deployed 4D-Var and the test configuration, at both window times, over the whole box and over observed-cloudy pixels alone. Both leads are fitted times, so this is a fit diagnostic.
Three things are worth reading off this, and the second is not in the variational analysis's favour.
The test configuration has the lowest error in the visible at both times, by a wide margin over cloud — 0.0487 against 0.0791 for the cloud analysis at the window start, and 0.0556 against 0.0843 at the analysis time — and the lowest infrared error over the whole box. Against the deployed analysis it is ahead on all eight comparisons.
But both variational analyses are behind the cloud analysis in the infrared over cloudy pixels, at both times. The deployed one by 24% and 19%, the test one by 4% and 7%. The cloud analysis is built to match the infrared cloud-top signal directly and it does that better than either analysis that starts from it, so the variational step is currently trading cloud-top temperature for the visible structure it gains. The test configuration gives away far less of it than the deployed one does, which is the direction of travel, but it is a real cost and it belongs in Table 7.2 alongside the gains.
At the window start the deployed analysis is behind the cloud analysis on all four measures. That is defensible given the cost weights the end of the window more heavily than the start, so the window start is the end the deployed configuration has least reason to protect; it is recorded here rather than left out because it bears on how these numbers should be read.
Both leads are fitted times and the table is a fit diagnostic, not a forecast result. The analyses are drawn to the observations at 1355Z and 1400Z, and the conserved pair with a placed vertical motion has more freedom in its control vector than the deployed species analysis, so fitting those times more closely is partly what the extra freedom buys. What the table establishes is that the freedom is spent on the observations rather than against them. Whether it survives to times that were never fitted is the separate question answered at the unassimilated leads. The visible margin also has a mundane component: the visible term carries the larger weight at the analysis time, so it is the term the minimiser works hardest on. The cloudy masks hold 5,744 infrared and 26,121 visible pixels at the window start.
As an independent check, the cloud analysis reports its own fit in its logs, computed over the full domain with its own forward operator rather than the one used here: 2.539 K infrared and 0.0250 visible at 1355Z, and 2.622 K and 0.0263 at 1400Z. Those are not comparable line-for-line with the table — a different domain, mask and operator move the visible figure in particular by more than the margins being discussed — but they confirm the cloud analysis is performing normally on this case rather than unusually well or badly.
Two of the three barriers to deployment have come down; the third has not.
Memory, which was the binding one, is resolved by tiling the analysis. The minimisation is the only phase that carries an adjoint tape — the forecast integrates without one — so only the analysis needs splitting, and its halo is set by how far information travels in one five-minute window rather than over a forecast. That distinction is worth about three times the overlap: roughly seven grid points of advection plus a five-point infrared footprint, against the thirty-four to sixty points a forecast-length cushion would demand. Tiling the 599×999 domain into eight pieces with a 24-point halo costs 23.6% extra work and peaks at 19 GB per tile against the 114 GB a single-piece analysis needs and the 20 GB the deployed system uses. One full-domain forecast then runs from the stitched analysis in 156 s at 20 GB, since it has no tape to hold. The complete cycle — eight tiles, stitch, forecast to +45 minutes — takes 22 minutes.
The analysis is better for being tiled, and the reason is physical rather than numerical. Several quantities the analysis treats as domain-wide are not uniform across it. The step normaliser is the largest gradient anywhere in the domain, and across a single 154×200 km box the two halves differ by a factor of 2.3, so every point in the weaker half is under-stepped by that factor. The steering wind varies 38% across the same box, and it sets both the amplitude of the inferred vertical motion, through the parcel transit time, and the direction of the dipole. Splitting the domain gives each region its own normaliser, its own line search and its own steering vector. Block-local normalisation recovers part of this without splitting anything, and it is now carried in the preset; the rest needs the split. Measured against a single tile that already has the block-local normaliser, on the wave case at 20 iterations with the error mask held 34 points inside the boundary:
| Lead | Run | IR RMSE (K) | VIS RMSE | VIS 9–14 km r |
|---|---|---|---|---|
| 0000, wave case | one tile | 1.841 | 0.0357 | 0.793 |
| eight tiles | 1.523 | 0.0329 | 0.853 | |
| 0005, wave case | one tile | 1.677 | 0.0314 | 0.918 |
| eight tiles | 1.666 | 0.0304 | 0.927 | |
| 0000, growth case | one tile | 3.206 | 0.0473 | 0.941 |
| eight tiles | 2.708 | 0.0442 | 0.948 | |
| 0005, growth case | one tile | 3.546 | 0.0580 | 0.898 |
| eight tiles | 3.292 | 0.0548 | 0.913 |
Table 7.3. Tiled against single-tile analysis at the two window times on both cases, with the error mask held 34 grid points inside the boundary. The single tile already carries the block-local step normaliser, so this isolates the split itself.
Splitting the domain did initially cost 3 to 8% in the infrared over cloudy pixels, and following that cost led to a defect with nothing to do with tiling. The cost was a tail rather than a general degradation — the worst one percent of cloudy points supplied half to two thirds of the increase in squared error and the top five percent more than all of it, so the other ninety-five percent were better — and those points sat in tile interiors, a median 73 to 85 grid points from the nearest boundary. That pointed at the line search rather than the stitching, and at the bound on its step, which had been left at the arbitrary default of three times the control scale. A five-point ladder on both cases puts the best value at one: averaged over both cases and both window times it improves the infrared by 4.5%, the visible by 1.9% and the wave band by 1.6%, and at the window start it is the best of the five on seven of eight measures. The bound is worth tightening whether the domain is split or not, and with it in place the cloudy-pixel cost of splitting is gone: at the window start the eight-tile analysis is now better than the single tile on that measure too, on both cases. One caution belongs with it: only the visible error at the window start varies monotonically with the bound. Everything else does not, because the bound changes which points the line search may move and so which solution it converges to, and a finer ladder would be measuring that rather than a trend.
Seams are a separate and smaller matter: the largest is a 0.24 K single-column step in temperature, and cross-fading the interiors over six grid points — free, since each tile already computes its halo — brings every seam below the typical column-to-column variation without measurably changing the scores. A fold of the cloud fields on the tile spacing shows no periodic imprint: fine-scale structure does run 16 to 35% higher within twelve points of a seam, but the untiled analysis shows a larger excess at the same columns, so that band is terrain and cloud rather than an artefact of splitting.
The iteration count is set by throughput rather than by the skill curve, and the shipping value carries a known cost. Measured end to end on both cases, the cycle takes 22.1 minutes of compute at twenty iterations and 16.0 at twelve. Against the single-tile reference across the sixteen window-lead measures, twenty iterations wins fifteen and twelve wins twelve; where twelve falls short it is the wave case at the analysis time, by 7.5% in the infrared and 13.9% over cloudy pixels. The preset carries twelve deliberately, the cheaper cycle having been preferred with that regression accepted and set aside for revisiting once more cases are available. Both configurations are archived and the choice is a single setting.
What sets the cycle time is contention, not the minimisation, which is why shortening the latter helps less than it appears. With the deployed five-minute cycle also running, the same work takes 31.5 minutes at twenty iterations and between 28.5 and 38.3 at twelve — the growth case at twelve took longer in wall clock than the wave case at twenty, on a third less compute. Per-tile times overlap almost entirely once queueing is counted, and the contention factor rises from about 1.45 to 1.78 as the compute shrinks, because fixed per-tile costs become a larger share. The levers that move the cadence are therefore standing the deployed cycle down while the tiled analysis runs, which gives 22 minutes, or an hourly update.
What remains is the radar term. It survives into the pair formulation but has almost nothing to act on, because rain and snow are taken from the background rather than carried in the control vector; switching radar off leaves the total-water and vertical-motion gradients unchanged to four significant figures. That is a loss of capability rather than a wiring detail, and it is the outstanding item.
Taken together over both cases, all four lead times and four measures, the configuration is better than a single tile on 26 of 32 comparisons and on 15 of the 16 at the two window times. The exception at the window times is cloudy-pixel infrared on the wave case at five minutes.
Three limitations belong with all of this. Everything above rests on two cases. The forecast leads are re-seeded from the written analysis rather than continuing in memory, which on the wave case costs about 10 to 13% of infrared error at ten and fifteen minutes while improving the visible by about 5% and the wave band by 3 to 5%; a control seeded from an untiled analysis through the same path loses the same amount, so the cost belongs to re-seeding rather than to splitting, and the growth case does not pay it at all. And the 22-minute cycle assumes the machine: with the five-minute deployed cycle also running, contention stretches it to 31.5 minutes, so a half-hourly update needs either a shorter minimisation, or the deployed cycle to stand aside on the half hour.
The case is archived with its verifying observations and its operational output. Nothing here rests on a single number, but all of it rests on a single case.
Kurzrock, F., S. Cros, F. Chane Ming, J. A. Otkin, A. Hutt, L. Linguet, G. Lajoie and R. Potthast, “A Review of the Use of Geostationary Satellite Observations in Regional-Scale Models for Short-term Cloud Forecasting”, Meteorologische Zeitschrift 27(4), 277–298 (2018). The review that frames this problem directly, and the closest survey of the family this system belongs to. Two points bear on the claims above: the field norm for the horizon over which satellite cloud assimilation demonstrably helps is measured in hours, against the 25–30 minutes quoted in section 4, so that horizon is a genuine weakness rather than a conservative statement; and the review identifies the treatment of vertical motion in cloud assimilation as an open gap, which is the gap section 7g addresses.
Satoh, M., B. Stevens, F. Judt, M. Khairoutdinov, S.-J. Lin, W. M. Putman and P. Düben, “Global Cloud-Resolving Models”, Current Climate Change Reports 5, 172–184 (2019), doi:10.1007/s40641-019-00131-0. Useful here for perspective rather than method. Global cloud-resolving models resolve nonhydrostatic motion at kilometre-scale meshes and abandon cumulus parameterization entirely, simulating cloud with explicit microphysics; the DYAMOND intercomparison it describes runs eight of them. They are solving a different problem from this one — global, free-running, and judged on convective and climate statistics over forty days, against a limited-area variational analysis judged on a sub-hour nowcast of a single orographic case. What transfers is the standard they set for the thermodynamics: at these scales the vapour/condensate partition is expected to follow from resolved dynamics and microphysics, not from an analysis increment placed by hand, which is precisely the defect identified in section 7g.
Khairoutdinov, M. F. and D. A. Randall, “Cloud Resolving Modeling of the ARM Summer 1997 IOP: Model Formulation, Results, Uncertainties, and Sensitivities”, Journal of the Atmospheric Sciences 60, 607–625 (2003) — the System for Atmospheric Modeling (SAM). The prior art for the change of variable proposed in section 7g: SAM advances liquid/ice water static energy together with total non-precipitating and total precipitating water, and diagnoses the vapour, liquid and ice split from those conserved quantities. The limit on what can be borrowed is that this system requires every operation to remain differentiable for the variational minimisation, so an iterative saturation adjustment is not directly usable.
Hakim, G. J., J. S. Whitaker, B. Huang and S. Frolov, “Long-window 4DVar for reanalysis using a differentiable weather model”, arXiv:2608.11515v1 (2026). The nearest published relative of the method used here: a variational analysis built on a differentiable model and minimised through automatic differentiation, with no separately maintained tangent linear and adjoint codes. Two of its choices bear directly on this system. It minimises with AdamW rather than conjugate gradient, on the argument that conjugate gradient assumes a quadratic cost surface and a nonlinear model does not provide one — a point that applies with force here, where saturation adjustment and the subgrid distribution switch condensation on and off and the total-water gradient carries a max/p90 ratio near 1200, and where the minimiser is still conjugate gradient. It also drops the background term altogether for windows of two to seven days, replacing it with a penalty on the first model step. That matters here because the same balance turns out to hold by accident rather than by design: measured on the conserved-pair configuration, the total-water background term sits around 10-8 of the total cost, so this analysis is already running with a background constraint that carries no weight, on a five-minute window. The distance is otherwise considerable — a global 2.8° reanalysis assimilating surface pressure alone over multi-day windows, against a limited-area cloud nowcast fitting satellite radiances over five minutes — so what transfers is the treatment of the optimisation, not the configuration.
All figures are single-case measurements unless stated otherwise. The combined figure of merit weights visible error roughly 2:1 against IR in practice. At +5, +10 and +15 the 4D-Var run is ahead in both channels separately, so those margins keep their sign under any reweighting. The analysis row does not: there the visible fit is marginally worse than the background's (by 0.004 reflectance) against a 1.40 K infrared gain, so a visible weighting more than about three times the present one would reverse that row alone.