Residual Machine Learning Compensates Turbidity Interference and Corrects In-Grid Saturation in On-Line Spectrophotometric Monitoring of Ammonia Nitrogen, Total Nitrogen and Total Phosphorus

High-resolution monitoring of ammonia nitrogen (NH$_4$-N), total nitrogen (TN), and total phosphorus (TP) underpins the management of drinking-water supply, wastewater treatment, agricultural runoff, and surface eutrophication [1], [2], [3]. Ultraviolet-visible (UV-Vis) spectrophotometry has become the dominant on-line technique for these analytes because it requires no reagents, returns results within minutes, and is sufficiently robust for unattended operation [4], [5], [6]. Its measurement principle, however, relies on the Beer-Lambert relation, which is violated by suspended particles that scatter incident light. Turbidity therefore introduces a baseline shift and a non-linear distortion in the analytical channel, which propagates through the calibration curve into the reported concentration [7], [8], [9]. The conventional response has been to add a turbidity-reference optical channel and to apply a fixed-coefficient correction: a single linear subtraction for TN, a multi-channel log-linear formula for TP, and a quadratic log-ratio embedded in the analytical expression for NH$_4$-N. The three baselines fail in characteristic ways. The TN linear formula reduces error at $\ge$100 NTU but worsens it at $<$100 nephelometric turbidity units (NTU) (13 of 16 low-turbidity samples lose accuracy after compensation). The TP embedded compensation raises the low-turbidity mean absolute error (MAE) from 0.013 to 0.022 mg·L$^{-1}$ (37 of 42 samples) while leaving the high-NTU mean absolute percentage error (MAPE) above 66%. The NH$_4$-N formula performs well between 0.43 and 1.35 mg·L$^{-1}$ but saturates at 2.35 mg·L$^{-1}$, returning a near-constant 0.46 mg·L$^{-1}$ that does not respond to turbidity. No fixed-coefficient adjustment can repair this last failure mode.

Machine learning offers an alternative to brittle fixed-coefficient compensation. Data-driven water-quality models have evolved from hydrological regression [10], [11] to tree ensembles and deep networks [12], [13], [14], [15], with applications spanning water-quality-index prediction [12], [13], [14], satellite-based nutrient retrieval [16], [17], [18], wastewater-effluent forecasting [19], [20], [21], [22], and on-line sensor calibration [23], [24], [25], [26]. eXtreme Gradient Boosting (XGBoost) in particular has become a standard reference for small-to-moderate tabular datasets that require interpretability and computational efficiency [27], [28], [29]. Within this broader literature, three methodological families have specifically addressed the turbidity-induced scattering problem in UV-Vis water-quality monitoring. The first family rests on classical chemometric corrections such as multi-wavelength absorbance subtraction, multiplicative scatter correction (MSC), orthogonal signal correction (OSC), and their direct or extended variants combined with partial-least-squares regression (direct orthogonal signal correction—partial least squares (DISC-PLS), extended multiplicative scatter correction (EMSC)), which restore an estimated clean absorbance spectrum by modelling the multiplicative and additive effects of particulate scattering before regression [30], [31]; hybrid partial least squares—artificial neural network (PLS-ANN) architectures have recently improved on this paradigm for chemical oxygen demand (COD) under turbidity [32], and rapid multi-wavelength turbidity-correction schemes have been demonstrated for COD in real water [33]. The second family applies wavelet-based multi-scale signal separation, in which a Symlet- or Daubechies-basis discrete wavelet decomposition isolates high-frequency instrument noise from the lower-frequency scattering signature in the 450–500 nm region, enabling targeted denoising or band-selective regression before downstream modelling [34], [35]. The third and most recent family applies one-dimensional convolutional neural networks (1D CNN) and U-Net architectures directly to the raw absorbance vector, with end-to-end deep learning shown to outperform OSC- and EMSC-style pre-processing for nitrate and COD quantification under random turbidity disturbance [36], and dedicated deep-learning turbidity compensators reporting substantial R2 gains for TN and TP under realistic field conditions [37], [38]. Each family carries distinct trade-offs in training-data requirement, interpretability, and on-device computational cost.

Recent UV-Vis-based application studies have further extended these methods to COD in industrial wastewater [39], nitrification monitoring [40], absorbance-based pathogen detection [41], and fluorescence excitation-emission spectroscopy [42]. Three gaps in this combined literature motivate the present study. First, methodological work has concentrated on the spectrum domain: classical chemometrics, wavelet decomposition, and 1D CNNs all act on the raw absorbance vector to estimate or denoise a clean spectrum that is then passed to a downstream regression step. The complementary strategy of retaining the certified analytical output and learning a correction in the concentration domain, which becomes particularly important when a regulated instrument's analytical formula must remain untouched, has received much less attention. Second, accuracy is usually reported only under random-fold cross-validation, which is optimistic when the data cluster by target concentration; leave-one-concentration-out (LOCO) testing, although better matched to deployment, is rarely used [43], [44]. Third, prior work on individual nutrients [2], [45], [46] has not been unified into a single compensation framework spanning NH$_4$-N, TN, and TP, so the intrinsic compensability of these analytes by AI is difficult to compare.

This study addresses the three gaps. A residual-learning architecture is proposed that retains the instrument output as a baseline and adds a learned data-driven correction. By construction, the corrected prediction reduces to the linear baseline as the turbidity proxy approaches zero, preventing degradation of already-accurate low-turbidity samples. Three model families of increasing complexity are evaluated within this architecture: a three-coefficient parametric polynomial as the interpretable baseline, a Ridge regression on engineered log-features as the production model, and a tightly regularised residual XGBoost with monotonic constraints. The framework is applied to NH$_4$-N, TN, and TP, three analytes whose baseline algorithms and turbidity sensitivities differ substantially, under controlled turbidity-addition experiments. Performance is reported under both random 5-fold and LOCO cross-validation, with metrics stratified by low- and high-turbidity groups. Trained models are further evaluated on preliminary TP field samples and off-grid intermediate NH$_4$-N concentrations to probe generalisation limits. The contributions are: (i) a deployable residual-learning compensation framework that applies to multiple nutrient analytes; (ii) a structured comparison of parametric, regularized-linear, and tree-ensemble residual models for turbidity compensation; (iii) evidence that gradient-boosted trees can correct an observed in-grid NH$_4$-N saturation case, while LOCO and off-grid tests show limited generalization to unseen NH$_4$-N concentrations; and (iv) LOCO-based accountability that identifies where AI compensation succeeds and where it cannot extrapolate, informing future data-collection priorities.

