npj Climate and Atmospheric Science Article · (2026) 9:157 Open Access · CC BY-NC-ND 4.0

Beyond satellite-based precipitation data: a novel soil-moisture physics framework with Green–Ampt and Bayesian optimization for rainfall estimation

Mohammad Saeedi1, Hyunglok Kim2✉, John Bolten3, John Eylander4, Venkataraman Lakshmi1
1Civil & Environmental Engineering, University of Virginia · 2Environment & Energy Engineering, GIST, Gwangju, South Korea · 3Hydrological Sciences Laboratory, NASA GSFC · 4Coastal & Hydraulics Laboratory, US Army ERDC
https://doi.org/10.1038/s41612-026-01425-5 · Received 5 Feb 2026 · Accepted 23 Apr 2026 · Published 7 May 2026
The thesis in one line
Replace SM2RAIN's empirical power-law drainage with a physically based Green–Ampt wetting-front infiltration term, then test how far down the soil column the rainfall signal survives — 0–20 cm, four layers, globally.
Best layer
0–5 cm
GA wins at 54% of stations
GA vs GPM · FAR
+50%
false alarms cut, 0–5 cm
GA vs GPM · NSE
+48%
continuous skill, 0–5 cm
01The gap why land rainfall constraints are still weak

Gauges are sparse and uneven

Networks are difficult to sustain at scale, radar coverage is geographically limited and needs calibration — heterogeneous, data-sparse regions stay imperfectly constrained.

Satellite products carry regime-dependent bias

Sampling, retrieval physics, surface and precipitation-type heterogeneity propagate into climate diagnostics and every workflow that uses precipitation as reference or forcing.

SM2RAIN's two bottlenecks

Inverting the soil-moisture water balance works — but skill is governed by (i) how infiltration/drainage physics is represented and (ii) the depth at which soil moisture is observed. Existing implementations fix literature parameters and use shallow layers only.

02Data & experimental design in-situ truth, four depths, one satellite benchmark

International Soil Moisture Network (ISMN)

Period2010 – 2023 Record filter≥ 5 years of daily data per station Native stephourly SM + precipitation, harmonised Referenceco-located in-situ rain gauges BenchmarkGPM IMERG (gridded satellite precipitation) SoftwareMATLAB R2023b · maps via Python Cartopy

Fig. 1 of the paper maps global station coverage, colour-coded by how many of the four soil-moisture depths each station reports (1→4). Coverage is dense over North America and Europe, sparse elsewhere.

Stratifications tested

Depth × 4Land cover × 5Köppen climate × 8

The depth ladder — where the rainfall imprint lives

The rainfall imprint in soil moisture is sharpest at the surface and becomes progressively attenuated, delayed and dispersed with depth — this is the paper's central control variable.

03The physics: from water balance to two inversions "soil as a natural rain gauge"

Common core — the SM2RAIN water balance (Brocca et al.)

EQ 1 Z(t)dt = P(t) − R(t) − ET(t) − G(t)
EQ 2–4 R(t) = P(t) θ(t)c   ·   G(t) = aθ(t)b   ·   ET(t) = ETpotθ(t)
EQ 5 — classic inversion P(t) = Z(t)dt + aθ(t)b 1 − θ(t)c

ET is assumed negligible during rainfall. Rainfall is inferred from the competition between storage tendency and loss terms — so errors in the drainage / infiltration term systematically bias retrieved rainfall and can make the inverse problem ill-conditioned.

Zdepth-related soil parameter (soil layer depth)
θ(t)relative / volumetric soil moisture
a, b, cempirical drainage & runoff exponents (calibrated)

Branch A SM2RAIN–NWF · analytical Net Water Flux

Sadeghi et al.'s NWF model, embedded in SM2RAIN by Saeedi et al. It reconstructs the effective vertical water flux from SM by analytically inverting the 1-D linearized Richards equation.

