Steve Albers' LAPS Adjoint-Based 4-D Variational Cloud and Hydrometeor Assimilation and Short-Range Forecast at Sub-Kilometer Scale

This page describes the variational cloud / hydrometeor assimilation and short-range forecast 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. Most of the testing so far covers leads to 20 minutes from the window start, the range of the real-time forecast; we are now extending the evaluation to 60-minute leads on archived cases, beginning with the 5 October case of section 5c.

At a glance

What it does. Every 30 minutes the system takes the latest LAPS cloud analysis and a short-range forecast from the HRRR (High-Resolution Rapid Refresh) model as its starting point, and adjusts the three-dimensional cloud, moisture and temperature field until a 5-minute run of a cloud model reproduces what the GOES-19 satellite and the radars observed at the start and the end of those 5 minutes. It then runs the same model forward to forecast the cloud: 15 minutes ahead in real time, and 60 minutes on archived cases. By day it runs on a 500 m grid over the Colorado Front Range; at night, when only infrared imagery is available, it runs on a 3 km grid over a wider area with a 10-minute window.

Why a model and not just a moving image. The usual short-range cloud nowcast moves the latest satellite image along its observed motion. That image is two-dimensional: each pixel gets one motion, the motion of whatever the satellite sees on top. Real cloud is three-dimensional. Layers at different heights move with different winds, a lower deck under cirrus moves on its own while the satellite sees only the cirrus, and a deep storm spans levels whose winds differ in speed and direction, so its anvil and its precipitation core drift apart (on 25 July 2026, 10.5 against 3.2 m/s; section 3a). The analysis keeps cloud at its own height in every column, including a lower deck the first guess holds under the cirrus, and moves each level with that level's wind, blended toward the storm's own motion in precipitation cores.

Growth, decay and terrain. A moving image also keeps the cloud's appearance, so the cloud can't grow, decay or stay anchored to terrain. Here the cloud is the condensed part of a moist air mass: it forms where air is lifted to saturation and evaporates where it sinks, so mountain-wave cloud stays in place while the air flows through it. On the Boulder lee-wave case the analysis matches the observed visible image with a correlation of 0.97 at the analysis time, against 0.91 for the moved image (Table 7.4).

Where it stands. The fit at the analysis time is close on every case evaluated so far. The forecast is the harder part: the steering winds come from the background model and stay fixed (only the inflow and outflow around rising air follow the evolving vertical motion), and the cloud life cycle is simplified, so a storm anvil that is decaying in reality can keep spreading in the model. Section 4 lists the limitations and section 8 the work under way.

Reading guide. Sections 1 and 2 describe the method and the physics as they run now; Table 1.1 is the one-page summary of the real-time configuration. Section 3 compares the system with moving the image, and section 4 lists its limitations. Sections 5 to 7 treat three problems in detail: the overshooting storm anvil, the ice-phase correction and mountain-wave cloud. Section 8 is future work, section 9 the case studies, section 10 related work and section 11 the references.

Terms used on this page

Table 0.1. Terms used throughout the page.

1. The technique

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 one set of control variables at the window start, and the adjoint returns the gradient of the total cost with respect to those controls.

1a. Control variable and cost function

The control variable is the thermodynamic state at the window start on the full 3-D grid: liquid-water potential temperature θ_l and total water q_t (the conserved pair), together with the vertical velocity w. Cloud liquid, cloud ice, vapour and temperature are not controlled directly; they are diagnosed from the pair by saturation adjustment, so condensate cannot exist without the water that makes it, and its latent heat is accounted for exactly. Rain, snow and graupel are taken from the background and carried through the model, and graupel contributes to the reflectivity operator. An earlier formulation controlled each hydrometeor species directly; several of the measurements in sections 3 to 6 were made with it, and they are dated so this is clear.

Writing x₀ for that initial state and x_b 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

• J_b sums the squared departure from the background separately for each control, each scaled by its own background error (0.5 K for θ_l, 1 m/s for w), so temperature, water and vertical motion move by different amounts.

• J_ir and J_vis compare rendered radiances against two channels of the GOES Advanced Baseline Imager (ABI): clean-IR window brightness temperature and visible reflectance factor. Both channels are parallax-corrected, so cloud-top radiances are mapped to the ground location they belong over rather than the sub-satellite line of sight, and both carry per-pixel observation errors. The visible term enters only where the sun elevation exceeds 11°.

• The time weight w(t) rises at the window end, since the analysis time is the deliverable. The visible term carries an additional end-time factor of its own, so in operational use the effective visible weight at the analysis time runs three times the IR weight.

• J_z compares modelled and observed reflectivity in Z-space, against the observed mosaic rather than against an analysed liquid-water field. A power transform on Z handles the enormous dynamic range between echo and no-echo, and the targets are held sparsely, over echo points only. This term carries its own weight rather than the time weighting, since the radar target sits at the end time.

• J_rh charges condensate held in air too dry to support it, weighted by the condensate-to-saturation ratio, with a temperature-dependent phase reference (section 2a). J_ph charges precipitation of the wrong phase for its wet-bulb environment. A moisture lid holds total water at its first-guess value above the tropopause (section 7c).

1b. Forward model

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 conserved pair and the precipitation are carried by semi-Lagrangian advection on the effective wind field described in section 3a, sub-stepped to respect the Courant limit, and moved vertically by the analysed w. After each step the cloud, vapour and temperature are diagnosed again from the pair, so condensation and evaporation follow the motion. The trajectory weights are computed once and shared, so cloud and precipitation in the same air move together. The free forecast beyond the analysis time runs the same model, with its physics. At each step the horizontal wind keeps the background's rotational part and takes its divergent part from the current vertical motion (section 2b). In the forecast, air arriving from outside the domain carries the analysed state at the edge it enters through, so a storm feeding the edge at the analysis time keeps feeding it; within the window, five minutes of travel stays inside the overlap between tiles. Section 1e lists what is switched on in real time and the recommended presets (Table 1.2) and where each data source acts (Table 1.3), and section 2 sets out every process and its status.

1c. Minimization and gradients

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 nine iterations, arranged in three blocks, one for each control: θ_l, q_t and w. Each block runs its own line search, because their gradient magnitudes differ by many orders of magnitude and a shared step length serves none of them. Within a block the step length is chosen region by region, from each 24 km region's own fit to the satellite (section 7d), and the cost function is compiled once per tile for speed. The domain is analysed in six overlapping tiles and the forecast runs once over the whole domain (section 7d).

1d. The first guess and its preparation

The first guess is the LAPS cloud analysis for the window start, re-read every cycle; the analysis never starts from its own previous output. Before the control variables are built, a short sequence of steps corrects errors in that first guess which the minimization cannot correct by itself, because the observation gradient that would do it either vanishes or would act at the wrong height. Each step changes the background and the starting state together, so the background term does not pull the analysis back toward the uncorrected field. In the order they run:

• Layer repartition. The observed cloud motion is measured locally from the satellite images, and a fraction of the condensate below about 8.5 km is moved to the level whose wind matches that motion, keeping the total mass exact and converting liquid lifted into cold air to ice. The motion is measured once over the full domain, which prevents tile edges from losing part of the field, and from the visible rather than the infrared image, with no thin-cloud limit, so optically thick cloud can be moved too. On the 24 September 2026 high-cloud case this places 23% of the visible optical depth above 450 hPa and improves the visible correlation of the forecast by 0.01–0.02 at 5–20 minutes, with the infrared unchanged. The visible is used because it moves with the layer that holds its optical depth, while the infrared moves with the cloud top.

• Infrared trim. Where the first guess renders more than 3 K colder than the observed brightness temperature, cloud at levels colder than the observation allows is removed. The minimization could only warm such a column by thinning it, which also dims the visible; the trim removes the cloud at the height the infrared rules out instead.

• Visible fit. Column by column, the first-guess condensate is rescaled so that the same visible operator used in the cost reproduces the observed reflectance at the window start. Where the first guess is clear but the visible shows cloud, new cloud is placed hanging from the highest level the observed infrared allows, over a 150 hPa depth, with the air brought to saturation so the conserved-pair diagnosis keeps it. The fitted optical depth is capped so that observations at the operator's low-sun brightness limit do not demand an unbounded column. This step exists because the visible gradient vanishes at both ends of the reflectance curve: in clear, subsaturated air, where a small moisture change produces no condensate, and in optically thick cloud, where reflectance no longer responds to optical depth. Because the conserved-pair analysis diagnoses cloud from total water and saturation, the fit sets the humidity as well as the condensate, point by point, so that the model's own saturation adjustment returns exactly the fitted cloud, which prevents the humidity pattern from overriding the fit. On the 24 September 2026 high-cloud case the analysis with this step matched the observed visible with a correlation of 0.91 at the window start, and its forecast beat persistence of the observed image at 5, 10 and 15 minutes (0.86, 0.80 and 0.71 against 0.84, 0.75 and 0.70), with a lower error than persistence at every lead to 20 minutes and a better infrared fit. It has run in real time since 26 September 2026, together with the observed motion that places relocated cloud.

• Wave-cloud seed and balanced start. Where the visible shows cloud standing still over rugged terrain, a vertical motion is placed from the observed crests, moisture is staged in the rising air, and the temperature and total water are displaced along the wave's streamlines so that advection and lifting cancel (section 7c).

• Conversion to the conserved pair. Liquid-water potential temperature and total water are formed from the prepared state, and during the minimization a moisture lid holds total water at its first-guess value above the tropopause, except where a tall radar echo and a very cold observed top indicate an overshooting storm.

