Neural flux stencils and oracle-assisted safeguards for soil-temperature reconstruction: a controlled negative result

Keywords: physics-informed neural networks; neural stencil emulation; operator error; soil temperature; moisture-dependent diffusion; reconstruction.

1 Introduction

A PINN penalises a solution field with a governing-equation residual [1, 2]. Replace one exact local term in that residual with a learned stencil and there is now a second approximation inside the penalty, so a low residual may mean agreement with the emulator and not with the equation. This paper asks when that extra approximation helps, harms or merely reproduces an exact numerical derivative, in an observed-temperature reconstruction task where the answer can be checked against a known oracle.

The reference flux map is known, differentiable and inexpensive, which makes it an oracle against which value and gradient error can be measured; no acceleration is claimed and no soil physics is discovered. An earlier pilot in the archive tested a learned replacement for an exactly known linear stencil, found no consistent benefit and is preserved as development history. The present study fixes a new protocol before scoring new windows, introduces a non-linear constitutive map, audits gradients on the states the PINN visits, adds a numerical reconstruction with the same information, and excludes the pilot’s January 2024 targets.

2 Related work and scope

Neural stencil emulation learns a discrete evolution term from simulator evaluations; designed stencil states reduce redundant training data [3], and learned discretisations [4] and STENCIL-NET [5] precede this work. Numerical spatial derivatives inside PINNs are established: CAN-PINN [6], meshless discretely trained PINNs [7], FE-PINNs [8], PhyGeoNet [9] and an auxiliary FD residual-gradient regulariser [10], which keeps the AD residual and adds a separate term and is therefore distinct from the replacement tested here. Because of these precedents the exact-flux FD arm is the essential control, since any effect of coupling neighbouring points must not be attributed to learning the operator. Moisture-dependent thermal properties and coupled hydrothermal PINNs are prior work [11]; the multiplier used here takes its form from a logistic conductivity family [12] and is a dimensionless effective-diffusivity assumption with an analyst-chosen shape, with no claim to be a validated law. Optimisation pathologies [13, 14, 15, 16], gradient-enhanced residuals [17] and generalisation theory under stability assumptions [18] are background; none establishes physical accuracy from a low training loss alone.

The contribution is a controlled observational evaluation, a pointwise operator-budget experiment and retained failure evidence. The error inequalities below are elementary decompositions that say what is measured.

3 Data and protocol fixed before scoring

3.1 Measurements and permitted information

Inputs are hourly USCRN files [19] for Newton (Georgia), Asheville (North Carolina) and Blackville (South Carolina); Crossville (Tennessee) was planned and supplied no eligible window. Each window has 96 hourly records at nominal depths 5, 10, 20 and 50 cm. Temperatures at 5 and 50 cm are supplied throughout; 10 cm every six hours; the first four-depth profile including 20 cm initialises every procedure; all later 20 cm temperatures are withheld and used only for scoring; moisture at all four depths is available throughout. What is simulated is a virtual withheld channel with auxiliary moisture rather than a sensor-failure incident, and full-window conditioning makes it retrospective reconstruction. Hourly averages are treated as discrete observations with piecewise-linear interpolation; no averaging operator or sensor-error model is fitted.

3.2 Selection, amendments and calibration

The plan crossed three stations with four 2024 start dates (20 January, 15 April, 15 August, 15 October); the first complete 96-hour block starting within seven days of each date is eligible, with temperatures in $[-50, 80]$ °C, moisture in $[0, 0.75]$ m³ m⁻³ and consecutive timestamps. Availability checks found no eligible Crossville block and no complete February 2023 calibration week; before any calibration or scoring an amendment added Blackville and selected the first available calibration candidate from 1 February, 1 March and 1 April (April was eligible). Unavailable blocks are recorded and never substituted. The final register has 16 planned combinations and ten eligible blocks at three stations, with winter coverage only at Asheville, so the design is far from a balanced factorial.

The April and July 2023 source weeks at Newton fit one global reference diffusivity $\kappa_0$ by minimising forward-solver 20 cm error, bounded in $[10^{-8}, 2\times10^{-6}]$ m² s⁻¹: $\kappa_0 = 5.18078\times10^{-7}$ for the variable-coefficient model and $2.39133\times10^{-7}$ for an independently calibrated constant-diffusion ablation. No 2024 target enters either. The ridge baseline is calibrated on the same source weeks. The protocol fixes architectures, seeds, caps, metrics, teacher support, gate thresholds, data rules and estimands; the scoring script refuses to run if the configuration and core-source hashes differ from the recorded control configuration. This is a local computational freeze with recorded history; public preregistration would be the stronger version.