EQ 6 — recursive flux FN = { Θ1(Z) e0.5Z−0.25TU(Z, ΔT),N = 1 FN−1 + ΘN(Z)e0.5Z−0.25TF1U(Z,T) − Σi=3N[Fi−1Fi−2]U(Z,(Ni+2)ΔT),N ≥ 2
EQ 7 — modified inversion P(t) = Z(t)dt + a f(t) 1 − θ(t)c
ΘN(Z)normalised SM anomaly, (θ(t) − θ)/θ; θ = background moisture
U(Z,T)analytical unit response function
Tdimensionless time, T = k²t/Di (k = hydraulic conductivity, Di = diffusivity)
f(t)drainage/flux term G(t) proxy inserted into the inversion
fmax, Di, kf, zfinfiltration parameters — fixed from literature in the original scheme

Limitation attacked: fixed literature parameters neglect site-specific soil-hydraulic variability → systematic bias.

NEW · Branch B SM2RAIN–GA · Green–Ampt wetting front

First application of a Green–Ampt infiltration calculation inside a bottom-up rainfall-inversion framework. A sharp wetting front separates a saturated zone (θs) above from the initial water content (θi) below and advances downward.

EQ 8 — infiltration rate dFdt = Ks(1 + Ψf ΔθF(t) )
EQ 9 — implicit cumulative form F(t) − F0ΨfΔθ ln ( F(t) + ΨfΔθF0 + ΨfΔθ ) = Ks(tt0)
Instantaneous rate f(t) = dF/dt from Eq. 8 replaces the NWF flux term inside Eq. 7.
F, F0cumulative infiltration depth at t and t0 (mm)
Kssaturated hydraulic conductivity — constrained to [0.1, 200] mm h⁻¹
Ψfwetting-front suction head — constrained to [20, 1000] mm
Δθθs − θi, constant moisture deficit across the front

Assumption & its cost: 1-D vertical flow with a sharp front is most defensible in shallow layers. Deeper, diffuse redistribution, lateral flow, root uptake, preferential flow and heterogeneity take over — so at 15–20 cm the GA term must be read as an effective first-order, physically constrained representation, not full mechanistic subsurface flow.

04Learning the parameters & probing their uncertainty Bayesian optimization + one-at-a-time stochastic sensitivity

SM2RAIN–NWFCalib — site-wise Bayesian optimization

The fixed-literature-parameter assumption is relaxed: all infiltration- and SM2RAIN-specific parameters are estimated independently for every station.

EQ 10–11 θ* = arg minθ∈𝒳 L(θ)   ,   L(θ) = −NSE(Psim, Pobs)
EQ 12 — GP surrogate L(θ) ~ 𝒢𝒫(m(θ), k(θ, θ′))
EQ 13 — Expected Improvement α(θ) = ξ + 𝔼[max(0, LbestL(θ))]

ξ > 0 encourages exploratory sampling in high-posterior-uncertainty regions — balancing exploration against exploitation with a limited budget of expensive model evaluations. All infiltration- and SM2RAIN-specific parameters are jointly calibrated; the four reported against elevation in Fig. 8 are fmax, Di, kf, zf.

One-at-a-time (OAT) stochastic sensitivity

Isolates the marginal effect of each calibrated parameter: perturb one, hold the rest at their BO optima, keep the infiltration term fixed. 1000 realizations × 8 representative stations spanning a broad elevation and geographic range.

EQ 14–16 — three perturbation blocks B1: a(i) = |a* + a*ξi|,  c(i)=c*,  Z(i)=Z* B2: a(i)=a*,  c(i) = |c* + c*ηi|,  Z(i)=Z* B3: a(i)=a*,  c(i)=c*,  Z(i) = |Z* + Z*ζi| ξi, ηi, ζi ~ i.i.d. 𝒩(0,1)  ·  perturbation applied as P = Popt + σP·ε
EQ 17 — ensemble statistics μP(t) = 1NΣi=1NP(i)sim(t)  ,  σP(t) = √ 1N−1Σ(P(i)simμP

Reported per time step: mean, σ, 2.5th/97.5th percentiles (95% CI) and the 25–75% IQR. Result: at almost all stations the deterministic BO estimates for Z and a fall inside the 5–95% Monte-Carlo envelope and near the distribution median — the calibration is robust to parameter perturbation.

05Evaluation protocol three products · five metrics · four depths · non-parametric testing

Products compared

SM2RAIN–NWF bottom-up, analytical flux SM2RAIN–GA bottom-up, Green–Ampt GPM IMERG top-down satellite benchmark

Continuous metrics

  • NSE — Nash–Sutcliffe efficiency (1 = perfect; <0 worse than the observed mean)
  • R — Pearson correlation
  • RMSE — root mean square error

Categorical (event) metrics

  • FAR — false alarm ratio (lower = better)
  • POD — probability of detection (higher = better)

Together they separate "getting the totals right" from "getting the storms right".

Statistics

  • Paired Wilcoxon signed-rank tests on station-wise metrics — chosen because pairing is natural at the station and distributions are skewed with outliers
  • Benjamini–Hochberg FDR correction; adjusted q-values reported
  • Rank-based composite score across all five metrics → per-station "winner"
06Results every panel below is interactive — hover any mark

Percent change relative to GPM, by soil-moisture depth (Fig. 3)

SM2RAIN–NWF SM2RAIN–GA positive = better than GPM

Table view — all metrics, all depths

Which product wins at each station? (Fig. 4)

Share of stations where each method ranks first on the joint rank-based score over all five metrics.

NWF GA GPM

GA dominates at shallow depth, but the winner distribution shifts systematically downward: GA 54%→43%, NWF 30%→20%, while GPM climbs 16%→37%. Consistent with attenuation of the rainfall imprint in deeper soil.

Significance of the inter-model comparison (Fig. 5)

Paired Wilcoxon signed-rank on station-wise metrics, FDR-adjusted q-values.

SM2RAIN–NWF vs GPM
SM2RAIN–GA vs GPM
SM2RAIN significantly better GPM significantly better Non-significant Significant, median paired difference ≈ 0

Read this way: FAR favours the bottom-up schemes at every depth for both methods — all eight cells green. NSE favours them everywhere except GA at 15–20 cm, where the difference is non-significant (q = 0.199). NWF loses correlation and POD in the two deepest layers; GA never loses POD anywhere — significantly better at all four depths.

Skill decay with depth — the headline pattern

Composite view: GA's advantage over GPM on event detection (POD ↑ and FAR ↓ simultaneously) is what separates it from NWF, and it survives deeper into the column.

Stratified by land cover, 0–5 cm (Fig. 6, Table 1)

NWF vs GPM GA vs GPM

Classes: 10 cropland (rainfed) · 70 tree cover, needle-leaved evergreen, closed-to-open (>15%) · 120 shrubland · 130 grassland · 190 urban. Mean improvement across classes — NSE +65% (NWF) / +53% (GA); FAR +20% vs +46%; POD +14% vs +49%. Replacing the drainage term with Green–Ampt is what buys transferability across surfaces.

Table view — land cover

Stratified by Köppen climate, 0–5 cm (Fig. 7)

NWF vs GPM GA vs GPM

Mean over classes — NSE +78% / +69%; R +14% / +17%; RMSE +29% / +24%; FAR +25% / +53%; POD +21% / +48%. Sign reversals matter: under NWF, POD collapses in Cfb (−55.3%) and Cfa (−7.6%); GA turns both positive (+11.7%, +10.5%). The one systematic exception is Cfa, where GA's RMSE change is slightly negative despite better FAR and POD — residual mismatch in infiltration timing.

Table view — climate

Fixed literature parameters vs site-wise BO (Fig. 8)

Optimized fmax, Di, kf and zf span broad ranges with no monotonic dependence on elevation — elevation is not a primary control, and much of what BO captures is idiosyncratic rather than a large-scale driver.

Central tendency and spread of R and RMSE largely overlap between the two configurations; FAR changes are minor; the one systematic deviation is POD, which is higher under the fixed setting. Holding (fmax, Di, kf, zf) fixed does not materially degrade skill for most stations — and it removes the computational burden of site-wise BO.

OAT sensitivity of Z and a (Fig. 9)

1000 stochastic realizations per station, 8 stations. Box = 25–75% IQR, whiskers = 5th–95th percentile, dashed = mean ± 1σ, green dot = the deterministic BO estimate.

At almost all stations the deterministic calibration sits inside the 5–95% uncertainty envelope and close to the median of the distribution — the calibrated coefficients are relatively robust to the imposed perturbations. (Schematic reproduction of the paper's Fig. 9 layout; box widths/whiskers illustrate the reported structure.)

07The fairness check: point gauge vs satellite grid how much of GPM's penalty is a sampling artefact?

SM2RAIN estimates are validated against point-scale rain gauges while GPM IMERG represents areal rainfall at grid scale. To quantify that penalty the authors ran a multi-gauge pixel sensitivity analysis on GPM pixels containing several gauges (Fig. S1).

Median P95
−2.99%
reduced by areal reference
Median P99
−6.38%
reduced by areal reference
GPM false alarms
15.9%
resolved by areal reference

Context from the literature

  • Standard error from a single gauge ≈ 33% for a 2 × 2 km pixel
  • Gauge-to-gauge differences within a 500 × 500 m pixel can reach 100% over several days
  • > 25 gauges may be needed to push area-average error below 20% at fine temporal scale
  • Sparse gauge networks tend to underestimate IMERG skill relative to dense-gauge references

Verdict: scale mismatch explains only part of the observed performance difference — it does not overturn the main findings.

08What to take away three contributions · three firm results · operational guidance

Three contributions

  • i. First integration of a Green–Ampt infiltration constraint into an SM2RAIN-type inversion — an interpretable route to reducing regime-dependent drainage/infiltration error.
  • ii. First depth-stratified, multi-layer assessment of SM2RAIN-family retrieval using in-situ SM across 0–20 cm — delineating where and why skill weakens with depth.
  • iii. Site-specific parameters via Bayesian optimization plus OAT stochastic experiments that quantify sensitivity to parametric uncertainty.

Three firm results

  • Depth matters. Performance peaks at 0–5 cm and degrades progressively as the rainfall imprint attenuates.
  • Green–Ampt > NWF flux closure. GA improves on both the NWF variant and GPM across NSE, R, RMSE, POD and FAR — with simultaneous POD gains and FAR reductions, i.e. real event discrimination, not just better totals.
  • Fixed parameters are a defensible default. BO's benefit is regime-dependent (mostly FAR in Cfa, Cfb, Dwb, Dfb), not universal.

Operational recommendations

  • Prioritise near-surface layers for inversion.
  • Adopt GA-based SM2RAIN as the default bottom-up configuration.
  • Tailor parameters by depth, land cover and climate.
  • For continental/global roll-out: start from the fixed parameter set (light, stable POD); reserve BO for sites whose post-hoc diagnostics show persistent volumetric bias or regime mismatch.
  • Ship OAT uncertainty alongside operational products.

Caveats the authors keep on the table

  • Point-scale rainfall and profile SM do not sample the same footprint — affects both correlation and event metrics.
  • Unmodelled processes decouple rainfall from the SM signal: vegetation, frozen soil/snow, irrigation, soil-moisture noise → inflated false detections or misses.
  • The GA sharp-front assumption weakens with depth; Bayesian optimization can only partially absorb the unresolved site-specific complexity.
  • OAT perturbs selected parameters individually; extending uncertainty treatment to observational and structural components would strengthen reliability claims.
  • Class-dependent trade-offs persist, most visibly in Cfa (humid subtropical) and Cfb (oceanic).

Where this goes next

  • Fuse multi-depth information rather than choosing one layer.
  • Relax the simplifying assumptions of the front-based formulation.
  • Assimilate additional hydrologic cues.
  • Extend benchmarking beyond GPM to multi-sensor products.

Why it matters for climate science

More defensible, observation-based rainfall constraints — with explicit uncertainty — for climate diagnostics, land–atmosphere coupling analysis, and benchmarking of gridded precipitation products in exactly the heterogeneous, data-scarce regions where gauges fail.