Measurements were obtained from two parallel sources. On-line spectrophotometric readings provided the raw photometric signals together with the instrument's own back-calculated concentration. Manual laboratory determinations, performed by certified analysts, provided the ground-truth reference. NH$_4$-N was measured by Nessler’s reagent spectrophotometry according to HJ 535-2009 [47], TN by alkaline potassium persulfate digestion followed by UV spectrophotometry according to HJ 636-2012 [48], and TP by ammonium molybdate spectrophotometry according to GB/T 11893-1989 [49]. The manual values constitute the regression target. The raw photometric signals and the instrument's back-calculated concentration serve as model inputs and as the baseline against which residual learning is evaluated. Turbidity additions were prepared with a formazine reference suspension and verified against a calibrated nephelometer before each experimental day. NH$_4$-N comprised four target concentrations (0.43, 0.82, 1.35, and 2.35 mg·L$^{-1}$) crossed with six turbidity levels (0, 100, 200, 300, 400, and 500 NTU) at three replicates per condition ($n$ = 72). TN comprised four target concentrations (0.5, 1.0, 2.0, and 4.0 mg·L$^{-1}$) crossed with seven turbidity levels (0, 50, 100, 200, 300, 400, and 500 NTU) on two photometric instruments ($n$ = 56). TP comprised four target concentrations (0, 0.05, 0.1, and 0.2 mg·L$^{-1}$) crossed with eight turbidity levels (0, 98, 181, 308, 407, 561, 691, and 1034 NTU) at multiple replicates per condition ($n$ = 219). Two external evaluation sets were retained for generalization testing: twelve TP field samples, comprising six raw-water samples and six matched aliquots spiked by adding 0.10 mg·P·L$^{-1}$, probing real-water matrix effects, and twenty-nine NH$_4$-N measurements at 1.6 and 2.0 mg·L$^{-1}$, retained as an off-grid intermediate-concentration set because their turbidity coverage was incomplete.

The three analytes share a common photometric front-end (analytical channel + turbidity-reference channel) but differ markedly in their existing turbidity-compensation algorithms. Conventional on-line monitoring practice converts the raw photometric signals into a reported concentration through a linear fit of the analytical absorbance (with at most a single fixed-coefficient subtraction of the turbidity-reference channel), and this approach is known to falter as soon as samples deviate from the narrow conditions under which the linear coefficients were calibrated, producing systematic biases at low turbidity, saturation at high concentration, and unstable behavior across different optical regimes. In the instruments studied here, the TN baseline applies an explicit single-coefficient linear subtraction of the turbidity-reference absorbance followed by a per-instrument linear calibration, the TP “new algorithm” produces a turbidity-aware concentration through a multivariate log-linear combination of all four signal/reference channels, and the NH$_4$-N back-calculation employs a quadratic-in-log-ratio expression in which the turbidity correction is embedded inside the analytical formula itself. Each represents a distinct engineering response to the same underlying problem, and each, as the Results below show, exhibits a different characteristic failure mode that the residual-learning correction is designed to address.

Rather than training a machine-learning model to predict concentration directly from raw signals, an approach that discards the considerable engineering effort embedded in the instrument's existing analytical formula and is brittle in deployment, we take the instrument's current concentration output as a baseline and learn an additive correction on top of it. The corrected prediction is therefore the sum of two terms: the unmodified baseline output, and a small learned residual whose magnitude depends on a vector of engineered photometric features. Three properties make this architecture safe for on-site use. It is additive, so any failure of the learned residual (e.g., a numerical outlier or a missing input) at worst reproduces the existing instrument output rather than introducing arbitrary error. It is trained against the deviation between the manual ground-truth concentration and the baseline output, a target with near-zero mean that is closer to a Gaussian regression problem than the raw concentration. And it satisfies a zero-turbidity safety property, enforced through feature design and, for the tree-ensemble variant, monotone constraints, whereby the correction collapses to zero whenever the turbidity proxy signal is small, preventing degradation of already-accurate clean-water samples [13], [25]. The features used by the residual model are constructed to expose photometric relationships in an additive form: log-transforms of the analytical, baseline and turbidity channels and their references; the log ratio that serves as the analyte-specific turbidity proxy; the instrument's own back-calculated baseline output; and analyte-specific intensity contrasts such as the analytical-to-baseline signal ratio (for NH$_4$-N and TP) and the turbidity-to-analytical contrast (for NH$_4$-N), each chosen to be scale-invariant to source-lamp intensity drift.

The first model family, Approach A (parametric polynomial), implements the residual as a three-coefficient function (a constant, a linear and a quadratic term) of a single analyte-specific turbidity proxy: the turbidity-reference absorbance for TN, the logarithm of the reference-channel ratio for TP, and the logarithm of the turbidity-to-reference signal ratio for NH$_4$-N. The three coefficients are estimated by Levenberg-Marquardt non-linear least squares with the no-correction point as the initial guess, so the model degenerates to the existing baseline whenever the data do not support a non-trivial correction. With only three parameters, this approach is the most interpretable of the three and serves as the physics-aware comparison baseline against which more flexible models are judged.

The second model family, Approach B (Ridge regression, the production model), uses the eight engineered log-features described above and passes them, after standardization, to a ridge-regression model whose L2 penalty $\alpha$ is selected from the candidate set {0.01, 0.1, 1, 10, 100} by 5-fold inner cross-validation [14], [50]. For every analyte, the chosen penalty lies in the moderate-to-strong range, indicating that the residual signal is small relative to the engineered-feature variance and benefits from substantial shrinkage on the few-hundred-sample training sets. The third model family, Approach C (residual XGBoost with monotone constraints), augments the eight engineered features with the underlying raw signals to give twelve to fourteen inputs depending on analyte, and trains a tightly regularized gradient-boosted tree with maximum tree depth fixed to two, between thirty and fifty estimators, an L2 leaf-weight penalty of ten, a minimum child weight of five, a learning rate of 0.1, the histogram tree-construction method and a single thread for reproducibility [27], [28]. To prevent the tree model from exploiting correlated features in a way that would violate the zero-turbidity safety property, monotone-non-decreasing constraints are imposed on every input feature that co-varies positively with the analyte-specific turbidity proxy. Across the three analytes, the model interfaces are deliberately identical (fit, predict, and a predict with baseline method that exposes the baseline, the learned residual, and the final non-negativity-clipped output), so that a single production inference class can load any of the three trained artefacts.

Two complementary cross-validation schemes were applied independently to every (analyte, model) pair to separate optimistic from deployment-realistic accuracy estimates [43]. Random 5-fold cross-validation, with a fixed random seed for reproducibility, distributes samples uniformly across folds and reproduces the conventional performance number reported in most water-quality machine-learning studies; five folds rather than ten were used so that the smaller TN ($n$ = 56) and NH$_4$-N ($n$ = 72) datasets retained enough training samples per fold without inflating per-fold variance. LOCO cross-validation holds out all measurements at a single target concentration as the test fold, forcing the model to extrapolate to an unseen concentration regime, a scheme that is informative because random folds inevitably leak concentration information whenever replicates of the same target are split across folds, flattering accuracy. Four error metrics were tracked per fold: MAE and root-mean-square error as the primary deployment-relevant quantities, MAPE computed only on rows with positive ground-truth concentration to avoid the singularity at zero, and the maximum absolute error to detect rare large failures. All four metrics were further stratified into low-turbidity ($<$100 NTU) and high-turbidity ($\ge$100 NTU) sub-groups in addition to the overall figure, exposing differential behavior across the two operationally meaningful regimes [14], [19]. A safety-floor check flagging any (analyte, model) pair whose MAE on an individual fold exceeded 1.25× the corresponding linear-baseline MAE was applied as a pre-deployment regression guard, and pair-wise per-sample absolute-error differences between Approach B and Approach C were tested with a paired Wilcoxon signed-rank test at $\alpha$ = 0.05, Bonferroni-corrected over the three analytes.