• Hydrometeors kept above ground. Following the rule of the LAPS cloud analysis, cloud is allowed only above the terrain and precipitation down to the first pressure level below it. The pair keeps any water at underground levels as vapour, so total water is unchanged, and precipitation that falls below that level leaves as surface precipitation (section 7g).

The first guess takes its horizontal placement from the satellite but its vertical structure from the cloud analysis's assumption that the top is opaque, and that is where it is weakest: on the 5 October case, thin cirrus was missing aloft in 70–75% of the pixels where the water-vapour and window channels together place a thin top, and where present it sat about 200 hPa too low. Section 8c describes a planned first guess that takes its vertical structure from the background model's own columns instead.

1e. What runs in real time

The table gathers the real-time configuration in one place. Processes that are switched on but have no effect on this path are listed as such, so a switch is never mistaken for a process that is acting.

Component	In real time
Controls	θl, qt and w at the window start; rain and snow from the background
Window	5 minutes, two observation times (10 minutes at night on the 3 km grid); visible wherever the sun is more than 11° above the horizon at each pixel, rendered at that pixel's own sun angle, and infrared at both, radar reflectivity in Z-space, weighted so that the radar term is comparable to the satellite terms where there is echo; the 6.9 µm water-vapour channel at the window start, in clear and cloudy air (16 October, section 5c)
First-guess preparation	infrared trim, visible fit with observed motion, wave-cloud seed and balanced start, moisture staged at the wave ascent (RH 0.80), a moist halo (RH 0.80) around the fitted wave cloud, below-ground mask; layer repartition acts on about 3% of columns. Convective vertical motion comes from a parcel formed from the lowest levels above the ground and lifted from the top of that layer; it rises dry below its condensation level, stops once its buoyant energy is spent, and reaches no higher than the cloud top the infrared shows
Condensation and evaporation	on, by saturation adjustment of the pair with a subgrid ramp beginning at 85% relative humidity, in the window and the forecast
Latent heating of the air	on: temperature is diagnosed from the conserved pair, so condensation warms and evaporation cools
Temperature	one analysed temperature (16 October): the window physics and the infrared render start from the temperature the analysed pair gives, and the forecast starts from the window-end temperature, vapour and condensate together (section 5c)
Vertical advection	on, vapour and condensate by the evolving w, with one anti-diffusive (MPDATA) correction under a non-oscillatory limiter (7 October), so no cell leaves the range of its vertical neighbours; nothing enters through the ground or the model top
Vertical motion	evolving, the same in the window and the forecast (9 October): a part anchored to the terrain and the standing wave, plus a convective part that is carried by the wind, relaxes toward the parcel ascent over 60 minutes, is projected onto mass continuity, and is capped at 3 m/s. The forecast starts from the first-guess vertical motion plus half the analysed increment (16 October, section 5c). The analysed increment is correlated in height over 1.5 km and, in cloudy columns, also takes the shape of the first-guess cloud layer, so ascent forms deep, coherent cells; the adjoint sees how the parcel ascent responds to the state it starts from (20 October)
Ice phase	the ice share of new condensate follows the current temperature, and ice is judged against ice saturation (5 October); cloud ice falls at a speed set by its content (9 October)
Sedimentation and autoconversion	on; cloud turns to precipitation above 0.5 g m-3, so low stratus does not drizzle where the radar sees no echo
Adiabatic lifting source, subsidence evaporation	switched on but bypassed: saturation adjustment already condenses rising air and evaporates sinking air, so they would count it twice
Accretion, precipitation evaporation	off
Winds	the rotational (steering) part from the background model, blended toward the radar steering vector in precipitation cores (section 3a) and held fixed through the forecast; the divergent part follows the current vertical motion through continuity, within about 20 km of the ascent, in the window and the forecast (20 October, sections 2b and 5c)
Domain edges	in the forecast, inflow takes the analysed state at the edge, held through the forecast (20 October); the window's five minutes of travel stays inside the tile overlap
Minimization	nine iterations, three blocks, regional step length chosen on the observation terms, which a block may take only when it does not worsen the fit to the observations (9 October), six tiles
Forecast	15 minutes beyond the analysis time (20 minutes from the window start), same model and physics; the simulated visible of each lead is rendered under the sun of its own valid time (20 October, section 2a); archived cases are now also run to 60 minutes for evaluation

Table 1.1. The real-time configuration as of 20 October 2026. "Bypassed" marks processes switched on whose effect the saturation adjustment already supplies.

Recommended combinations of settings. The analysis program has several hundred options, and the useful combinations are kept as named presets, selected with --preset. The configuration in Table 1.1 is the preset realtime: the real-time cycle runs it unchanged, adding only launch settings (tiling, logging and cadence). A test meant to match real time starts from --preset realtime and adds just the setting under test, so any difference in the result comes from that setting alone. Table 1.2 lists the presets and when to use each.

Preset	What it is	Use it for
realtime recommended	The real-time configuration of Table 1.1: pair-standing plus the first-guess preparation, regional step length and physics safeguards added from 22 September to 20 October 2026	Every test that should match real time, and any rerun of an archived real-time case
pair-standing	The conserved pair (θl, qt) with w, the standing-wave seed and the lee-wave dipole; realtime is built on it, so a change here reaches real time as well	Comparisons with the conserved-pair line before the September additions
operational-zspace	The earlier species minimizer, with radar in Z-space	Reproducing results from before the conserved pair
operational, radar-zspace, radar-off	The July 2026 species-minimizer configurations, with and without radar	Reproducing July results
subset-quick	A small cloudy subset, in seconds	Quick checks while developing code

Table 1.2. Presets: the recommended combinations of settings. Settings given explicitly on the command line take precedence over the preset, and each run prints the settings its preset supplied.

Where each data source acts. The observations enter the analysis at two stages. In the preparation of the first guess, each source places, trims or seeds the starting state directly. In the minimization, the analysis adjusts the conserved pair (θ_l, q_t) and w so that the rendered observations match at both window times. The minimization is well suited to refining amount, height and small displacements, where the cost has a useful gradient. Creating cloud where the first guess has none, placing a standing wave, and removing condensate far from its observed position give the minimization little gradient to work with, so the preparation handles them. Table 1.3 shows how each source divides between the two stages.

Source	Preparation of the first guess	Minimization
Background model (HRRR or NAM, via LGA)	temperature, humidity and winds for the window; rain and snow taken as they are; rotational winds held fixed, blended toward the radar steering vector in precipitation cores	background term: the analysis stays close to the prepared first guess
LAPS cloud analysis	starting cloud field	cloud constraint at the end of the window
Infrared	removes first-guess cloud that renders too cold; sets the layer allowed for wave cloud; limits convective ascent to the cloud top the infrared shows	rendered brightness temperature fitted at both window times
Visible (by day)	fits the first guess column by column; the image sequence supplies the motion field; places and seeds the standing wave and its moist halo	rendered reflectance fitted at both window times
Radar	moisture lid from the echo top; steering winds in precipitation cores	reflectivity at both window times, with echo targets and an upper bound where the radar sees no echo, weighted to be comparable with the satellite terms where there is echo; rain and snow themselves are not controls in real time (8b)
Water vapour (6.9 µm)	none	rendered brightness temperature fitted at the window start only, clear and cloudy, with the error allowed growing with the infrared cloud signal (section 5c)
Terrain	below-ground rules; terrain-forced vertical motion	below-ground rules in the model

Table 1.3. How each data source divides between the preparation of the first guess and the minimization in the real-time configuration. The free forecast uses no observations; it carries the analysed state and vertical motion forward, and takes the analysed state at its inflow edges.

1f. Correction of cloud motion

Cloud in the model moves with the wind at its own level. Every level is carried by the three-dimensional background wind, in the window and in the forecast, and precipitation cores are blended toward the radar steering vector (section 3a). Getting the motion right therefore means two things: moving cloud has to sit at the level whose wind moves as the cloud is seen to move, and standing wave cloud has to be held by standing lift while the air passes through it.

What the minimization contributes. The cost compares the rendered satellite and radar images at both ends of the window, and the model carries the state between them. The adjoint of that transport takes a misfit at the analysis time back along each level's own wind to the window start, so the minimization corrects amounts, edges and displacements of a few pixels in the natural 4D-Var way. The vertical motion is itself a control variable: seeded where the observed crests stand still, it is adjusted so that condensation on the ascent and evaporation on the descent keep each crest in place while the air flows through it (section 7a).

Moving cloud to a different level is harder for the minimization. In principle the same adjoint signal carries it: cloud at the wrong level ends the window displaced, and the misfit points to the level whose wind would have carried it correctly. In practice the signal is the wind difference between levels times five minutes, often a few kilometres and comparable to the visible footprint, and moving condensate to a new level means creating cloud in unsaturated air, where the gradient of the saturation adjustment vanishes. Measured on the moving high cloud, the minimization changes the share of condensate in each layer by at most two percentage points, and that small change is downward rather than toward the faster winds above.

The first-guess preparation therefore chooses the height by the same visible fit the cost function measures. The observed cloud motion is measured from the two visible images of the window (infrared at night) over the full domain, with a confidence for each column. Where the motion is confident and the cloud's present drift misses it by 3 m/s or more, the cloud is placed in a layer centred on the level whose wind matches the observed motion, kept within the range the infrared allows, and the humidity is set so the model keeps it there. This is the choice that minimizes the analysis-time visible misfit under advection, made directly rather than by following a gradient. Where the visible stays put over rugged terrain, the column is treated as standing wave cloud and receives standing lift instead. The minimization then refines amount, edges and small displacements around the height the preparation has set.

