Integrate refractivity with shape-preserving quadrature and use virtual temperature

[jlmaurer]

Two independent corrections to how zenith delay is computed from a weather model, plus the fixture regeneration they require.

1. Hydrostatic refractivity from virtual temperature

_get_hydro_refractivity used k1 * P / T. Every model in the package sets k2 = 0.233 K/Pa, which is the compensated coefficient k2' = k2 - k1 * Rd/Rv (Bevis/Davis: 0.716 - 0.776 x 0.6221 = 0.2333). That coefficient pairs with k1 * Rd * rho = k1 * P / Tv, not with k1 * P / T; pairing it with P / T double-counts part of the water-vapour contribution.

The fix uses the exact identity k1 * P / Tv = k1 * (P - (1 - Rd/Rv) * e) / T.

This lowers the total delay by roughly 2% of the zenith wet delay — a few mm in typical conditions — so it is not merely a re-split between the hydrostatic and wet components. Absolute delays change for every weather model.

2. Refractivity integrated by shape-preserving quadrature

Every quadrature rule is exact integration of some interpolant. The trapezoid rule integrates a piecewise-linear reconstruction, and refractivity is convex in height, so it systematically overestimates the delay. On the coarse 32-level pressure-level z-grid that bias reaches ~9 mm of zenith delay.

_getZTD now integrates a shape-preserving piecewise cubic Hermite (PCHIP) reconstruction instead, evaluating its antiderivative analytically. PCHIP rather than a natural cubic spline because it introduces no new extrema and cannot overshoot — a natural spline undershoots into negative wet refractivity on the majority of grid columns, which is unphysical.

Validation

ERA-5 across three climate regimes, 596 grid columns, comparing the 32-level pressure-level grid against the 145-level native grid:

quadrature	mean abs. ZTD error	max
trapezoid on 32 levels	8.90 mm	18.47 mm
PCHIP on 32 levels	2.09 mm	13.23 mm
either rule on native 145 levels	0.32 mm	0.39 mm

The control that matters: on the native grid the two rules agree to 0.32 mm, so the reference is essentially rule-independent rather than flattering PCHIP.

The improvement is climate-dependent, and the wet tropics are the weak case:

	trapezoid	PCHIP	gain
Alaska (arctic, dry)	8.78 mm	0.73 mm	12x
Mexico (subtropical)	11.86 mm	3.00 mm	4x
Brazil (equatorial, wet)	6.76 mm	4.28 mm	1.6x

In humid equatorial columns PCHIP still leaves ~4 mm (max 13 mm): 32 levels cannot resolve the near-surface water-vapour profile, and no quadrature rule recovers data that is not there. Brazil is also the only case where the trapezoid bias is not one-signed (87% of columns overestimate, versus 100% elsewhere).

Native model-level runs are effectively unaffected, so this matters specifically for pressure-level (pl) configurations.

Memory

Building one PCHIP interpolator over a whole cube costs ~20x the cube in float64 temporaries — several GB for a CONUS-sized grid. cumulative_integral_from_top therefore processes columns in chunks against a fixed budget. Peak is now 2x cube + ~90 MB rather than ~20x cube (measured: a 15.4 MB float32 cube peaks at 129 MB instead of ~320 MB; extrapolated to a CONUS HRRR grid, ~850 MB instead of ~7.5 GB). Columns are independent, so the result is bit-identical for any chunk size — pinned by a test that compares one-column-per-chunk against a single chunk with array_equal.

3. Regenerated test fixtures and golden values

The processed weather-model fixtures store precomputed wet, hydro, wet_total and hydro_total, and the tests that assert delay values read those stored cubes rather than recomputing them. They therefore encoded the old physics, and would have kept passing unchanged — silently asserting superseded numbers.

Eight fixtures were regenerated by recomputing the derived cubes from their own stored t/p/e. Every step that produces t/p/e is untouched by this PR, so this reproduces exactly what the current pipeline would write for the same raw input, without re-downloading.