Hyper-parameter selection followed a nested design that prevents test-set leakage. The three polynomial coefficients of Approach A were estimated by Levenberg-Marquardt non-linear least squares with the zero-correction solution as the initial guess, so the model degenerates to the existing baseline whenever the data do not support a non-trivial correction. The L2 penalty of Approach B was selected, on each outer training fold, from the grid {0.01, 0.1, 1, 10, 100} by an inner 5-fold cross-validation against MAE; the selected values clustered in the moderate-to-strong range (100 for TN, 0.1 for TP, 1.0 for NH$_4$-N), reflecting differential noise-to-signal ratios across the three analytes. Inputs to Approach B were z-score standardized on each training fold, with the standardizer parameters never updated from the held-out fold. The XGBoost hyper-parameters of Approach C were fixed across all three analytes at the conservatively regularized values described in Section 2.2, with the analyte-specific monotone-constraint vectors that enforce the zero-turbidity safety property documented in the manuscript repository. SHapley Additive exPlanations (SHAP) attribution used the tree-path-dependent TreeExplainer with the full training fold as the background distribution. After cross-validation, models trained on the complete grid were applied without any further refitting to two independent withheld evaluation sets: twelve TP field samples (six unspiked, six 0.1 mg·L$^{-1}$ spike) probing matrix-effect robustness on real water rather than formazine-spiked deionized water; and twenty-nine NH$_4$-N measurements at 1.6 and 2.0 mg·L$^{-1}$ probing off-grid concentration generalization within the calibrated NH$_4$-N span.

The full training and inference pipeline is implemented in Python 3.12 using pandas and NumPy for data processing, SciPy's optimize module for Approach A, scikit-learn 1.4 (Ridge, KFold, StandardScaler) for Approach B, XGBoost 3.2 for Approach C, SHAP 0.45 for feature attribution, Joblib for model serialization, and Matplotlib 3.8 for figure generation [51]. The pipeline is organized as eight modules per analyte under a unified package layout, with seventy-six-unit tests covering data ingestion, cross-validation, model contracts, and the production inference API. End-to-end execution from raw Excel files to all trained artefacts and figures takes approximately three minutes on commodity workstation hardware. Random seeds, train/test splits, hyper-parameter grids and the deterministic-XGBoost settings are version-controlled in the manuscript repository, enabling bit-exact reproduction of every result reported in Section 3.

The TN linear formula A3 = A1 + 0.35·A2 of Qi et al. [9] reduces the overall MAE from the uncompensated 0.149 mg·L$^{-1}$ to 0.041 mg·L$^{-1}$. The reduction is non-uniform across the turbidity range. At $<$100 NTU, the MAE rises from 0.032 to 0.036 mg·L$^{-1}$, and 13 of the 16 low-turbidity samples lose accuracy after compensation; at $\ge$100 NTU, the MAE falls from 0.225 to 0.043 mg·L$^{-1}$ (Figure 1a). Closing this non-uniformity is the design objective of residual learning. Cross-validation results for all four methods are summarized in Table 1 and Figure 2a. The Ridge production model (Approach B) attains a random-fold overall MAE of 0.0237 mg·L$^{-1}$, a 42.5% reduction from the linear baseline; the residual XGBoost (Approach C) attains 0.0226 mg·L$^{-1}$ (45.2%). Both reductions are significant against the linear baseline at the Bonferroni-corrected $\alpha$ = 0.05 level on the paired Wilcoxon signed-rank test of Section 2.3. The gains do not come at the expense of low-turbidity accuracy: the low-NTU MAE falls from 0.036 mg·L$^{-1}$ (linear) to 0.025 mg·L$^{-1}$ (Approach B) and 0.024 mg·L$^{-1}$ (Approach C), and the high-NTU MAE falls from 0.043 to 0.023 and 0.022 mg·L$^{-1}$, respectively. The residual models therefore satisfy the criterion of preserving low-turbidity accuracy, which no fixed-coefficient adjustment in the existing TN compensation literature [8], [9] has been able to meet. Under the stricter LOCO scheme, both models retain robust gains (Ridge 0.029 mg·L$^{-1}$, 30.4% reduction; XGBoost 0.026 mg·L$^{-1}$, 36.9%). The 7 percentage-point drop relative to random-fold is comparable to the optimism gap reported by Valladares Castellanos et al. [43] for spatial-LOCO in ecosystem-scale water-quality calibration.

The parametric Approach A under-performs the linear baseline (random-fold MAE 0.054 mg·L$^{-1}$ vs 0.041 mg·L$^{-1}$; LOCO 0.063 mg·L$^{-1}$). The cause is structural rather than algorithmic. Approach A relies on a per-instrument calibration recovered from the existing baseline values. Even with the recovery tuned to its analytical minimum (root-mean-square calibration residual of 0.039 mg·L$^{-1}$, comparable to the per-fold standard deviation of the manual reference itself), the residual acts as an irreducible noise floor that absorbs the gain the polynomial $k$(A2) would otherwise generate. The limitation differs from those reported for parametric machine learning (ML) water-quality correctors elsewhere [13], [50], where parametric form-mismatch rather than calibration-recovery noise has been the dominant constraint. The learned $k$(A2) curve (Figure 3a) nevertheless retains diagnostic value. The empirically optimal coefficient is $\approx$ 0.16 at the lowest A2 values (clean water) and rises to $\approx$ 0.30 at A2 $\approx$ 0.2, consistently below the conventional fixed 0.35 of Qi et al. [9]. The concentration-aware variation gives a physically interpretable explanation for the low-NTU bias observed above: the fixed 0.35 over-corrects clean-water samples, whereas the data-driven $k$(A2) would dampen the correction precisely where it is most damaging. The curve is therefore best read as evidence that the direction of correction in the existing analytical formula is correct, while the magnitude calibration would benefit from a turbidity-dependent rather than fixed coefficient. The Ridge and XGBoost residual models implement this refinement implicitly through their richer feature sets.

