Ray-traced delays still use trapezoid quadrature, so they disagree with ZTD by ~9 mm on pressure-level models

[jlmaurer]

Summary

After #816, the zenith and ray-traced delay paths integrate refractivity with different quadrature rules, so the two outputs will disagree by roughly 9 mm of zenith delay on pressure-level weather models.

That PR replaced the trapezoid rule with exact quadrature of a shape-preserving (PCHIP) reconstruction in WeatherModel._getZTD, because refractivity is convex in height and the trapezoid rule therefore overestimates the integral on the coarse 32-level pressure-level z-grid. The ray-tracing path was deliberately left alone, and this issue records why that is not yet consistent.

Where the inconsistency lives

The ray-traced path integrates along the ray by finely subdividing each segment and applying the trapezoid rule:

https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/delay.py#L313

# Trapezoidal integration with scaling
wt = 0.5 if findex in [0, fracs.size - 1] else 1.0
wt *= ray_lengths[zz] * 1.0e-6 / (nparts - 1.0)

The values it samples come from a RegularGridInterpolator over the refractivity cube:

https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/delayFcns.py#L55

ifWet = Interpolator((ys_wm, xs_wm, zs_wm), wet, fill_value=np.nan, bounds_error=False)

which is linear in z by default. Subdividing finely and applying the trapezoid rule to a piecewise-linear-in-z field converges to exactly the piecewise-linear integral — that is, to the same answer the old trapezoid code gave. Increasing nparts does not help, because the bias is in the reconstruction between model levels, not in the step size along the ray.

So the fine sub-sampling makes the ray integration numerically accurate with respect to a reconstruction that is itself the biased one.

Consequence

For the same weather model and target, on a pressure-level model:

• raider.py ZTD output (and the std output projected via cos(inc)), which reads wet_total / hydro_total, now carries roughly 0.3 mm of quadrature error;
• ray_trace: true output, which reads the pointwise wet / hydro cubes, still carries roughly 9 mm.

Measured on real ERA-5 data across three climate regimes (596 columns; equatorial Brazil, subtropical Mexico, arctic Alaska), 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 trapezoid bias is one-signed (an overestimate) in 100% of columns in the subtropical and arctic cases and 87% in the wet equatorial case, so it does not average away.

Native model-level runs are largely unaffected either way — at 145 levels the rules agree to 0.32 mm — so this matters specifically for pressure-level (pl) configurations.

Why this is not a quick fix

Making the ray path consistent is not just swapping a rule. The ray integration samples refractivity at arbitrary points along a slanted path, so there is no single vertical axis to integrate a 1-D reconstruction along; the fix has to decide how a shape-preserving vertical reconstruction composes with horizontal interpolation and with the ray geometry. Plausible directions, none obviously correct without testing:

1. Build the interpolator over a PCHIP-in-z reconstruction (for example, resample the cube onto a finer z-grid via PCHIP before constructing the RegularGridInterpolator), so the linear-in-z sampling approximates the shape-preserving profile. Cheap to try; costs memory.
2. Integrate analytically per model-level crossing along the ray rather than by uniform subdivision.
3. Reformulate the ray integral in terms of the already-integrated wet_total / hydro_total where geometry permits.

More broadly, the ray-tracing path has not been fully vetted across the package, and it interacts with other open questions — notably the height-datum inconsistencies, where the zenith and ray paths currently interpret the same z axis under different datums. It deserves its own PR with dedicated validation rather than being folded into a quadrature change.

Suggested scope for a fix

• Decide the reconstruction the ray path should integrate, and make it explicit rather than an emergent property of the interpolator's default method.
• Add a test that pins zenith and ray-traced delays against each other for a vertical look direction, where they must agree to within the quadrature tolerance. No such test exists today, which is why this drift is invisible.
• Validate on a pressure-level model, where the effect is largest; a model-level-only check will not show it.

Related

• Integrate refractivity with shape-preserving quadrature and use virtual temperature #816 introduced the zenith-path change and deliberately left the ray path alone.