3.3 Outcomes and inferential units

Twenty-centimetre RMSE over hours 1 to 95 is the primary score, computed per seed, averaged over three seeds within a block, then averaged equally over blocks. Seeds measure optimisation variation; hours and seeds are not station replicates. The primary contrast is neural-flux minus exact-flux FD; gated minus neural and neural minus AD are secondary; gated minus FD and the constant-diffusion ablation are exploratory. A station-cluster bootstrap resamples the three station labels with replacement (10,000 draws, seed 2041) for percentile intervals on paired block-mean differences; leave-one-station-out means and block win counts are also kept. With three clusters there are only ten distinct resamples, so the intervals are coarse descriptive devices. Three clusters and uneven seasons preclude population-level significance claims, and no result-based tuning, rerun or block replacement is permitted.

4 Operator, models and safeguards

4.1 Assumed equation

The effective one-way heat equation is $$\partial_t T = \kappa_0\,\partial_z\!\left[g(\theta)\,\partial_z T\right], \qquad 0.05 \le z \le \SI{0.50}{m}, \label{eq:heat}$$(1) with $\theta$ the observed volumetric water content and $$g(\theta) = 0.25 + 1.75\,\sigma\!\left(\frac{\theta - 0.20}{0.04}\right), \qquad \sigma(a) = \frac{1}{1 + e^{-a}}.$$(2) Non-linear in moisture, the map leaves the equation linear in $T$ given moisture. A constant effective heat capacity is absorbed into $\kappa_0$; advection, latent heat, freeze-thaw and texture are omitted. With $\xi = (z - 0.05)/L$, $L = \SI{0.45}{m}$, $s = t/86400$ and $u = (T - \mu)/5$ for $\mu$ the mean initial profile, $\partial_s u = \alpha\,\partial_\xi[g(\theta)\partial_\xi u]$ with $\alpha = 86400\kappa_0/L^2$. The continuum AD operator carries both $g(\theta)u_{\xi\xi}$ and $g'(\theta)\theta_\xi u_\xi$.

4.2 Matched fields and fitting

All five main neural arms share a three-hidden-layer, width-32 tanh field with inputs rescaled time, depth and daily sine-cosine features, and boundary lifting $$u_\psi(s, \xi) = (1 - \xi)\,b_0(s) + \xi\,b_1(s) + \xi(1 - \xi)\,N_\psi(s, \xi),$$(3) so boundary values are hard. The initial profile is linear interpolation of the four initial temperatures at 19 depths. The objective averages squared errors for the six-hourly 10 cm data, the 19-point initial profile and the physics residual with unit weights; the data-only arm drops the physics term. Each fit gets 800 Adam steps at learning rate 0.002 and an L-BFGS cap of 100 iterations on a seeded 512-point collocation set over $h \le \xi \le 1 - h$, $h = 0.05$, refreshed every 50 Adam steps. Seeds 11, 29 and 47 index paired teacher and field initialisations, so neural-flux seed variation mixes teacher approximation with field optimisation. The constant-diffusion ablation uses seed 11 only and is compared with the AD arm’s seed-11 fit. All runs use float64 on one CPU thread.

The AD arm uses the continuum derivative, and the FD arm exact face fluxes: with $d_\pm$ the forward and backward differences of $u$ over $h$ and $m_\pm$ the corresponding face-mean moistures, the face numerator is $q(d, m) = g(m)\,d$ and the spatial term is $[q(d_+, m_+) - q(d_-, m_-)]/h^2$. The neural-flux arm substitutes a frozen learned numerator $\hat q_\phi$; temporal derivatives stay automatic in every physics arm.

4.3 Teacher and gate

A two-hidden-layer, width-24 tanh network learns $q(d, m)$ from 4,096 uniform pairs in $d \in [-0.25, 0.25]$, $m \in [0.05, 0.35]$, with the odd construction $\hat q(d, m) = [N(d, m) - N(-d, m)]/2$ giving exact sign antisymmetry and zero flux at zero difference. Each field seed has its own teacher (2,000 Adam steps, L-BFGS cap 150), which must pass frozen value and Jacobian RMS thresholds on a disjoint 4,096-pair draw before scoring; accepted weights are frozen with input gradients enabled.