How well it works. In the real-time cycles since 26 September the height the analysis sees mostly agrees with the level whose wind matches the observed motion (section 7f); on days of fast, extensive high cloud the harder limit is the wind, with no level matching the observed motion in up to two thirds of cloudy columns. On the optically thick high deck of 24 September the analysed cloud top, near 300 hPa, agrees with both the infrared (about 315 hPa) and the level whose wind matches the observed motion (275–300 hPa). What still slows that deck in the forecast is its depth, already present in the first guess: a third of its condensate lies below 400 hPa, where the wind is 5 to 7 m/s against the 12 m/s the tops show, and the visible sees deep enough into the deck to follow that slower cloud. A longer window, a control on each column's cloud height, or a term comparing the observed motion with the wind at the level the visible sees would let the minimization take on more of this itself.

2. Physical processes in the analysis and the forecast

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. The diagnostic physics is active every cycle. Since the conserved pair went into real time, the core prognostic processes are active too: condensation and evaporation, latent heating, vertical transport, sedimentation, autoconversion and the fall of cloud ice, and since 9 October the vertical motion evolves through the window and the forecast. Since 20 October the divergent part of the horizontal wind follows that vertical motion. Growth processes for ice and graupel, and a momentum equation for the winds, remain absent. The three subsections below follow that split, Table 1.1 gathers what runs in real time, and section 7e measures the physics against a physics-off run and simple image advection.

2a. Diagnostic processes in the observation operator

• Mixed-phase partition. A temperature-dependent ice water fraction mirroring the LAPS cld_iwf formulation exactly (ramp between −40 °C and −1 °C), used to split condensate between liquid and ice for both radiative and reflectivity purposes.

• Phase-aware saturation. Saturation vapour pressure blends the liquid formula with Goff–Gratch ice saturation by that same ice fraction, so the humidity term stays physically referenced at cirrus temperatures. Detail below.

• Variable effective radius, Rayleigh scattering, and a convective cloud-fraction diagnosis in the visible reflectance operator, so sub-grid brightness departs from plane-parallel.

• Instrument footprint convolution on both channels. The visible channel is rendered at 500 m but the instrument resolves about 1.2 km, so the rendered field is convolved to the footprint with three cascaded box passes, matching the treatment the IR channel has always had. This prevents the simulated visible from showing sub-footprint structure no instrument can see, such as isolated dark pixels, and the fit to observations improves steadily with the number of passes. The observed field must be convolved too when the two are compared: the stored visible observation is on the 500 m grid and carries sub-footprint resampling content, so convolving only the model flatters the sharper field and makes a correctly convolved model look over-smoothed.

• Anisotropic reflectance factor. The visible operator applies the same bidirectional reflectance correction the LAPS cloud analysis applies, so the two systems render identical cloud identically. This prevents a difference in angular reflectance from being mistaken for a difference in the cloud field.

• The sun of each image's own time. The cost renders the visible at each pixel's solar altitude, above the 11° limit of the visible term. The forecast product renders every lead under the sun of its own valid time, down to the true altitude, with the Rayleigh, ozone, aerosol and anisotropic-reflectance terms recomputed for that sun, so the forecast reflectance factor is quantitatively comparable with the observed one through sunset. On the 5 October case the mean simulated reflectance at 60 minutes is 1.24 times the observed, against 2.07 when every lead used the analysis-time sun. The operator is column by column and carries direct sunlight only; the diffuse sky light and cloud shadows that matter at low sun are future work (8c).

• Temperature-dependent Z–M relations for snow, and for graupel/hail a temperature-dependent weight following the LAPS dry/wet distinction. A single warm-line relation renders cold anvil snow roughly 8 dB too bright.

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, which prevents a kink that would stall 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. On the current real-time configuration its gates keep about 3% of columns, mostly because of the terrain ceiling described below; the visible fit's use of observed motion (section 1d) now does most of this work. 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.

2b. Prognostic processes in the model

All of the following are implemented as elementwise operations so the adjoint flows through them, and each has its own switch. The real-time status of each is given in the second column, and summarised in Table 1.1.