Mechanistic interpretation of Approach C is provided by the SHAP feature attribution in Figure 4a. The analytical-channel intensity signal1_digest and the analytical-to-baseline contrast signal3_base/digest carry the largest mean absolute contributions (0.0085 and 0.0061 mg·L$^{-1}$, respectively), more than four times the contribution of any single absorbance feature. The model therefore extracts its predictive power from the absolute photometric intensities rather than from the pre-computed absorbance values that the linear baseline uses. This feature-engineering pattern is consistent with the SHAP results reported by Makumbura et al. [44] for explainable water-quality regressors and contrasts with the absorbance-centric attributions of Lyu et al. [6] for UV-Vis nitrogen retrieval, suggesting that the relative importance of intensity and absorbance depends on the structure of the embedded baseline algorithm. Residual distributions (Figure 5a) reinforce the per-fold MAE summary. The linear-baseline high-NTU residuals are biased negative (median = $-$0.034 mg·L$^{-1}$), indicating systematic under-prediction at high turbidity, whereas the Ridge and XGBoost residual histograms are tightly centered on zero (median $\lvert \cdot \rvert$ $<$ 0.005 mg·L$^{-1}$) with comparable inter-quartile spread (0.030 mg·L$^{-1}$). The residual-learning framework therefore removes the systematic component of the baseline error rather than merely scaling it down.

The TP-embedded multivariate log-linear compensation reduces the overall MAE from the uncompensated 0.137 to 0.051 mg·L$^{-1}$. Turbidity-stratified results (Figure 1b) reveal the same low-NTU regression observed for TN. At $<$100 NTU, the embedded compensation raises the MAE from 0.013 to 0.022 mg·L$^{-1}$ (37 of the 42 low-turbidity samples worsen), whereas at $\ge$100 NTU it cuts the MAE from 0.113 to 0.058 mg·L$^{-1}$. With $n$ = 219 grid samples, the largest of the three analyte datasets and on the same order as the sample sizes used in recent TP-targeted ML studies [2], [29], the TP problem supports the most flexible AI correction of the three. Table 2 and Figure 2b show that all three residual models outperform the embedded baseline (random-fold overall MAE reductions of 78.7% for the parametric Approach A, 75.6% for the Ridge production model B, and 87.2% for the XGBoost variant C; paired Wilcoxon p $<$ 0.001 against the linear baseline after Bonferroni correction for each model). The XGBoost random-fold MAE of 0.0066 mg·L$^{-1}$ is about a factor of three below the best phosphorus-prediction MAE reported by Wang et al.[29] for their tree-based stream-phosphorus model (MAE $\approx$ 0.020 mg·L$^{-1}$). The gap is consistent with the residual-learning task being easier than the direct-prediction task addressed by Wang et al., since the residual model corrects a known baseline rather than reconstructing the concentration mapping from scratch.

An unexpected feature of the TP results is that the simplest three-coefficient parametric model is the most LOCO-robust of the four methods (LOCO MAE 0.0112 mg·L$^{-1}$, indistinguishable from its random-fold value of 0.0110 mg·L$^{-1}$). The ordering inverts the TN case, in which the parametric model was the worst, and reflects a single structural difference: TP does not require per-instrument calibration recovery (the embedded algorithm's concentration output is used directly as the baseline y_linear), so the 0.04 mg·L$^{-1}$ calibration-noise floor that limited Approach A on TN is absent. The fitted polynomial (Figure 3b) acts as a quadratic correction in the reference-channel asymmetry, an indirect turbidity indicator, whose magnitude grows rapidly as the asymmetry departs from zero. Approach C delivers the lowest random-fold MAE (0.0066 mg·L$^{-1}$) but loses 63% of that advantage under LOCO (0.0177 mg·L$^{-1}$ vs Approach A's 0.0112 mg·L$^{-1}$). The extrapolation gap is consistent with the within-vs-across concentration generalization pattern reported by Valladares-Castellanos et al. [43] for tree-ensemble water-quality models. The SHAP summary for TP (Figure 4b) identifies log10 (signal_digest) (mean $\lvert \text{SHAP} \rvert$ = 0.0042 mg·L$^{-1}$) and the existing baseline output y_linear (mean $\lvert \text{SHAP} \rvert$ = 0.0029 mg·L$^{-1}$) as the dominant features, consistent with the residual-learning interpretation that the model applies a moderate concentration-dependent correction on top of the embedded algorithm rather than re-deriving the concentration mapping.

A preliminary field-sample check on twelve independent TP samples yielded a small but informative discrepancy. The six matched aliquots spiked by adding 0.10 mg·P·L$^{-1}$ were predicted at mean values of 0.180 ± 0.012 (linear baseline), 0.179 ± 0.013 (Approach A), 0.206 ± 0.015 (Approach B), and 0.200 ± 0.014 (Approach C) mg·L$^{-1}$. All four predictions exceeded the added 0.10 mg·P·L$^{-1}$ spike amount, and this positive offset was already present in the unmodified baseline. Because background TP in the raw water was not independently certified, the offset cannot be interpreted as a model-specific error; plausible contributors include background TP, matrix-effect amplification of the spike [7], [42], or a spike-recovery issue introduced during sample preparation, rather than a model error. Repetition of the field-spike test with a freshly prepared standard and a known-recovery internal reference is recommended before any conclusion is drawn about absolute accuracy on field water. The six unspiked raw-water samples, although their true concentration is unknown, exhibit tight method-internal reproducibility (per-method standard deviation 0.002–0.005 mg·L$^{-1}$, below the random-fold MAE of any of the four methods), indicating output consistency on these field samples. Reference measurements or controlled recovery tests are needed before drawing conclusions about external accuracy or concentration tracking in real water.

The NH$_4$-N quadratic-in-log-ratio formula embedded in the instrument performs well across the calibrated working range: at 0.43, 0.82, and 1.35 mg·L$^{-1}$, the MAE stays below 0.12 mg·L$^{-1}$ at every turbidity level tested. The formula saturates at 2.35 mg·L$^{-1}$, where the back-calculated mean concentration is 0.46 ± 0.04 mg·L$^{-1}$, a 1.89 mg·L$^{-1}$ systematic under-estimation that is independent of whether the added turbidity is 0 or 500 NTU. The invariance is mechanistically diagnostic. A turbidity-induced error would vary with NTU, whereas the observed near-constant output is the signature of a formula that has exhausted its calibrated range and returned its asymptotic ceiling value. The NH$_4$-N failure at 2.35 mg·L$^{-1}$ is therefore a calibration-range saturation problem rather than a turbidity-compensation problem. The distinction does not appear to have been drawn explicitly in the on-line spectrophotometric NH$_4$-N literature [2], [3], and it has direct consequences for the choice of AI correction architecture.