The gated arm evaluates the exact numerator online and accepts a face only if it is inside teacher support, $|\hat q - q| \le \epsilon_h$ and the componentwise Jacobian error is at most 0.02 in normalised units, where $\epsilon_h = b h^2 / (2\alpha)$ with $b = 0.02$ normalised temperature per day. Rejected faces use $q$; the branch mask is detached from the gradient. With $e_\pm$ the selected-minus-exact face errors, the triangle inequality gives at every evaluated stencil $$\alpha\left|\hat D^{\text{gate}}_h u - D_h u\right| \le \alpha\,\frac{|e_+| + |e_-|}{h^2} \le b, \label{eq:budget}$$(4) a bound of $5b = \SI{0.10}{\celsius}$/day in physical units on the evaluated spatial-residual perturbation relative to the exact discrete operator on the same field. It bounds nothing else. Truncation error, gradient damage elsewhere, solution error and the mismatch between Equation 1 and the soil all lie outside it. The gate is an oracle-assisted diagnostic; deployment would need an oracle-free version.

4.4 Numerical and statistical comparators

A forward finite-volume solver uses 37 nodes, five-minute backward-Euler steps, exact coefficients and observed boundaries, and does not assimilate the later 10 cm data. A matched-information FD reconstruction on 19 nodes and hourly steps solves a sparse least-squares system with backward-time flux residuals, the same six-hourly 10 cm data and the same initial profile, each row group scaled by the square root of its size to match the averaged neural loss, hard boundaries, LSQR damping $10^{-8}$, tolerance $10^{-10}$ and a 20,000-iteration cap. Ridge regression uses contemporaneous 5 and 50 cm temperatures, interpolated 10 cm readings and four moisture channels, standardised on 2023 data, penalty 0.01. Initial 20 cm persistence is the last baseline.

5 Verification before scoring

The control chain checks analytic continuum derivatives against an independently written manufactured source, smooth-field FD refinement, a sine-decay solution, a stationary affine sparse reconstruction, teacher values and Jacobians, and deliberately corrupted gate inputs. The analytic operator mismatch is $8.88\times10^{-16}$; FD refinement orders run from 1.994 to 2.000; the finest sine-decay check has RMSE 0.00025 °C; the affine reconstruction has maximum error $3.98\times10^{-8}$ °C; correct oracle faces pass and doubled fluxes are rejected. An end-to-end manufactured case, $T = 20 + 3\xi + 2(1 - e^{-s})\sin\pi\xi$ with $\theta = 0.20 + 0.06\sin(2\pi s/4)\cos\pi\xi$ and an analytic source, requires the three physics controls to reach target RMSE below 0.10 °C (Table 1). Two control executions returned identical field RMSEs with unchanged teachers.

6 Results

6.1 Block-level reconstruction

The operator variants have similar mean errors, and the gated variant has higher mean error. Table 2 gives equal-block means and medians of the within-block seed-mean RMSE, Table 3 the paired contrasts and Appendix 10 every block. The neural-flux arm has higher mean error than exact FD (0.591 against 0.590 °C) and the gated arm higher than raw neural flux (0.601). The matched-information numerical reconstruction, with no network, equals FD at 0.590 in a twentieth of a second. Lowest of all among complete procedures is the constant-diffusion AD ablation at 0.567, which motivates a separate assessment of the prescribed moisture multiplier.

At 0.0013 °C the primary difference has an interval spanning zero; the neural arm has lower error in three of ten blocks (two strictly at the printed precision, one tie). Only the gated-minus-FD interval excludes zero, favouring exact FD. The bootstrap is limited to three stations and one winter block.

6.2 Visited inputs, gradients and gate acceptance

Audits after Adam updates 0, 100, …, 700 save field weights, collocation coordinates and face inputs for every neural-flux and gated fit, and post-fit diagnostics use a new $97 \times 31$ grid. The mean fraction of faces inside teacher support rounds to 0.999 and oracle-gate acceptance to 0.823 in both families. The raw learned-minus-exact spatial-operator RMS discrepancy averages 0.025 °C/day. Gated maxima are 0.0969 on neural-flux states and 0.0975 on gated states, inside the 0.10 °C/day budget of Equation 4 on every sampled state. That comparison holds the field fixed under two operators; it does not show the gated field is closer to the measured target, and a thresholded objective can change the optimiser’s trajectory even when every selected value meets its budget. A spacing ladder $h = 0.10, 0.05, 0.025, 0.0125$ on each saved field, without refitting, measures the amplification channel: face errors enter the divergence over $h^2$, but changing $h$ also changes the face inputs, so the ladder need not scale as $h^{-2}$.