Process	Switch, and real-time status	Form	Note
Sedimentation	--precip-sink on (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 on, window and forecast	Vapour and condensate carried by the analysed vertical velocity in flux form, with one MPDATA anti-diffusive pass	The w it uses is a control variable, seeded from the observed wave crests where they stand still
Cloud condensation / evaporation	--sat-adjust, --fcst-condense on, window and forecast	Saturation adjustment of the conserved pair, with a subgrid ramp that begins condensing at 85% relative humidity. An older relaxation form remains for the species formulation	This is the sink that erodes a spreading anvil edge
Adiabatic lifting source	--lift-source, --k-lift on but bypassed by saturation adjustment	dqc/dt = −w ∂qsat/∂z for saturated ascent, applied as a direct source so it scales linearly with w	The linear form prevents saturation at small ascent rates, so stronger terrain-scale ascent produces proportionally more condensate
Subsidence evaporation	--oro-couplet, --oro-k-subs on but bypassed by saturation adjustment	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
Evolving vertical motion	--w-prog, --w-prog-continuity on, window and forecast	A convective part carried by the wind and relaxed toward the parcel ascent (60 min), projected onto mass continuity by an isobaric Poisson solve, added to a part anchored to terrain and standing waves	One routine serves the window and the forecast, so the analysis is fitted with the motion the forecast then uses
Divergent wind	--div-wind, --div-wind-screen-km 20 on, window and forecast	Per level, the wind keeps its rotational part and adds the gradient of a potential chosen so that the horizontal divergence balances the vertical motion; a screening length of 20 km confines the correction to the storm scale	Rising air drives outflow aloft and inflow below, and both weaken as the ascent does; in the window the adjoint carries the vertical motion through to the outflow it drives
Cloud-ice fall	--ice-sediment on	Fall speed from the ice content (Heymsfield and Donner 1990)	Removes ice that would otherwise only advect
Autoconversion	--precip-sink, --k-auto on	Cloud → precipitation above a content threshold, separate thresholds for liquid and ice	—
Accretion	--accrete, --k-accr off	Collection of cloud water by falling precipitation, Kessler-scaled (Kessler 1969)	At a defensible rate it is statistically indistinguishable from being off
Precipitation evaporation / sublimation	--precip-evap off	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, with the switch for each, its real-time status and its form. "Bypassed" means switched on but superseded by the saturation adjustment.

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 a switch is only live if the term or process it feeds is on the path actually being run: the conserved-pair path bypasses the lifting source and subsidence evaporation, and each run prints a line for every such switch so a bypassed one is never mistaken for a tested one.

2c. Prognostic processes still to add

• Graupel and hail generation. Graupel appears in the reflectivity operator and advects as a passive tracer, while its creation, removal and fall all remain absent. This is the leading candidate for the complete absence of core intensification documented in section 4, since riming is what builds a convective core's reflectivity.

• Melting and the bright band. No melting layer, so the radar bright band is neither represented nor usable as a constraint.

• Aggregation, depositional growth, and ice nucleation as distinct processes. Ice mass is carried, while the pathways that would grow or convert it stay absorbed into the autoconversion and condensation terms above.

• Latent heating feedback on the wind field. The rotational winds are prescribed and held fixed through the forecast. Condensation warms the air, through the conserved pair, the vertical motion evolves, and the divergent wind follows it, so a storm's inflow and outflow weaken as its ascent does. There is still no momentum equation: the winds respond to the vertical motion only through continuity, and heating, pressure gradients and rain-cooled downdrafts don't act on them (section 5). The steering-wind blend is explicitly not a mass-consistent wind, and is labelled as such in the code to keep it distinct from dynamics.

Section 8f describes the path from this kinematic forecast toward evolving dynamics.

3. How this goes beyond a 2-D advection nowcast

The measurements in sections 3 to 6 were made in July and August 2026 with the earlier formulation that controlled each hydrometeor species directly, and with less physics in the forecast. The mechanisms they establish still apply; the numbers belong to that configuration. Section 7e repeats the central comparison, against image advection, with the current one.

3a. The advection is a hybrid: 2-D through storm cores, 3-D at the periphery

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.

3b. What the analysis adds

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.

3c. Comparison against a forecast started from the unassimilated cloud field

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.

The two forecasts share one code path and differ only in their starting state: the same steering-blended winds, five-minute frames and condensation physics, and each lead stamped with its true valid time. That keeps the comparison controlled, so the difference is the assimilation and nothing else. The seeding switch that stamps the valid time correctly is off by default and has to be set for this comparison.

4. Current limitations

• The advantage narrows with lead. Section 3c resolves this by lead: the assimilation earns its cost at the analysis time and is at its largest at +5 minutes, the range that matters most for very short-range forecasting, then narrows steadily through +10 and +15 as the initial-condition advantage is carried away by a forecast that was, at the time of that measurement, advection alone. A mean over +5/+10/+15 averages that structure away, so the per-lead figures are the ones to quote. The narrowing trend is the real limitation.

• Terrain-anchored cloud is translated downstream. The forecast integrates condensation physics and carries the diagnosed mountain-wave vertical velocity, yet cloud held in place by the terrain is still carried downstream while the observation keeps it anchored. This was tested directly by scaling the wave forcing over an eightfold range (0.25 to 2.0). The pattern displacement between the two window times stayed at 3.4 km at every rung — the full wind — while amount and peak improved monotonically and placement against the observation degraded. The amplitude of the vertical motion is therefore an amount lever, not an anchoring lever, and the translation is set by horizontal advection alone. Measured on the 2 September 2026 wave case, the forecast falls behind persistence over high relief by 1.44 K at +15 minutes while beating it by 3.19 K over flat ground in the same domain at the same lead. The effect is not generic to mountains: on an ordinary convective case the ordering over high relief is the other way round and the margin grows. Section 7a gives the measurements and the missing term. Section 7c adds a direct test: cloud placed on the observed crests at the start of the window is carried off them within five minutes at every wave phase and amplitude, although lift and descent do make and evaporate cloud as they should.

• The analysis carries almost no structure below about 3 km. Comparing radial power spectra against the visible observation, the rendered field holds 64% of the observed power at 10 km, 22% at 4 km, 4% at 3 km and essentially none below 2 km, while the observation carries real power to 1 km. Part of this is the instrument footprint convolution of section 2a, whose cascaded box kernel reaches an effective width about 39% broader than the footprint it is meant to represent, because the physically correct width falls between the odd pixel counts the kernel can express. But the deficit is larger than that kernel alone accounts for, so the increment itself is spatially smooth and the kernel is not the whole cause. A root-mean-square error cannot see this at all, which is why it went unnoticed: it is a distribution of variance across scales, not a magnitude.

• Analysis fit is the wrong acceptance metric for anything that moves cloud. Measured repeatedly through August 2026: configurations that improve the analysis figure of merit can be worse at every forecast lead, and configurations that cost analysis fit can win at +10 and +15 minutes. On one stratiform case the variant with the best analysis fit was the only one harmful at every lead. Lead-resolved scoring costs no additional runs and is now required before any change that redistributes condensate is accepted.

• Convective cores fail to intensify. On 28 July 2026, observed ≥40 dBZ area grew 47% in 20 minutes (2,304 → 3,378 px) while the model went from 2,931 to 2,897 px, a 1% decline, and the increment was inert there (±50 px). The model also starts about 27% over-strong in core area and ends under. This is a missing-physics problem rather than an analysis or radar-weighting one: tuning the observation weights can only redistribute intensification the model already represents.

• Moderate echo under-grows. Observed ≥20 dBZ area grew 61%; the model grew 39–43%.

• Cloud-top temperature drifts warm with lead. IR bias runs −2.35, −1.00, +0.40, +1.71 K across the four leads: the analysis carries too much cold cloud and the +15 minute forecast too little. Anvil area spreads while the tops fail to stay cold and high, so area and depth are decoupled.

• The free forecast still trails image advection. It now carries condensation, evaporation, vertical transport and sedimentation (Table 1.1), and the physics helps on terrain-driven cloud, yet on all three cases of section 7e simple image advection is better beyond the window. On the first cases studied the main cause was cloud height, with the visible cloud below the level whose wind carries it. In the real-time cycles since 26 September the height mostly agrees with the wind, and on days of fast, extensive high cloud the limit is often the wind itself, with no level matching the observed motion in up to two thirds of cloudy columns (section 7f).

• The advection scheme fails to conserve mass. See the next section; this is a significant limitation of the present forecast.

• The analysis over-dries clear air aloft. With no observation of humidity in clear air, the minimization can lower it freely, and in the clear upper troposphere it does so to zero. On the 5 October case the water-vapour channels simulated from the analysis run 10–14 K too warm, while the background model is within 1 K of them. Assimilating the 6.9 µm channel corrects this (section 5c), but the forecast does not yet use the corrected humidity well.

• Clouds in the forecast do not decay as observed. On the 6 October growth case the satellite loses 61% of the analysis-time area colder than −20 °C within the hour, while the forecast keeps 75% of its own. A budget of the 5 October forecast above 450 hPa shows horizontal transport adding condensate where it should add almost none, and deposition adding more.

• About a third of the observed cloud at 60 minutes enters from outside the domain. On the 5 October case only about 10% of the domain is refilled by inflow within the hour, at upper-level winds near 8 m/s, yet that inflow area holds 37% of the observed cloud colder than −20 °C at 60 minutes, which no analysis of our domain can contain. The forecast takes the analysed state at its inflow edges, so a storm feeding an edge at the analysis time keeps feeding it; cloud still outside the domain at that time stays out of reach until the inflow comes from the background model's own forecast (8f).

• The rotational winds are held fixed through the forecast. The divergent part follows the vertical motion, and the steering flow a storm moves in stays as analysed. On the 5 October case this matters less than expected: the background model's own upper-level winds change by only about 1 m/s over the hour, a displacement of 1–2 km. Section 8f sets out the steps toward winds that evolve.

4a. How the prognostic processes came to be enabled

Section 2b lists the prognostic processes built into the model. The record below explains how each reached its present status in Table 1.1.

• Now on: sedimentation and autoconversion. Sedimentation had the strongest a priori case of anything on this page. The magnitude argument is first-order rather than a correction, since rain at 5 m/s falls about six model levels per 5-minute observation interval and roughly 9 km over 30 minutes. Its adjoint has been verified to work, with the gradient peaking at the source level that falls to the ground, exactly the vertical lever a reflectivity observation needs on the initial state. Both run in real time with the conserved pair. Precipitation evaporation remains untested.

• Tested and roughly neutral: accretion. At a defensible Kessler-scale rate it is statistically indistinguishable from being off on a deep-convection case.

• Now on: cloud condensation, through the conserved pair. Condensate is diagnosed from total water rather than controlled directly, which prevents a free condensate control from out-competing the humidity pathway, so humidity-driven condensation acts in the window and the forecast.

• Superseded: the adiabatic lifting source and subsidence evaporation. The lifting source improved the IR by 12–13% on a convective case and 16–33% on an orographic one with the species formulation. With the conserved pair, rising air condenses and sinking air evaporates through the saturation adjustment itself, so both are bypassed rather than added a second time.

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 caution from section 5 still applies: the horizontal semi-Lagrangian scheme is not monotone, so rates tuned against it may partly compensate a numerical error.

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.

5. Why the storm anvil overshoots

Forecast anvils spread faster than observed. On 28 July 2026 the observed ≥5 dBZ area grew 27% in 20 minutes while the model's grew 41%; on the decaying 5 October storm the forecast area colder than −20 °C reached more than three times the observed at 60 minutes in the configuration before 16 October (1.85 times now, section 5c). Three mechanisms contribute, each with its own remedy.

5a. Missing sinks

An anvil that spreads with no evaporation of cloud and no fall of precipitation keeps all its condensate at the leading edge. Evaporation through saturation adjustment, sedimentation and the fall of cloud ice now run in the forecast (Table 1.1). The numerical effect of section 5b comes first in the order of work: sinks tuned before it is fixed would be tuned to cancel a numerical error.

5b. Mass created by the advection scheme

With the minimizer disabled, so the forecast is pure advection of a conserved scalar, total condensate should stay constant. It grows instead:

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.

The growth is numerical. It sits in the interior (more than 40 grid points from any edge, 87% of the mass, +12.4% between 1800 and 1815) rather than at the boundary (+2.9%), and it is concentrated aloft where the anvil is: +82.6% above 200 hPa, +20.5% at 200–400 hPa, +8.8% at 400–600 hPa and +5.1% at 600–800 hPa. The semi-Lagrangian scheme's cubic interpolation dips below zero on the clear side of every cloud edge, and resetting those values to zero adds material each step: 0.85% of the upper-level condensate per one-minute step on the 5 October case, and 21% of that anvil's mass over 55 minutes of pure advection. A per-level mass fixer removes the growth in tests. A second effect remains with mass held fixed: repeated interpolation spreads cloud features one or two grid cells, as much as the 5th-order, flux-form, positive-definite scheme WRF uses, at a tenth of its cost. Tracking condensate as particles keeps those features intact and conserves mass exactly, and is the candidate if the spreading turns out to matter.

5c. An analysis that adjusts cloud through vertical motion

Humidity aloft. Visible-clear columns are dried until they render no cloud, and no observation measured humidity in clear air, so the analysis could dry the upper troposphere far too much: on 5 October the simulated 6.2, 6.9 and 7.3 µm water-vapour channels (RTTOV) ran 14.4, 10.8 and 6.3 K too warm in clear air while the background model was within 1 K, and the Grand Junction sounding gave ice relative humidity of 0.3–0.5 at 200–350 hPa against 0.00–0.09 analysed. The analysis then made the observed cloud with extra ascent, and the forecast kept it.

The water-vapour channel. The 6.9 µm channel is now part of the cost function, as a third channel of the same radiative transfer, with its vapour absorption varying with pressure and fitted to RTTOV 14 (0.3–0.8 K in clear air, peak sensitivity within 25 hPa of RTTOV's). It takes the clear-air error from about +11 K to +1.6 K. It acts at the window start only: fitted at both window times, the minimization met it with ascent over about 70% of the domain, and every constraint on large-scale vertical motion that stopped this also removed the real mountain-scale motion over the Front Range.

One temperature. The window physics, the infrared render and the forecast start from the temperature the analysed pair gives, together with the analysed vapour and condensate. Before, they started from the background temperature, so the liquid-water potential temperature they carried was not the analysed one. The forecast improves at most leads on five cases.

The compensating pair. The analysed vertical motion still carries two compensating parts in every case: broad descent of 0.08–0.14 m/s between 225 and 600 hPa, where the background is near zero, and small ascent cells (rms 1.2–1.9 Pa/s) between 450 and 700 hPa. The cells grow cloud and the descent damps it. Within the window the two nearly cancel, so the fit to the images is good; in the forecast each acts on its own. Started from the cells alone, the 5 October forecast reaches a 60-minute frequency bias of 7.0 for cloud colder than −20 °C; from the descent alone, 0.6; from both, 3.4.

The real-time configuration. The forecast starts from the first-guess vertical motion plus half the analysed increment, with its broad and small parts kept in proportion: a weight on the initial state, in the spirit of incremental analysis updates. With the water-vapour channel and the single temperature, this raises the forecast ETS by 0.067 over the configuration before 16 October, with 89% of 36 leads better across five cases. On 5 October the amount of cold cloud improves: the 60-minute bias falls from 3.65 to 1.85 and the ETS rises from 0.065 to 0.072 (Figure 5.1; the two scores are explained in section 8e); placement, measured by the frequency-matched ETS (Figure 5.2), improves slightly at 0 °C (by 0.02–0.04 from 10 to 45 minutes) and falls more at −20 °C, where the earlier configuration places the cold cores better at every forecast lead (by up to 0.056 at 35 minutes; 0.435 against 0.382 at 25). The gain on this case is in the amount of cold cloud, and the placement of the cold cores is the open part. Both growth cases improve (60-minute ETS 0.19 to 0.21 and 0.27 to 0.43), and so does Boulder at the analysis time (−20 °C ETS 0.73 to 0.83). The forecast's vertical motion then differs from the window's, and the growth case falls short of the observed cloud area after 25 minutes (bias 0.45 at 60 minutes). Weights of 0.25 and 0.75 trade ETS against amount, and 0.5 scores best overall.

Stopping the pair at its source. Two images five minutes apart observe vertical motion mainly where cloud stands on terrain. Allowing the full 1 m/s error inside the stationary-wave zone and 0.2 m/s outside it removes the descent aloft (5 October, 300 hPa mean +1.38 to +0.17 Pa/s), and the forecast can then take the whole increment: the 60-minute bias reaches 1.30. The same constraint lowers the Boulder wave cloud at the analysis time (−20 °C ETS 0.83 to 0.74), since the zone covers 1–12% of each tile, and it widens the growth case's shortfall. A looser error outside the zone, or a wider zone, comes next.

Coherent ascent, its outflow, and the domain edges (20 October). Three changes went into real time together. The analysed increment of vertical motion had lobes two to four levels thick, its vertical gradient changing sign between 41% of adjacent level pairs, because nothing tied one level to the next; it is now correlated in height over 1.5 km and, in cloudy columns, also takes the shape of the cloud layer, scaled by moist stability and condensate path. The horizontal wind then takes its divergent part from that ascent, screened at 20 km: an unscreened solve balanced the level-mean vertical motion with a domain-wide wind change of 5–10 m/s (median) and lost the cold cloud, while at 20 km the change is 0.5–1.4 m/s and a convective cell keeps 93% of its outflow at 4 km. Third, the forecast takes the analysed state at its inflow edges: on 5 October, 410 of the 421 core pixels colder than −40 °C lie within 5 km of the southern edge, and the analysed inflow keeps that edge supplied as observed. Against the 16 October configuration, over five cases and 126 lead and threshold pairs, the mean distance of the frequency bias from 1 falls from 0.240 to 0.232 (63% of pairs better) and the frequency-matched ETS rises from 0.534 to 0.552 (82% better). The largest gains are at Boulder (−20 °C frequency-matched ETS 0.58 to 0.70) and on the decaying storm of 5 October at −20 °C (distance from 1, 0.36 to 0.30; frequency-matched ETS 0.35 to 0.41). The 5 October cores colder than −40 °C score lower (frequency-matched ETS 0.13 to 0.07), and the calm growth case holds slightly too much cloud at 0 °C. Figures 5.1 and 5.2 show the configurations through 16 October.

Figure 5.1. Infrared frequency bias by lead, 5 October, at 0 and −20 °C: the amount of cloud through the hour (1 is the observed amount). Old deployed, THWT (whole w increment), new deployed (half), the outside-zone w constraint; HRRR, persistence.

Figure 5.2. Frequency-matched infrared ETS by lead for the runs of Figure 5.1: each forecast thresholded to the observed area, so it scores where the cloud is, apart from its amount. New deployed: slightly higher at 0 °C, lower at −20 °C.

6. The ice-phase saturation correction

Section 2a describes the correction: in cold air, relative humidity is judged against ice saturation rather than liquid. Its measurement, made with the earlier per-species formulation (section 3), also shows two things worth knowing when judging any change to this system.

6a. It acts through cloud fraction

The expected route was the humidity penalty. The larger effect came through the observation operator instead, because the same relative humidity feeds the convective cloud-fraction diagnosis: ice-aware humidity reads higher in cold air, so the diagnosed cloud fraction rises (Table 6.1). The state itself barely changes. Inside the analysed region the two configurations agree on condensate to within 0.3% (125,629 against 125,430 cloudy points), while their visible fits differ by 24%: the correction changes how a state renders.

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.

6b. How large the effect is in each domain

By the combined figure of merit, the correction improves the cloudiest 200×200 subset by 16.5% at the analysis time and leaves the full domain about 1% behind (Table 6.2). Two effects account for the difference, and both apply to any change:

• A cloudy subset overstates a change that scales with cloudiness. Subsets settle whether a change acts at all; full domains settle by how much.
• The minimizer absorbs a background-side gain. Without the correction the background is clearly worse (infrared 3.63 against 3.20 K), yet both converge to about the same analysis (2.51 against 2.53 K). A real improvement to the forward operator can show up almost entirely at the background.

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.

Over the forecast leads the full domain stays about 1% behind, and a daytime convective case came out about 0.7% ahead, so the combined score is indifferent to the change. It runs in real time because it's the physically correct formulation, it stops charging legitimate cirrus for existing, and it improves the infrared fit at the window start (from 4% worse to 7% better). It should also keep condensate near 125–175 hPa, where the air is close to ice saturation, and thin it at ice-subsaturated levels; the 28 July 2026 case, with real cirrus, is the planned test.

7. Wave cloud, cloud height and the tiled analysis

This section covers how the system handles mountain-wave cloud (7a–7c), how the analysis is split across the domain (7d), what the physics adds over moving the image (7e), how cloud height is checked (7f) and the rule that keeps hydrometeors above the ground (7g). Table 1.1 summarises the resulting real-time configuration, and section 8g collects the problems still open.

7a. Standing crests and moving cloud

Two regimes. Mid-level cloud over the mountains often has two parts moving at different speeds: wave crests locked to the terrain, and a cloud envelope that advects through them. A forecast that carries everything at one velocity reproduces the envelope and carries the crests away with it. Written as an equation, what the observations require is

     dq/dt + U . grad(q) = S(x)

The evidence. On the wave case of 2 September 2026 the forecast falls behind persistence over high terrain, by a margin that widens with lead, while over flat ground in the same domain, with the same model, winds and scheme, it stays ahead and its margin grows (Table 7.1). On an ordinary convective case the forecast leads persistence over high relief by a growing margin, so the effect belongs to terrain-locked cloud rather than to mountains in general.

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. Infrared RMS margin of the forecast over persistence (K) on observed-cloudy pixels, by terrain zone and lead; positive means the forecast is ahead.

The two components also separate by scale, measured from the observations alone. Band-passing the two window-time visible fields into octaves and finding, for each band, how much the correlation improves when one field is shifted along the wind, the 4–8 km band (the crest scale) gains nothing, −0.023, while every other band gains: +0.039 at 2–4 km, +0.039 at 8–16 km and +0.057 at 16–32 km. The crest band stands and everything larger and smaller advects through it. The figures subtract the same statistic computed across the wind, where advection can't contribute, to remove the bias of choosing the best of 21 trial shifts. A growth case (23 August 2026) shows no standing band at all, so a single standing scale is a property of the wave case rather than of the method.

What runs in real time. Cloud is carried by the analysed wind everywhere by default. Standing wave cloud is added only where the observations show it and at the heights where it can exist: in columns where the visible cloud moves less than 3 m/s between the two window images and the terrain within 20 km varies by at least 300 m, in a layer from 1 km above the local terrain up to the highest level consistent with the observed infrared brightness temperature (7c, 7f). Elsewhere cloud advects, and first-guess cloud the infrared rejects is removed in the same way everywhere.

7b. The conserved pair

The earlier formulation analysed each hydrometeor species directly, and it had two defects. Descent did not act as a sink: the correlation between vertical motion and cloud change reached only −0.14, and cloud grew on average even in sinking air. And condensate was placed without debiting the vapour that formed it: the small-scale total-water perturbation correlated with condensate at +0.32 to +0.39 through the wave-cloud layer, where a conserved formulation gives about zero.

The remedy follows cloud-resolving models such as SAM (Khairoutdinov and Randall, 2003): analyse conserved quantities and diagnose vapour, liquid and ice from them, so the inconsistency can't be represented. Total water and liquid-water potential temperature are the analysed pair; the forward model advects the pair and re-diagnoses temperature, vapour, liquid and ice at every step. Every operation has to stay differentiable, so a smooth saturation adjustment with an analytic derivative stands in for an iterative one. With the pair in place, descent evaporates the crest cloud, removing 84–91% of it in five minutes, and rising columns double their cloud: a vertical-only test with the model's own transport and saturation operators produces as much condensate as an exact parcel displacement.

The spatially fixed source S(x) is a vertical motion placed from the observed visible imagery. Its shape is the along-wind derivative of the observed cloud, zero-mean by construction so that ascent balances descent, and its amplitude comes from the saturation inversion.

7c. Holding wave cloud in place

The 22 September 2026 case (section 9) has lee-wave crests over Boulder County with a sky camera directly beneath them, so it tests height as well as placement.

Height of the seed. The LAPS cloud analysis placed the crest cloud near 450–475 hPa; the sky camera puts it near 600–700 hPa, and so does the infrared: a synthetic crest cloud laid on the observed visible pattern and passed through the operator's footprint reproduces the observed 1.4 K crest-to-trough contrast only between about 650 and 750 hPa. The wave-cloud seed is therefore confined to the layer described in 7a, and first-guess cloud the infrared rejects is removed before the minimization.

A moisture lid. Moistening is blocked above the tropopause unless a tall radar echo and a very cold observed cloud top indicate an overshooting storm. On this case stratospheric moisture would otherwise saturate into ice pockets and render the infrared 10 K too cold over northwest Grand County.

A balanced start. The observed crests are a steady balance, with air condensing as it enters a crest and evaporating as it leaves. Cloud simply added on the crests drifts 1–1.5 km downstream within the five-minute window and loses most of its contrast. In the wave zone the analysis therefore starts from a balanced state: the conserved temperature and total water are displaced along the streamlines of the seeded vertical motion, crests and troughs together, so that horizontal advection and lifting cancel. The drift over the window falls to about 0.5 km, one grid point, and the visible misfit falls by 74% at the window start and 39% at the analysis time.

Which columns. Every edge of visible cloud could seed a wave zone, so the zone needs two conditions, both measured from the observations. Between the two window images the Boulder crests move about 1–2 m/s while the cloud over northwest Grand County moves 6–10 m/s with the wind; a low deck over the flat southeast corner is just as still as the crests, so terrain supplies the second condition (100–150 m of relief within 20 km there, 600–1000 m under the Boulder crests). With both, the balanced start covers about three quarters of the Boulder cloud, an eighth of the northwest Grand County cloud and none of the southeast deck.

Result. Against the configuration it replaced, the Boulder visible correlation rose from 0.72 to 0.77 at the window start and from 0.73 to 0.80 at the analysis time, the wave-band correlation from 0.72 to 0.81 at the analysis time, and the full-domain figure of merit improved in every region at both times (2.56 to 2.34 and 3.41 to 2.79). Some smaller crests are still missing, and the cloud cover at the analysis time is under half of that observed.

7d. The tiled analysis and a step length per region

Tiles. The full 599×999 domain would need about 114 GB for the minimization, since it carries the record of every operation for the adjoint. The forecast doesn't, so only the analysis is split: into six overlapping tiles, with a 24-point overlap set by how far information travels in one five-minute window (about seven grid points of advection plus the five-point infrared footprint), and one forecast then runs over the whole stitched domain. Cross-fading the seams over six grid points brings every seam below the typical column-to-column variation.

Why tiles help. Several quantities the analysis would otherwise treat as domain-wide vary across the domain. The steering wind varies 38% across a 154×200 km box, and it sets both the amplitude of the placed vertical motion and the direction of its ascent-descent pattern; the largest gradient, the scale for the step, differs by a factor of 2.3 between the two halves of that box. Each tile has its own steering wind, and on two cases the tiled analysis, measured with eight tiles, beats a single tile on 15 of 16 measures at the two window times (Table 7.2).

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.2. 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.

A step length per region. Within a tile, the line search chooses a separate step for each 48-point (24 km) region, from a parabola through that region's observation cost, and accepts it only where the total cost falls. Choosing from the observation part matters: the background part rewards small steps everywhere, and a region's fit to the satellite is what the step is meant to improve. Larger tiles then lose nothing to a compromise step, so the domain runs in six tiles rather than eight, and both window times improve on all three cases (Table 7.3). On a night case the infrared error fell from 6.08 to 3.89 K at the window start and from 6.28 to 5.38 K at the analysis time.

Case	Lead	8 tiles, one step per tile	6 tiles, regional step
Boulder lee waves	0000	1.60	1.19
	0005	2.46	2.22
Moving high cloud	0000	2.03	1.53
	0005	3.66	3.50
Orographic	0000	3.31	2.92
	0005	5.42	5.19

Table 7.3. Figure of merit, (IR RMSE / 2 K)² + (VIS RMSE / 0.03)², over the full domain at the two window times; lower is better. The regional step with fewer, larger tiles improves both window times on all three cases.

7e. What the physics adds

The physics in Table 1.1 is moist thermodynamics on a given flow: condensation and evaporation through the conserved pair, latent heating, vertical transport, sedimentation and autoconversion, in the window and the forecast alike, with prescribed winds. A fair test of what it contributes needs three forecasts. Image advection moves the observed satellite image along its own measured motion, the standard no-physics nowcast. Physics off runs our analysis, minimization and three-dimensional winds but only advects between observation times. Physics on is the real-time configuration. Comparing the second and third isolates the physics; comparing the first and second isolates what the analysis and the three-dimensional winds add.

Figure 7.1. Visible correlation with GOES against lead; shading marks the window. Image advection starts from the observed 0005 image, so its line jumps at lead 10. Both model forecasts use the real-time configuration of 28 September 2026 (Table 1.1).

Table 7.4. Visible correlation at the window start and the analysis time (0000 / 0005), real-time configuration of 28 September 2026. The orographic case has no image ten minutes before the window, so image advection cannot be formed there.

In the window, where the system is meant to add value, the model beats image advection in visible correlation on both cases with a reference, and the physics adds to it on the two terrain-driven cases. On the orographic case it raises the analysis-time correlation from 0.804 to 0.868: condensation and evaporation following the terrain-forced ascent and descent. On the moving high cloud, with no terrain forcing, the physics costs a little visible correlation (0.009 and 0.004 at the two window times) while lowering the infrared error. By the figure of merit the physics improves the analysis on all three cases: at the analysis time from 3.49 to 2.19 on the Boulder case, 4.28 to 3.15 on the high cloud and 11.40 to 5.52 on the orographic case (Figure 7.3).

In the free forecast image advection keeps its lead on every measure. It starts from the observed image itself, while the model starts from its own analysis, a fit to that image with a correlation of 0.86 to 0.97. The physics keeps its lead over physics-off throughout on the Boulder case; on the other two it leads at ten minutes, is level at fifteen and falls behind at twenty.

Standing crests. The physics is what keeps wave crests in place. Tracking the analysed cloud itself, above ground only, the physics-off forecast carries it downwind at 8.3 m/s, the wind at 625–650 hPa; with physics on the pattern stays put, because air condenses where the standing lift raises it and evaporates where it sinks. A crest that brightens just east of the main Boulder crests shows the difference (Figure 7.2): the observed reflectance rises from 0.103 at 15:25 to 0.200 at 15:30 and 0.240 at 15:35. Image advection persists whichever image it starts from, and the physics-off crest dims, while the physics-on forecast follows the rise and reaches 0.220 at 15:35, the closest of the three. The physics-off panels look steadier than the drift implies because the simulated visible includes the fixed brightness of the ground.

Figure 7.2. Boulder visible at 15:25, 15:30 and 15:35, observed and in the three forecasts, with the crest that brightens after 15:25 circled. Image advection barely moves these standing crests, so it persists its starting image.

Figure 7.3. Full-domain figure of merit (top, lower is better), visible correlation (middle) and infrared RMSE (bottom) against lead, for image advection and the model with physics off and on (the real-time configuration of 28 September 2026).

7f. Cloud height

Moving cloud. Moving cloud belongs at the altitude whose wind carries it at its observed motion. Every real-time cycle checks this: for each cloudy column the background winds are searched for the level whose motion matches the observed cloud motion, and that level is compared with the level the visible channel sees in the analysis. In the cycles checked the two agree within 50 hPa, and the analysed level moves with the observed motion, within 3 m/s, in 55–80% of cloudy columns. On days of extensive fast-moving high cloud, 37–65% of cloudy columns have no background wind level that matches the observed motion; there the limit is the wind, or cloud that develops rather than moves.

What each channel sees. The infrared responds within the top visible optical depth of about 2; visible reflectance saturates slowly, following the two-stream albedo bτ/(1 + bτ), with the backscatter fraction b of 0.15 for liquid and 0.19 for ice, so an albedo of one half falls at τ of 6 to 7. For an optically thick high deck the visible-weighted level therefore sits below the top, while the top itself agrees with the infrared and with the wind-matching level. The infrared footprint is also larger, so a small cloud partly filling it reads warm, and therefore low; the operator's footprint convolution accounts for this.

Standing crests. Wave crests don't move with the wind at any level, so the infrared, corrected for the transparency of thin cloud, gives their height. With a visible optical depth near 2 and an infrared emissivity near 0.5–0.6, solving for the cloud temperature that gives the observed brightness temperature puts the tops about 2 km above the terrain: near 550 hPa over the Continental Divide and 625–650 hPa over the Boulder foothills. The analysed standing vertical motion peaks at 625–650 hPa (0.4 to 0.8 m/s), so cloud placed much higher has no lift to hold it: placed at the humidity maximum, it rose to 450–500 hPa, drifted downwind and the visible correlation fell from 0.86 to 0.75.

A moist halo. Wave cloud fitted as a thin saturated layer in air near 45% relative humidity mixes away within the window. Around the fitted wave cloud the analysis raises the humidity to 80%, still below the condensation threshold, within 1.5 km across and one level above and below, in stationary wave columns only and never into columns the visible sees as clear. Over the full domain this raises the Boulder wave-band correlation from 0.934 to 0.944 at the window start and from 0.788 to 0.851 at the analysis time, with the infrared error about the same. Over the Divide the share of cloudy columns still holding cloud at the analysis time rises from 30 to 64%, at 575 hPa and 1.6 km above the ground: most of the way to the height the infrared indicates. The remaining step is to raise the cloud at the window start, where the visible fit still sets its height.

Figure 7.4. ATOC sky camera, Boulder, at 15:25 and 15:30 UTC, and simulated all-sky images from the original real-time analysis, an eight-tile rerun and the 28 September 2026 configuration; r is each sky's correlation with the camera.

The sky camera shows both effects. The original real-time analysis put a false low deck across the southern and western sky; the current configuration removes it and correlates best with the camera at 15:25 (0.456, against 0.410 and 0.443). The wave cloud near the western horizon is present in all three simulations, lower and less complete than in the camera.

7g. Hydrometeors kept above ground

Over the mountains the lower pressure levels of the grid continue beneath the ground. Those points have no real atmosphere, so a rule following the LAPS cloud analysis keeps cloud above the terrain and precipitation down to the first pressure level below it. Underground water is kept as vapour, so total water is unchanged. The rule also keeps false low cloud from forming underground and being carried up into the lowest levels above the ground (Figure 7.5).

Figure 7.5. ATOC sky camera at 15:25 and 15:30 UTC, and simulated skies from the analysis without and with the below-ground rule. Without it, spurious low cloud lines the southern and western horizon at 15:25.

Monitoring. Every real-time cycle is scored once its verifying images arrive, against GOES and against image advection over the same leads, together with the height check of 7f. The results accumulate per domain, so the configuration is judged across many cycles rather than a handful of cases, alongside visual checks of the imagery.

8. Future work

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.

8a. Numerics: a conservative advection scheme

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.

8b. Physics still to enable

Sedimentation, autoconversion, condensation and evaporation, and vertical transport now run in real time (Table 1.1), and the conserved pair supplies what the lifting source was meant to. What remains:

• Precipitation evaporation (--precip-evap), the leading-edge sink on the precipitation side. It is built and off: on the September cases the model is already short of precipitation inside radar coverage, and the sink lowers the radar skill further.

• Rain and snow as controls (--precip-control). The radar term now acts at a weight comparable with the satellite terms, which raises the 20 dBZ equitable threat score at +5 minutes from 0.27 to 0.32 on a night case, but rain and snow still come from the background. Carrying them in the control vector is built and under test: it lowers the infrared error at the window start from 2.85 to 2.11 K, with the radar threat score about unchanged. The rain and snow gradients run about 25% below the true slope, so the reflectivity path still has a non-smooth piece to find.

• Ascent that sustains the observed precipitation (--radar-precip-w), in the manner of the LAPS cloud omega. It is built and off, awaiting a case with widespread moderate echo.

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. For convective cores this is the largest identified gap between the forecast and the observations; for wave and moving cloud it is cloud height (8g).

8c. The observation operator and cloud placement

Most changes measured through August traded IR against visible; the September changes to the first guess and the per-pixel sun angle improve both. 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.

Visible reflectance at low sun. With each lead rendered under its own sun (section 2a), two effects the operator leaves out show up in the observations near sunset. First, diffuse light from the sky: on 5 October, pixels where the observed cloud is moderately bright (reflectance 0.15–0.40) render 0.03 too dark at 5 minutes and 0.12 too dark at 60, as the sun drops toward 6.5°, while clear pixels hold steady. A skylight term would add diffuse illumination at the cloud top, larger for low cloud under more of the atmosphere's scattering. Second, cloud shadows: the observed clear ground shows them, the simulated doesn't. A shadow falls a distance h/tan(altitude) away from the sun, about 17 km for a 2 km cloud top at 6.5° and about 70 km for an 8 km top. The operator would attenuate the direct beam reaching the ground along the slant path through each layer toward the sun. The two share the diffuse part of the illumination, so they are best built and calibrated together. The downward shortwave irradiance of the LAPS derived products includes a twilight term, and serves as a check on the size of the diffuse light.

Two refinements to the infrared term were added in September 2026. The second of them runs in real time; the first is off by default. The first corrects the liquid-cloud absorption at 11 µm. The operator had taken the absorption efficiency as one at every droplet size, which fixes the ratio of infrared to visible optical depth at one half. A Mie calculation for water at 11 µm gives an absorption efficiency of 0.84 at the 8 µm effective-radius floor and 0.71 at 6 µm, so thin small-droplet cloud, such as orographic wave cloud, is more transparent in the infrared, relative to its visible brightness, than the operator assumed: by 1.2 times at 8 µm and 1.4 times at 6 µm. A fit to the Mie result, accurate to about 1% from 4 to 30 µm, now replaces the fixed value when enabled. The second weights the infrared misfit more tightly wherever the simulated brightness temperature is colder than observed. Pixels whose observed temperature lies within 10 K of the surface are otherwise treated as clear and given a looser error, and that classification also covers warm low cloud, where a too-cold simulation indicates cloud the satellite does not see.

A first guess from the background model's columns. Our first guess has its placement from the satellite but its vertical structure from the opaque-top assumption (section 1d). A model background has the opposite strength: physically consistent vertical structure, but features in the wrong places. We plan to combine them by column matching, in the manner of database retrievals such as GPROF: for each observed pixel, take the vertical profiles of the background columns nearby whose simulated infrared, visible, water-vapour and split-window signals best match the observation. The database is whatever the LAPS background ingest provides in the lga, so it serves the HRRR here and another model elsewhere; the lga carries the background's cloud liquid, ice, rain, snow and precipitating ice. Because each column carries all of its layers, low and high clouds moving in different directions need no motion estimate. Candidate columns would come from similar terrain, and where no cloudy column lies nearby the cloud-analysis first guess remains. A related approach aligns the background's fields to the satellite with displacements estimated layer by layer, holding terrain-locked cloud in place. Both will be sized offline first, against the cloud-top heights retrieved from the water-vapour and window channels together.

A longer window. Two images five minutes apart cannot show whether a storm is growing or decaying. Earlier images, which we archive, would show that trend, and in sheared flow they do more: as the upper cloud slides over the lower, successive images expose different parts of the low cloud, which the window's level-by-level advection can then assemble.

All 16 channels. The real-time cost function fits two of the 16 ABI channels, the 0.64 µm visible and the 10.3 µm window, with the 6.9 µm water-vapour band and the 12.3 µm split-window band in testing (section 5c). The remaining bands carry information the analysis currently infers indirectly: the near-infrared bands (0.86, 1.6 and 2.2 µm) hold particle size and phase, 7.3 µm the mid-level humidity, 8.4 µm the phase again, and 13.3 µm the cloud height. Two published routes use all 16 channels (section 10): a retrieval of cloud-top products assimilated as a humidity proxy (Otsuka et al. 2021), and a learned mapping to water-content profiles on fine vertical layers (C3DIR, White et al. 2026). The plan here is to take the channels to the model levels directly: a C3DIR-style mapping as an engine for the first-guess vertical structure, and the cost function extended band by band, each new operator calibrated against RTTOV as the water-vapour band was. The infrared bands come first, being straightforward; the near-infrared reflectances need size-dependent operators, and 3.9 µm, mixing reflected sunlight with thermal emission, comes last. To our knowledge, direct assimilation of all 16 channels has yet to be attempted, so this extension would open new ground.

8d. Minimization

The step length is now chosen region by region (section 7d), which let the domain run in fewer, larger tiles and improves both window times on all three test cases. The conserved-pair formulation can add rain and snow blocks (8b); the rain block needs a smaller starting step than the others to descend. Choosing the step from each region's observation fit may also allow fewer iterations, shortening the cycle. Hakim et al. (2026; section 10) minimise a similar differentiable 4D-Var with AdamW, an adaptive per-element step suited to a cost surface as far from quadratic as ours. Adam directions fed into our line search descended more slowly, since the two then compete for control of the step; the test still to run is Adam in its own regime, as the step control over a fixed iteration budget.

8e. Evaluation of results

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 every substantive change is compared against simple image advection and against the same model with its physics off (section 7e), while every real-time cycle is scored against both GOES and image advection as it verifies. Kurzrock et al. (2018; section 10) place the useful horizon of satellite cloud assimilation in hours, so as the window lengthens and the dynamics evolve (8f), the verification will extend from the present 60 minutes toward that range.

Amount and placement are scored separately, since a forecast can get one right and the other wrong. The frequency bias measures how much cloud the forecast holds at each lead, and so its growth and decay: the life cycle is judged on it, at −20 °C and at −40 °C for storm cores, together with the share of the analysis-time cold pixels that are still cold at each lead, model against observed. The frequency-matched ETS measures where the cloud is, independent of its amount. A change is deployed only when both improve over the 0–60 minute forecast: the bias moves closer to 1 and the frequency-matched ETS rises, across the evaluation cases. A bias corrected by thinning cloud in the wrong places fails. Figures 5.1 and 5.2 show the two axes for the same runs.

8f. Longer term: evolving dynamics in place of prescribed motion

This remains a short-range forecast 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, and it sits upstream of the advection routine, so improved dynamics modify the same wind field. Part of that path is now in place. The forecast carries cloud and moisture vertically, and takes its vertical motion from three diagnosed sources: terrain forcing, the standing lee-wave field placed by the analysis (7a), and convective ascent from a parcel buoyancy calculation capped at the infrared cloud top. The conserved θ_l–q_t pair gives the latent heating of condensation and evaporation consistently. Since 9 October the convective part of the vertical motion also evolves, in the window and the forecast alike: it is carried by the wind, relaxes toward the parcel ascent of the current state, and is projected onto mass continuity, in the spirit of the pressure step of an anelastic model such as SAM. Since 20 October the divergent part of the horizontal wind follows it.

What the established models have and this one lacks is momentum. WRF, the HRRR, SAM and ICON all carry the winds prognostically, so a storm's inflow and outflow respond to its heating and fade when it weakens, and rain-cooled downdrafts end it. Here the rotational winds are held at their analysis values for the whole forecast, and their divergent part follows the diagnosed vertical motion through continuity rather than a momentum equation; section 5 shows the cost on a decaying storm. Three steps lead toward evolving winds, each taken only as far as it pays:

1. Background winds that change in time. Interpolate between the background model's winds at the analysis time and its next hourly output, as a limited-area model takes its boundaries from a parent, and take the inflow at the domain edges from the same hourly output. The forecast already takes the analysed state at its inflow edges (20 October); the time-varying winds and inflow come next, first for the analyses at the top of the hour, whose next hourly output is written early enough to archive. The background is read through the LAPS background file, so the same path serves the HRRR now and another regional model, or a nest of our own, later.
2. A divergent wind that follows the vertical motion (in real time since 20 October). The analysis wind keeps its rotational, steering part, and its divergent part comes from the current vertical motion through continuity, screened to about 20 km so that the storm-scale circulation responds while the larger-scale divergence stays the background's. Convergence then weakens as the convection does instead of holding it up (section 5c).
3. A prognostic vertical velocity driven by buoyancy, condensate loading and evaporative cooling, with the projection acting as the pressure solve. This is the first step at which a storm's decay would follow from its own physics.

Advection then 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.

8g. Open problems

• Wave-cloud height and survival. Emissivity-corrected infrared puts the Boulder wave tops about 2 km above the terrain, 50 to 75 hPa above the analysis, and most so over the Continental Divide (7f). The moist halo now in real time keeps most of that cloud through the window and lifts it to about 1.6 km above the ground (7f); raising the cloud needs its standing lift raised with it, since placing cloud above the lift made it drift and evaporate.

• Moving cloud: height and wind. In real time the analysed height mostly agrees with the level whose wind matches the observed motion (7f). What remains is the depth of thick decks, whose lower part moves with slower winds, and, on days of fast, extensive high cloud, columns where no level of the background wind matches the observed motion at all.

• Crest drift. The balanced start holds most crests within a grid point over the window (7c), but smaller crests and the cloud cover at the analysis time remain short of the observations.

• Rain and snow controls (8b), and a monotone horizontal advection scheme (8a).

9. Case studies

• 28 July 2026, 1800 UTC (archive key 262091800, 1755 reftime) — the richest of the four: cirrus above 200 hPa, deep ice cloud, and live precipitation including graupel, with an active radar cost. Archived complete with its operational forecast output, so it can be re-run. All section 3b, 4 and 5 numbers come from here.

• 26 July 2026, 1530 UTC (archive key 262071530) — mid-level stratiform capped by 63%-ice upper cloud, precipitation-free; the IR-limited complement.

• 25 July 2026, 1800 UTC (archive key 262061800) — daytime convective, visible-limited; source of the CA-initialized comparison in section 3c.

• 23 August 2026, 1600 UTC (archive key 262351600) — thin cirrus over the Weld County plains with a large height error: observed drift 12.2 m/s against an analysis that had the cloud far too low. The case the height-inversion control of section 2a was developed and tuned on.

• 2 September 2026, 1400 UTC (archive key 262451400) — mid-level mountain wave cloud, partly anchored by terrain and modulated by standing waves. Archived for the forecast dynamics test of section 7a.

• 22 September 2026, 1530 UTC (archive key 262651530, 1525 reftime) — lee-wave crests over Boulder County and a false low deck in the southeast corner, with the ATOC and NREL sky cameras and their simulated all-sky images archived for 1400–1615 UTC. The case section 7c is based on, and the one on which the wave-cloud configuration now running in real time was chosen. Also the source of the sky-camera comparison in section 7f.

• 21 July 2026, 1830 UTC (archive key 262021830) — orographic cloud over the mountains with a clear plains contrast; the terrain case of section 7e.

• 24 September 2026, 2325 UTC (archive key 262672325) — high cloud moving steadily across the domain with no terrain forcing; the case on which physics costs a little and the visible cloud sits furthest below its wind-matching level (section 7f).

• 5 October 2026, 2300 UTC (archive key 262782300, 2255 reftime) — a late-afternoon storm complex over the eastern plains near Akron whose cores were already weakening; the anvil grew for about 20 minutes, then shrank. The case of section 5c, with 60-minute forecasts and the HRRR's own forecast archived, and the GOES-19 water-vapour imagery (bands 8–10) and the 00 UTC Grand Junction sounding kept with them under extra/262782300.

• 25 July 2026, 0050 UTC (archive key 262060050) — mature convection with strong 3-D radar; source of the steering-wind measurements in section 3a.

10. Related work

Kurzrock et al. (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 7a addresses.

Satoh et al. (2019). 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 forecast 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 7a.

Khairoutdinov and Randall (2003), the System for Atmospheric Modeling (SAM). The prior art for the change of variable described in section 7b: 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 et al. (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 assimilation and forecast fitting satellite radiances over five minutes — so what transfers is the treatment of the optimisation, not the configuration.

Otsuka et al. (2021). The first assimilation of cloud products derived from all 16 channels of a geostationary imager. The Himawari-8 optimal cloud analysis (an optimal-estimation retrieval with RTTOV, after Hayashi and Watts 2016) uses every channel of the Advanced Himawari Imager to retrieve cloud phase, top pressure, optical thickness and effective radius, with the visible, near-infrared and window channels each constraining one of these. The authors assimilate the result in JMA's mesoscale 4D-Var (JNoVA) as a pseudo-relative humidity: a cloud layer from the retrieved top down to a base estimated from the liquid water path at a fixed 0.3 g m^-3, set to about 90% RH, and a climatological clear-sky value where the first guess is cloudy but the satellite sees none. Over two of three heavy-rain cases the 3-h rainfall threat scores improve, with the humidity increments driving wind changes through the 4D-Var tracer effect where low cloud traces inflow from the ocean. Three of their findings bear directly on this system. Their analysis kept the hydrometeors outside the control vector, since the retrieved water paths differed widely from the model's; here the condensate is analysed through the conserved pair and fitted to the radiances themselves. The retrieval describes the uppermost cloud layer, so low cloud beneath high cloud remains hidden to it; section 8c describes how successive images in a longer window can expose that cloud. And the vertical extent of the cloud is where the most is to be gained: the estimated base correlates at 0.30 with ceilometers, with an 890 m high bias, and the authors point to the statistical cloud-base method of Noh et al. (2017), used in the NOAA JPSS cloud products, and to variables closer to the cloud itself as the way forward. That is the direction taken here: to use the full channel set to map the satellite observations directly to the hydrometeors and water vapour at each model level, carrying the vertical structure that the channels hold into the analysis itself.

White et al. (2026). A retrieval that draws on all 16 ABI channels (0.47–13.3 µm, pooled to a common 2 km grid), with the equivalent full channel sets of SEVIRI and AHI handled by the same network through a stem layer per instrument. Trained against EarthCARE cloud profiles matched along the imager's slant line of sight, and given a GFS short-range forecast of temperature and humidity, it returns vertical profiles of ice, cloud-liquid and rain water content on 0.25 km layers, together with cloud base, cloud top and multi-layer detection. The present system fits the 0.64 µm visible and 10.3 µm infrared bands, with the 6.9 µm water-vapour band (section 5c) and the 12.3 µm split-window band added in testing, each through its own forward operator inside the cost function. C3DIR shows how much vertical structure the full channel set carries, and its profiles are a natural candidate for the column-matching first guess of section 8c, where the vertical distribution of hydrometeors matters most: thin cloud aloft and cloud layered beneath it. Taken together, these two studies show the full channel set feeding a retrieval, and retrieved cloud products feeding an assimilation. Direct assimilation of the radiances themselves has so far covered the 0.6 µm visible band (operational at DWD since 2023) and the water-vapour and window infrared bands; to our knowledge, direct assimilation of all 16 channels has yet to be attempted. The present cost function already fits visible and several infrared bands together, making it a natural place to take that step.

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.

11. References

• 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).
• Hayashi, M. and P. Watts, “Two-layer cloud retrieval using visible to infrared bands of Himawari-8”, EUMETSAT Meteorological Satellite Conference, Darmstadt, Germany (2016).
• Heymsfield, A. J. and L. J. Donner, “A Scheme for Parameterizing Ice-Cloud Water Content in General Circulation Models”, Journal of the Atmospheric Sciences 47(15), 1865–1877 (1990).
• Kessler, E., “On the Distribution and Continuity of Water Substance in Atmospheric Circulations”, Meteorological Monographs 10(32), American Meteorological Society, Boston, 84 pp. (1969), doi:10.1007/978-1-935704-36-2_1.
• 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).
• 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).
• Noh, Y.-J., J. M. Forsythe, S. D. Miller, C. J. Seaman, Y. Li, A. K. Heidinger, D. T. Lindsey, M. A. Rogers and P. T. Partain, “Cloud-Base Height Estimation from VIIRS. Part II: A Statistical Algorithm Based on A-Train Satellite Data”, Journal of Atmospheric and Oceanic Technology 34(3), 585–598 (2017), doi:10.1175/JTECH-D-16-0110.1.
• Otsuka, M., H. Seko, M. Hayashi and K. Koizumi, “Data Validation and Mesoscale Assimilation of Himawari-8 Optimal Cloud Analysis Products”, Journal of Atmospheric and Oceanic Technology 38(2), 223–242 (2021), doi:10.1175/JTECH-D-20-0015.1.
• 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.
• White, C. H., Y.-J. Noh, J. M. Haynes and I. Ebert-Uphoff, “C3DIR: A Deep Learning 3-Dimensional Cloud Property Retrieval Scheme for Passive Satellite Imagers”, arXiv:2607.16929 (2026), doi:10.48550/arXiv.2607.16929.