Cross-validation results and the per-concentration breakdown are presented jointly in Table 3. The predicted-vs-true panel for NH$_4$-N appears as Figure 2c, and per-NTU curves stratified by target concentration appear in Figure 6. Approach C achieves an overall random-fold MAE of 0.125 mg·L$^{-1}$, a 75.6% reduction relative to the 0.511 mg·L$^{-1}$ baseline (Bonferroni-corrected paired Wilcoxon p $<$ 0.001). The aggregate figure is dominated by the 2.35 mg·L$^{-1}$ samples, where the embedded baseline formula shows saturation failure. The most consequential single statistic in this study is the reduction in 2.35 mg·L$^{-1}$ MAE from the baseline’s 1.885 to 0.055 mg·L$^{-1}$ under Approach C, a 97.1% reduction that is accompanied by improved accuracy at the lower in-grid NH$_4$-N concentrations of 0.43, 0.82 and 1.35 mg·L$^{-1}$ (MAE falling from 0.059, 0.038 and 0.062 to 0.039, 0.016 and 0.020 mg·L$^{-1}$, respectively). Approach C therefore corrected the in-grid high-concentration saturation case without degrading accuracy at the lower in-grid NH$_4$-N concentrations. To our knowledge, previous on-line water-quality ML studies have rarely evaluated whether residual machine-learning models can correct an in-grid saturation failure in an embedded analytical formula. Earlier spectrophotometric water-quality ML studies [1], [5], [6], [42] concentrated their samples inside the analyser’s calibrated range, so saturation behaviour could not be surfaced by their experimental designs.