The reconstruction was validated before anything was overwritten: recomputing wet (whose formula is unchanged) reproduced the stored values bit-exactly, as did recomputing hydro under the old k1*P/T formula — confirming both that the stored t/p/e are the true inputs and that the fixtures were generated with the old physics.

Updated golden values:

test	old	new	change
test_HRRR_ztd (hydro)	2.2622863	2.2611293	-1.16 mm
test_slant_proj	2.333865144	2.3324788251164725	-1.39 mm
test_ray_tracing	2.97711681	2.9756965061161926	-1.42 mm
test_gnss_intersect	2.34514	2.343549	-1.59 mm

All four fixtures here are native model-level grids, so the shift is dominated by the virtual-temperature correction rather than the quadrature change. It scales with water vapour as expected — the surface hydrostatic change ranges from -1.2 mm over dry HRRR CONUS to -4.95 mm over equatorial Brazil.

Not regenerated: the four 33.7 MB HRRR inputs under test/gunw_azimuth_test_data/weather_files/. They back test_azimuth_timing_interp_against_center_time_interp, which compares azimuth-time against center-time interpolation and is self-consistent either way, so no golden value depends on them. Regenerating would add ~135 MB of binary churn to this PR, and #815 proposes untracking them entirely.

4. Unrelated CI fix: retired HRRR mirror

test_hrrr_availability_check_using_gunw_ids was failing in CI with a socket timeout. This turned out to be unrelated to the changes above and pre-existing.

Herbie tries each HRRR source it knows until one has the file. The University of Utah Pando archive is being wound down in favour of NOAA's Registry of Open Data on AWS, and pando-rgw01/pando-rgw02.chpc.utah.edu now resolve in DNS while accepting no connections. Verified directly:

source	host	result
aws	noaa-hrrr-bdp-pds.s3.amazonaws.com	TCP OK, 0.07 s
nomads	nomads.ncep.noaa.gov	TCP OK, 0.06 s
google	storage.googleapis.com	TCP OK, 0.08 s
azure	noaahrrr.blob.core.windows.net	TCP OK, 1.08 s
pando	pando-rgw01.chpc.utah.edu	TCP timeout
pando2	pando-rgw02.chpc.utah.edu	TCP timeout

A dead mirror costs nothing when the file exists — but check_hrrr_dataset_availability asks exactly the question that walks the entire list. For an unavailable time Herbie reached both dead hosts and blocked on each until TCP timeout:

🤝🏻⛔ Bad handshake with pando? Am I able to move on?
ConnectTimeout HTTPSConnectionPool(host='pando-rgw01.chpc.utah.edu', port=443):
Max retries exceeded with url: /hrrr/nat/20160809/hrrr.t16z.wrfnatf00.grib2

That is why only the negative assertion failed, and why it was intermittent: the positive case (2023, available) stops at AWS immediately.

Measured on the exact case the test uses: 151.9 s and an exception before, 1.1 s and the correct False after. The previously failing test now passes in 25.9 s.

The retired entries are filtered out of the source list the model's own Herbie template declares, rather than hard-coding the survivors: new mirrors are picked up automatically, and each model keeps Herbie's ordering — hrrr tries AWS first while hrrrak tries NOMADS first, so a single hard-coded list would silently reorder one of them.

Note this is not fixed by upgrading Herbie — current Herbie main still lists both Pando entries, and no upstream issue tracks it. Worth reporting there separately, since every Herbie HRRR user hits this on any unavailable date.

Known limitation: the ray-traced path

The ray-tracing path still integrates the trapezoid reconstruction, so zenith and ray-traced outputs will now differ by roughly 9 mm on pressure-level models. This is deliberate — see #817 for the analysis and why it needs its own PR.

Tests

Full suite (-m "not long"): all pass except test_hrrr_availability_check_using_gunw_ids, which hits live HRRR servers and fails with a socket timeout on this network, unrelated to these changes.