6.3 Parameter-gradient audit

After the outcomes were known, an additional diagnostic was applied to all 60 neural-flux and gated fits (480 training snapshots and 60 final states) without refitting or changing any score or threshold. At each identical state the exact-FD, neural-flux and gated objectives were differentiated with respect to the field parameters $\psi$, with data and initial terms identical and the gate at a fixed branch, and the relative discrepancy $\rho_M = \|\nabla_\psi \mathcal{L}_M - g_{\text{FD}}\|_2 / \max(\|g_{\text{FD}}\|_2, 10^{-12})$ recorded with absolute differences and cosines (Table 4). On neural-flux states the gate lowers the median full-objective discrepancy from 1.207 % to 0.508 % and raises the median cosine from 0.999936 to 0.999989: closer sampled gradient agreement with FD at fixed states, and no explanation of why the separately optimised gated arm has higher error. Separating adaptive updates, trajectory effects, branch switching and model mismatch needs a different experiment.

6.4 Residuals, repeatability and cost

Mean new-grid continuous residual RMS values are 4.330 (AD), 4.226 (FD), 4.146 (gated) and 4.133 (neural flux) °C/day. They describe remaining equation mismatch under the assumed model, share the trained field, equation and AD library, and certify nothing physical. For neural-flux fields the discrete-operator discrepancy of about 0.025 and the continuous residual of about 4.133 are different quantities and do not partition fitting from physical error. A fresh Python process repeated all five seed-11 fits for the first eligible block with a maximum RMSE difference of 0 against a frozen tolerance of $10^{-10}$. The gate performs exact-map and Jacobian checks online, so no acceleration is inferred; any cost claim must include offline fitting, amortisation, guard overhead, inference and accuracy on controlled hardware.

7 Discussion

The matched-field comparison isolates the spatial operator and shows similar mean reconstruction errors. Neural flux, FD and AD sit within 0.004 °C of each other, and the gate adds 0.01 °C. A sparse linear reconstruction with the same information matches the exact-flux PINN at the reported precision. Improvement over the data-only field does not establish an advantage over the conventional reconstruction.

Jointly non-linear in moisture and temperature difference, the learned map approximates an oracle that is linear in the difference given moisture, so the generic odd approximation can introduce deviations from that conditional structure, and its held-out box error does not determine behaviour on fitted-field states outside the sampled moisture range. The input, Jacobian and fixed-state gradient audits probe possible perturbations without identifying a cause. Oracle access is at once the guard’s experimental strength and its practical limit. Exact values permit a transparent rejection test and an evaluated budget, and they also remove the reason to replace the operator. An oracle-free successor needs a separately justified enclosure and must be scored on new windows, since these targets have now informed interpretation.

Physical-model uncertainty is separate. The shape parameters were not measured at each site, the constant heat capacity is an approximation, and the Asheville winter block may involve omitted processes; the constant-diffusion ablation checks one choice and the fact that it scores best is a reason to revisit the multiplier before revisiting the operator. The study establishes no coupled hydrothermal dynamics, Caribbean transfer, forecast skill or projection. What it does supply is evidence about local coupling, residual design and reproducible evaluation for any face-flux architecture, and a reminder that a satisfied operator budget is not a solution certificate.

8 Conclusion

On these ten blocks a learned moisture-to-flux map does not improve the exact-FD PINN, and an oracle-assisted gate does not improve the learned map, while the gate’s evaluated budget is met everywhere it was sampled. The frozen protocol, before-scoring controls, retained snapshots, matched-information comparator and fresh-process repeat make the scope of that conclusion inspectable. The pointwise budget bounds the evaluated operator error, not the reconstruction error.

A Coverage register

Requested start dates were 20 January, 15 April, 15 August and 15 October 2024 at four stations. Eligible blocks: Newton 15 April, 15 August, 16 October (12:00); Asheville 20 January, 15 April, 15 August, 16 October (16:00); Blackville 15 April, 15 August, 15 October. Unavailable: Newton January; all four Crossville dates; Blackville January. Each eligible block supplies 96 hourly records and 95 scored targets.

B All block comparisons

Data and code availability

USCRN hourly files are public. The archive (protocol and configuration hashes, teacher and field weights, the 480 training snapshots, audit grids, score tables and the control-chain outputs) is held by the author and is available on request.