Approach A and Approach B under-perform the linear baseline within the 0.43-1.35 mg·L$^{-1}$ working range (working-range MAE: A 0.359 mg·L$^{-1}$, B 0.341 mg·L$^{-1}$, baseline 0.053 mg·L$^{-1}$). The mechanism is a global fit captured by the large 2.35 mg·L$^{-1}$ residual, which is 30 times larger in magnitude than the residuals at lower concentrations. The regression therefore shifts the overall prediction surface to flatten the 2.35 mg·L$^{-1}$ error, at the cost of a working-range bias of comparable absolute magnitude. Approach C, by contrast, partitions the input space into concentration-aware leaves and applies a large correction only when the raw signals identify the saturated regime, leaving the working-range predictions almost unchanged. The SHAP summary for NH$_4$-N (Figure 4c) is dominated by log_ratio_digest (mean $\lvert \text{SHAP} \rvert$ = 0.158 mg·L$^{-1}$) and log_x (the embedded analytical formula's log-ratio), the analytical log-ratio that drives the embedded formula (mean $\lvert \text{SHAP} \rvert$ = 0.121 mg·L$^{-1}$). The XGBoost model has rediscovered the same photometric quantity that the original analytical chemists used and is using it to detect saturation rather than to re-derive the concentration. Residual histograms (Figure 5c) confirm the picture: only Approach C achieves near-zero residual centering at both low and high NTU (median $\lvert \text{residual} \rvert$ $<$ 0.04 mg·L$^{-1}$), whereas Approach B's residuals at 2.35 mg·L$^{-1}$ remain biased by 0.4 mg·L$^{-1}$. The practical implication for instrument designers is that complex calibration-range failures require flexible non-linear correctors, and the choice of model family should be guided by the structure of the baseline failure (smoothness for TN/TP, discontinuity for NH$_4$-N) rather than by training-set size alone.

LOCO results for NH$_4$-N expose a clear generalization limit to unseen concentration levels. All three residual models degrade substantially under LOCO testing: Approach A 0.803, Approach B 1.838, and Approach C 0.680 mg·L$^{-1}$, compared with the random-fold values of 0.541, 0.472, and 0.125 mg·L$^{-1}$, respectively. Approach B's LOCO MAE exceeds the 1.25× linear-baseline safety floor (0.511 × 1.25 = 0.639 mg·L$^{-1}$) by a factor of 2.9. This regression flag is reported transparently here in line with the deployment-realistic cross-validation practice advocated by Valladares-Castellanos et al. [43] and Makumbura et al. [44]. Approach C remains best of the three (LOCO 0.680 mg·L$^{-1}$) but is itself 444% above its random-fold MAE, an optimism gap larger than the 30–60% LOCO degradations reported for water-quality-index models on more uniform concentration distributions [13], [14]. The principal cause is the concentration density of the training grid. With four target concentrations spanning a 5.5× dynamic range, holding out 2.35 mg·L$^{-1}$ leaves the remaining three concentrations clustered between 0.43 and 1.35 mg·L$^{-1}$, a range from which no model can learn the saturation behaviour at the held-out point. Denser sampling between 1.35 and 2.35 mg·L$^{-1}$ would directly address the limit.

Off-grid intermediate-concentration evaluation on the twenty-nine measurements at 1.6 and 2.0 mg·L$^{-1}$, which lie within the 0.43–2.35 mg·L$^{-1}$ NH$_4$-N span but were excluded from the main factorial grid because their turbidity coverage was incomplete, shows a reversal of the in-grid model ranking. The parametric Approach A delivers the lowest MAE (0.189 mg·L$^{-1}$), followed by Approach B (0.260), the linear baseline (1.233), and Approach C (1.137 mg·L$^{-1}$). The tree-based model, which was strongest inside the training grid, generalizes poorly to these off-grid intermediate concentrations. The smooth polynomial $\delta$($t$) (Figure 3c) interpolates more stably across the gap between the 1.35 and 2.35 mg·L$^{-1}$ grid levels, and the rigid baseline reverts to its saturated near-constant value. The in-grid versus off-grid intermediate-concentration reversal (XGBoost $>$ Ridge $>$ Parametric inside the grid; Parametric $>$ Ridge $>$ XGBoost outside) is the trade-off between expressivity within the training distribution and graceful extrapolation outside it that Khatri et al. [25] and Wang et al. [29] identified as a central design tension in deployed water-quality ML, and that Valladares-Castellanos et al. [43] characterized quantitatively across spatial-LOCO scenarios. The present results extend that discussion in a chemically meaningful way: the choice of correction model is not only a bias-variance trade-off but also an operational-range trade-off. Instrument designers must therefore know in advance whether the deployment will encounter concentrations beyond the densest part of the training grid.

Cross-analyte synthesis. Aggregating the three analytes, four patterns emerge that have direct implications for residual-learning deployment in spectrophotometric water-quality monitoring. (i) The residual-learning architecture yields substantial random-fold MAE reductions for every analyte: 45.2% (TN, $n$ = 56), 87.2% (TP, $n$ = 219), and 75.6% (NH$_4$-N, $n$ = 72) for the XGBoost variant, and 42.5%, 75.6%, and 7.6% for the production Ridge model. The gains are of comparable magnitude to the best ML water-quality-index results reported by Banda et al. [14] and Makumbura et al. [44] on similarly sized datasets, but achieved here on the more demanding task of correcting a structured baseline error rather than predicting concentration from scratch. (ii) The optimal model family is determined by the structure of the baseline error rather than by dataset size: TN favors XGBoost marginally over Ridge (smoothly turbidity-dependent error, both models adequate); TP favors the parametric model under LOCO (smoothly turbidity-dependent error in a calibration-noise-free baseline, simplest model generalizes best); NH$_4$-N requires the flexibility of XGBoost (concentration-dependent saturation, only piecewise-flexible models can isolate the failure regime). (iii) LOCO reporting changes the conclusion for NH$_4$-N more than for TN or TP. The Ridge model that ranks second on random-fold becomes the worst on LOCO, confirming the concern raised by Valladares-Castellanos et al. [43] that random-fold-only reporting can substantially overstate field accuracy. (iv) Data size and baseline-error structure interact: with the largest dataset (TP), the simplest model wins; with the smallest dataset (TN), the moderate-complexity model wins; and with the most pathological baseline (NH$_4$-N), only the most flexible model achieves the largest gains, at the cost of the weakest LOCO performance.

Quantitative feature attribution from Ridge coefficients and XGBoost gain. The Ridge solution exposes the linear weight that each engineered log-feature receives, and the XGBoost gain and weight statistics expose how often, and with what split-quality benefit, each feature is used in the tree ensemble. The two intrinsic rankings, complementary to the SHAP attributions discussed below, give a quantitative answer to the question of whether the baseline analytical output or the bare turbidity proxy drives each model. For TN, the Ridge coefficient budget is concentrated on the bilinear turbidity term A1·A2 (31.9% of total $\lvert \text{coefficient} \rvert$) and its higher-order companions A22 (18.1%) and A2 (12.8%), with the embedded baseline output y_linear carrying 17.8% and the residual quartile distributed across A1 and log(A2). The XGBoost gain ranking is led by the raw analytical-channel intensity signal1_digest (51.2% of total gain), followed by the analytical-to-reference channel pair signal3_base and A1 (12.4–12.5% each). For TN, the residual model is therefore driven by the turbidity proxy and its photometric correlates rather than by the baseline analytical output. For TP, the picture inverts. The Ridge solution distributes its coefficient budget across log_signal_digest (29.7%), signal_ratio (28.9%), y_linear (19.3%), and absorbance_new (18.4%), while the bare turbidity proxy log_ref_ratio receives only 2.9%. The XGBoost gain ranking places y_linear first (31.4%), followed by the dual-reference logarithms log_ref1 (21.5%) and log_ref2 (17.8%), with log_ref_ratio again at only 2.7%. For TP, the residual model relies primarily on the embedded baseline output and the individual reference channels rather than on their ratio. For NH$_4$-N, the two model families adopt distinct but internally consistent strategies. The Ridge solution gives 51.7% of its coefficient budget to the analytical-channel ratio analytical_ratio, with the next two contributors being intensity_contrast (12.7%) and the turbidity proxy log_ratio_turb (10.5%). The XGBoost gain ranking, by contrast, concentrates 99.3% of its predictive weight on the embedded formula's own log-ratio quantities log_ratio_digest (67.0%) and log_x (32.3%), with all other features below 1%. Quantitative attribution therefore yields three qualitatively different answers across the three analytes: TN is driven by the turbidity proxy and its interactions, TP is driven by the baseline analytical output and the dual reference channels, and NH$_4$-N is driven by the embedded formula's central log-ratio.

Cross-analyte feature attribution. The three SHAP summaries (Figure 4) reveal a consistent picture of what the XGBoost residual model has learned for each analyte, and of why the same architecture selects different photometric quantities depending on the structure of the baseline error. For TN (panel a), the two largest mean absolute SHAP contributions are signal1_digest (0.0085 mg·L$^{-1}$) and the analytical-to-baseline intensity contrast signal3_base/digest (0.0061 mg·L$^{-1}$), which exceed any single pre-computed absorbance feature. The model has therefore extracted its predictive power from absolute photometric intensities rather than from the absorbance values the linear baseline already uses, a pattern consistent with the SHAP-based feature rankings of Makumbura et al. [44] for explainable water-quality regressors. For TP (panel b), log10(signal_digest) (0.0042 mg·L$^{-1}$) and the existing baseline output y_linear (0.0029 mg·L$^{-1}$) are dominant, consistent with the residual model adding a moderate concentration-dependent correction on top of the embedded algorithm rather than re-deriving the concentration from raw signals. For NH$_4$-N (panel c), the analytical log-ratio log_ratio_digest (0.158 mg·L$^{-1}$) and the embedded formula's own log_x term (0.121 mg·L$^{-1}$) dominate by more than an order of magnitude over any other feature. The model has rediscovered the photometric quantity that the analytical chemists embedded in the original formula and is using it as a saturation detector rather than as a concentration regressor, a mechanistic insight that could not be inferred from MAE numbers alone. The SHAP results therefore reinforce the model-selection narrative of Section 3 (TN and TP show a smooth dependence on intensity contrasts; NH$_4$-N shows a discontinuous dependence on the calibrated log-ratio) and provide a transparency layer consistent with the explainability requirements increasingly demanded for AI in environmental measurement [42], [44], [46].

The most informative pattern emerging from the cross-analyte results is not any single accuracy figure but a structural asymmetry: the residual-learning framework that closes the turbidity gap for TP survives LOCO testing essentially intact, whereas the analogous framework for NH$_4$-N collapses once the training grid loses information about the over-range regime. The asymmetry is not a property of the residual algorithm or the cross-validation protocol. It is a predictable consequence of three physicochemical factors operating jointly: the position of each analytical reaction on the Beer-Lambert linearity curve, the spectral character of the turbidity interference at the wavelengths the two methods use, and the geometry of the laboratory concentration grid relative to the analyzer’s calibrated dynamic range.

The first factor concerns analytical chemistry. TP is determined here by the molybdenum-blue reaction, in which orthophosphate combines stoichiometrically with molybdate and a reductant to form a heteropoly complex with a sharp absorbance maximum near 700 nm and a molar absorptivity of roughly 2.5 × 104 L·mol$^{-1}$·cm$^{-1}$. At the concentrations of interest for drinking water and surface water, this places the analytical channel comfortably below an absorbance of 0.5 in the 1-cm cell used by the analyzer, well inside the linear region of the Beer-Lambert relation. Within that region, the analytical absorbance is, to a good approximation, an affine function of phosphorus concentration, and the turbidity interference manifests as a separable, multiplicative scatter term that the existing log-linear compensation already captures approximately. NH$_4$-N, by contrast, is measured here through direct UV-Vis absorption without a colorimetric amplification step, which forces the analyzer to work at higher absorbance values in order to achieve adequate signal-to-noise. At the upper end of the experimental range, the analytical absorbance approaches the photodetector's linear ceiling, where stray-light leakage from the reference channel, photon-shot-noise floor, and analog-to-digital quantization jointly compress the apparent signal. The transition into this regime is not gradual; it is a piecewise non-linearity that activates only once the analytical absorbance crosses a threshold determined by the photometric front-end. The baseline distortion at the saturated concentration is therefore a detector-side artefact rather than a turbidity artefact, and it is mathematically a different object from the smooth multiplicative perturbation that TP experiences.

The second factor concerns the spectral character of particulate scattering at the two analytical wavelengths. The 690–710 nm region used for the TP measurement sits well into the visible part of the spectrum, where suspended particulates produce a Mie/Rayleigh scattering response that is approximately proportional across the analytical and reference channels for the range of particle sizes generated by formazine and by the natural particulates that formazine simulates. A smooth turbidity proxy therefore captures the TP scatter accurately, and a low-order polynomial residual closes the remaining gap. The UV and near-UV wavelengths used for NH$_4$-N lie at the steeply rising edge of the scattering cross-section spectrum, where Mie/Rayleigh scaling is much more sensitive to particle size distribution, and where chromophoric dissolved organic matter introduces a wavelength-selective absorbance contribution that interacts with the ammonia signal itself. The composite error has both multiplicative and additive components, and these components scale differently with concentration. A scatter term that is approximately separable for TP becomes only partially separable for NH$_4$-N, and the residual model must therefore learn an interaction surface that is local to each concentration band.

The third factor, and the one that ultimately determines LOCO behavior, is the geometry of the laboratory concentration grid relative to the analyzer’s calibrated range. The four TP target concentrations are arranged so that adjacent levels are separated by a small fraction of the calibrated upper limit, and so that holding out any single concentration always leaves training points on both sides of the held-out value. The TP LOCO test is therefore an interpolation problem in the strict sense: the held-out value lies inside the convex hull of the training distribution. The four NH$_4$-N target concentrations do not enjoy this geometry. The gap separating the highest concentration from the next-lower one is wider than the spacing among the three lower levels combined, and the highest concentration is also the only one in which the detector enters its non-linear regime. Holding out the highest concentration therefore strips the training set of both its upper geometric anchor and the only sample exhibiting the qualitative behavior the residual model would need to extrapolate. LOCO for NH$_4$-N is consequently not an interpolation problem at all; it is a true out-of-distribution extrapolation that requires the model to predict a regime it has never seen.

These three factors combine into a single mechanistic explanation. TP residuals are smooth, separable, and concentration-stationary because the analytical chemistry stays inside the Beer-Lambert linear region, the scattering is approximately multiplicative at the analytical wavelength, and the concentration grid keeps every LOCO fold inside the convex hull of the training data. NH$_4$-N residuals are piecewise, non-separable, and concentration-non-stationary because the analytical method drives the detector into its non-linear ceiling, the scattering interaction adds a wavelength-dependent component, and the concentration grid forces LOCO into genuine extrapolation. The TP/NH$_4$-N LOCO contrast is therefore not a statistical artefact of the cross-validation protocol or a deficiency of any particular residual model. It is the predictable signature of the joint analytical chemistry, photometric physics, and experimental geometry under which the two analytes were measured.

The physicochemical reading developed above clarifies why a residual architecture is the natural choice across both regimes. Where the analytical chemistry sits inside its linear region, and the embedded compensation has approximately the right functional form, the certified analytical output already encodes most of the relevant information, and a residual layer on a few engineered photometric features supplies a smooth, physically interpretable refinement. Where the analytical chemistry drives the detector into its non-linear ceiling and the embedded formula fails outside its calibration band, a residual layer can still preserve the formula's correctness in the working range while learning a separate, piecewise correction that activates only when the raw signals identify the failure regime. The same architectural decision, namely to retain the certified analytical output as a baseline and to constrain the learned residual to vanish at zero turbidity, therefore behaves as a smooth refinement in one regime and as a discontinuity-aware backstop in the other, without any structural change in the inference pipeline.

This contrasts with the two alternative paradigms most prominent in the recent spectrophotometric water-quality literature. Spectrum-domain deep learning architectures such as 1D CNNs and U-Nets reconstruct a denoised absorbance vector from the raw signal before passing it to a downstream regression step [36], [37], [38], and chemometric pre-processing schemes (MSC, OSC, EMSC, DOSC-PLS) achieve the same end through explicit multiplicative- and additive-scatter modelling [30], [31], [32]. Wavelet-based decomposition is sometimes layered on top of either family to separate high-frequency instrument noise from lower-frequency scattering contributions [34], [35]. All of these methods operate on the assumption that the analytical formula downstream of the cleaned spectrum will behave correctly once the spectral interference has been removed. The TP/NH$_4$-N contrast shows that the assumption holds only when the analytical chemistry stays inside its linear region. When detector non-linearity drives the baseline into saturation, no amount of spectrum-domain cleaning can recover information that the saturated analytical formula has already discarded; the residual must be learned in the concentration domain on top of, rather than in place of, the certified analytical output. The two paradigms are therefore best regarded as complementary: spectrum-domain methods address scattering and noise in the optical front-end, while concentration-domain residual learning addresses non-linearity and saturation in the analytical back-end. A hybrid pipeline that combines wavelet- or CNN-based spectral denoising with a downstream residual learner is a natural direction for future work, particularly for analytes in which both error sources are operationally relevant.

Practical model choice within the residual paradigm follows the same physicochemical reasoning. Smoothly turbidity-dependent residuals admit a parametric or Ridge-regularized representation whose extrapolation behavior is bounded by the analytic form of the function and is therefore well behaved at concentrations outside the training grid. Piecewise residuals driven by detector non-linearity admit no such bound and require flexible tree-based models that can carve the input space into concentration-aware regions, but those same models lose graceful extrapolation outside the densest part of the training distribution. The reversal of model rankings between in-grid and off-grid evaluation for NH$_4$-N is not paradoxical; it is the bias-variance trade-off expressed in the analytical chemistry of the measurement, with a direct implication for the choice of model family at deployment.

The Ridge coefficients and the XGBoost gain rankings reported quantitatively in Section 3 admit a clean physico-chemical interpretation that is itself analyte-specific. For TN, the certified analytical formula reduces to a single-coefficient linear subtraction of the turbidity-reference absorbance from the analytical-channel reading, and the baseline output therefore encodes only that one fixed-coefficient operation. The residual model has access to no information that the baseline has not already passed through, so it redistributes its predictive weight onto the bare turbidity-reference absorbance, its quadratic, and its interaction with the analytical channel, exactly the photometric quantities that a turbidity-aware linear formula leaves unexploited. The baseline output y_linear is therefore not the central feature for TN, because the chemistry encoded in y_linear is too thin to dominate. For TP, the embedded multivariate log-linear compensation already absorbs the dominant log-multiplicative coupling between the analytical and reference channels into a single concentration estimate. The baseline output y_linear is consequently a dense, information-rich summary of the analytical chemistry, and the residual model rationally treats it as the principal feature, with each reference channel contributing independent, non-collinear corrections rather than entering as a ratio. The turbidity proxy log_ref_ratio is correspondingly marginal here, because the proxy information has already been absorbed into the baseline. For NH$_4$-N, the embedded quadratic-in-log-ratio expression places the analytical log-ratio at the centre of the analytical computation, and both the Ridge and XGBoost residual models converge on that same quantity. The Ridge model expresses this convergence through a large linear weight on the analytical-channel ratio; the XGBoost model expresses it through almost exclusive reliance on the embedded log-ratio and the formula's own log-domain variable. The unifying principle across the three analytes is that residual learning rewards the analytical chemistry the baseline has already performed: the richer the chemistry encoded in the baseline output, the more the residual treats that output as a single dense feature; the thinner the chemistry, the more the residual falls back onto raw photometric quantities to reconstruct what the baseline never computed. The question of whether the baseline analytical output or the bare turbidity proxy dominates a residual model therefore admits no general answer; it is itself a probe of how much analytical work the certified formula was doing before the AI layer was attached.

The physicochemical reading also dictates a chemistry-aware data-collection strategy that goes beyond standard ML-deployment heuristics. The granularity of the laboratory concentration grid should be selected to match the geometry of the deployment range relative to the analyzer’s calibrated dynamic range, not to fill the calibrated range uniformly. Where the analyzer is known to remain in the Beer-Lambert linear region throughout the operational range (the situation that applies to colorimetric methods such as the molybdenum-blue TP determination, and that produces the smooth, separable residual surface observed here), a sparse but uniform grid is adequate, because the residual structure is smooth enough for cross-concentration interpolation and the residual model can rely on the analytical formula's internal linearity. Where the analyzer is expected to be driven into detector non-linearity at the upper end of the operational range (the situation that applies to direct-UV ammonia determination and to any spectrophotometric channel operating near absorbance unity), the grid must include explicit over-range samples. The absence of such samples is not a sampling oversight but a structural blind spot: no quantity of in-range data can compensate, because the saturation behavior does not exist within the linear region and cannot be inferred from it.

This requirement reframes calibration-range extension as a chemistry-aware exercise rather than a generic ML one. A thin residual-learning layer can extend the effective working range of an existing instrument when, and only when, the deployment campaign includes samples that exercise the same detector non-linearity that the field instrument will encounter. The certified analytical formula remains untouched, preserving regulatory traceability and bypassing the months-long re-certification cycle required for changes to the analytical method itself [19]. The additive residual architecture's safety property anchors the learned correction to the unmodified analytical output, reducing the risk of unconstrained black-box behavior relative to direct concentration prediction from raw spectra [7], [25]. The qualifier, namely that tree-based residual models extrapolate poorly to concentrations entirely absent from the training distribution, whereas smoother parametric residuals extrapolate more gracefully, is itself a direct consequence of the bias-variance trade-off translated into the analytical chemistry of the measurement, and it implies that the choice of residual model family should follow the operational range expected at deployment rather than the in-grid accuracy alone. The mechanism is in principle transferable to other analytical techniques whose embedded formulas saturate outside their calibration range, including direct-reading sensors for pH, electrical conductivity, and dissolved oxygen, and other reagent-free spectrophotometric parameters such as COD and nitrate [23], [39], [40], [46], provided the same chemistry-aware data-collection discipline is observed.

The present study has four limitations that should temper the broader interpretation of its findings. (i) The NH$_4$-N and TP grids were collected on a single instrument and TN on only two; cross-instrument generalization has therefore not been assessed, and a model deployed on a third device would, on present evidence, require re-training or transfer-learning from the existing artefacts [25]. The TN comparison between the two instruments in the dataset already shows 15% MAE variation across folds attributable to inter-instrument photometric differences, suggesting that an instrument-agnostic residual model is not directly available without a calibration-transfer step. (ii) The controlled-turbidity samples used a formazine reference, whose scattering phase function differs from the natural particulates (clay, humic colloids, biofilm fragments) encountered in surface and drinking water. The TP field-sample evaluation in Section 3.2, in which every method including the unmodified baseline over-estimated the 0.10 mg·L$^{-1}$ spike by 0.08–0.10 mg·L$^{-1}$, is consistent with this matrix-effect concern [7], [42]. (iii) The concentration grids are sparse (four target levels per analyte), which makes the LOCO scheme simultaneously stringent and noisy. Held-out folds for NH$_4$-N differ from the training folds by more than a factor of two in concentration, an extrapolation gap larger than typical deployment conditions, whereas held-out folds for TP differ by only 0.05 mg·L$^{-1}$, an interpolation step that may flatter LOCO performance. (iv) All measurements were collected within a single experimental campaign of approximately three weeks, so long-term temporal drift (light-source ageing, optical-window fouling, electronics warm-up) has not been characterized, leaving the question of model durability over months of deployment open [25], [43].

Future work will address these limitations along four corresponding axes. Multi-instrument data collection, ideally across three or more nominally identical analyzers from at least two manufacturers, will quantify the inter-instrument generalization gap and permit explicit calibration-transfer or transfer-learning experiments. Real-water matrix studies using particulates representative of the deployment environment (kaolinite for drinking water, mixed-organic colloids for surface water, return-activated sludge for wastewater) will probe whether the residual-learning gains translate from formazine-spiked to natural turbidity. Denser concentration grids (six to ten target levels per analyte, with at least one above the embedded calibration band for each saturating analyte) will tighten the LOCO testing and provide enough information for tree-based in-grid saturation correction to generalize beyond the densest part of the calibration grid. Longitudinal drift quantification, with monthly data collection over at least one calendar year, will determine whether the residual model itself requires periodic re-fitting and, if so, on what cadence and with how much new data. Taken together, these four extensions would convert the in-grid residual-learning gains documented here into a deployable, certifiable AI-augmented analytical method for nutrient monitoring in real drinking-water and wastewater systems.

We have shown that a unified residual-learning AI framework, applied with the same skeleton across three nutrient analytes, substantially reduces the turbidity-induced error of existing on-line spectrophotometric water quality sensors. The production Ridge regression model delivers MAE reductions of 42.5% (TN), 75.6% (TP), and 7.6% (NH$_4$-N) under random 5-fold cross-validation. The XGBoost residual variant achieves further gains where appropriate (45.2%, 87.2%, and 75.6%, respectively), and corrects the observed calibration-range saturation in the NH$_4$-N baseline at 2.35 mg·L$^{-1}$ (97.1% MAE reduction). The residual-learning architecture's safety properties (additive on the certified analytical formula, with $\hat{\Delta}$→0 at zero turbidity) make on-site deployment robust to model failure. LOCO validation reveals real generalisation limits, particularly for NH$_4$-N, that are invisible to conventional random-fold reporting; we recommend that both schemes be reported as a baseline standard in this domain. Future work should test cross-instrument transferability, extend the data collection to denser concentration grids, and quantify long-term drift to enable autonomous recalibration.