1.1 Background

Machine learning models are increasingly used in hydro-climate science to predict streamflow, estimate drought indices, map flood risk, and project evapotranspiration under changing climate scenarios. These models — random forests, gradient boosting machines, deep neural networks — often outperform traditional physically-based models in predictive accuracy (Breiman 2001). However, they operate as black boxes: their internal decision logic is opaque to practitioners, regulators, and policy makers who must trust model outputs for operational water management and climate adaptation decisions.

Explainable Artificial Intelligence (XAI) addresses this gap by providing mathematically rigorous tools to attribute model predictions to individual input features. XAI transforms a black-box model into an interpretable system without sacrificing predictive performance.

1.2 The xaiHydro package

xaiHydro provides a unified, hydrology-aware R workflow for applying three complementary XAI families to any trained hydro-climate model:

All three methods require only a trained model and its training data — no additional field data collection, no satellite downloads, no pre-processing pipelines. This makes xaiHydro uniquely accessible for researchers who already have a fitted model and want to explain it.

1.3 Scope of this vignette

This vignette serves as both:

1. A package tutorial — covering every function in xaiHydro with reproducible code.
2. The book chapter case study — the outputs produced here are the figures and analyses referenced in the accompanying Springer chapter.

3.1 Daily streamflow dataset

df_flow <- sim_streamflow_data(n = 365, seed = 2026)
head(df_flow, 6)

str(df_flow)
#> 'data.frame':    365 obs. of  9 variables:
#>  $ precipitation     : num  9.22 0 6.27 4.53 0 0 0 0 6.37 1.72 ...
#>  $ temperature       : num  23.8 17 20.7 23.2 22.8 ...
#>  $ evapotranspiration: num  5.49 4.85 5.56 6.31 5.72 ...
#>  $ soil_moisture     : num  0.34 0.307 0.264 0.227 0.283 ...
#>  $ snow_water_equiv  : num  7.17 3.83 2.07 11.25 5.94 ...
#>  $ ndvi              : num  0.496 0.505 0.468 0.528 0.459 ...
#>  $ lag1_precip       : num  0 9.22 0 6.27 4.53 0 0 0 0 6.37 ...
#>  $ lag3_precip       : num  0 0 5.16 3.6 3.6 1.51 0 0 2.12 2.7 ...
#>  $ streamflow        : num  11.1 12.9 10.5 15.7 11.8 ...

The function generates a synthetic daily hydro-meteorological dataset mimicking realistic seasonal dynamics:

The response is generated as a non-linear function of predictors with added noise, including a precipitation × NDVI interaction term to capture vegetation-modulated runoff.

3.2 Monthly drought index dataset

df_drought <- sim_drought_data(n = 240, seed = 2026)
head(df_drought, 6)

This dataset simulates 20 years of monthly hydro-climatic anomalies with SPEI (Standardised Precipitation-Evapotranspiration Index) as the response, driven by precipitation anomalies, temperature anomalies, SPI at 3- and 12-month scales, PDSI, NDVI anomaly, and ET deficit.

4.1 Random Forest — streamflow prediction

We train a Random Forest model on the simulated streamflow dataset. The same workflow applies to any regression model supported by predict().

library(randomForest)

# Separate predictors and response
y_flow <- df_flow$streamflow
X_flow <- df_flow[, setdiff(names(df_flow), "streamflow")]

set.seed(42)
rf_model <- randomForest(
  x         = X_flow,
  y         = y_flow,
  ntree     = 300,
  mtry      = 3,
  importance = TRUE
)

print(rf_model)
# Mean of squared residuals: ~X.XX
# % Var explained:            ~XX.X

4.2 XGBoost — drought index prediction

library(xgboost)

y_spei   <- df_drought$spei
X_drought <- df_drought[, setdiff(names(df_drought), "spei")]

xgb_mat   <- xgboost::xgb.DMatrix(as.matrix(X_drought), label = y_spei)
xgb_model <- xgboost::xgboost(
  data      = xgb_mat,
  nrounds   = 150,
  max_depth = 4,
  eta       = 0.1,
  objective = "reg:squarederror",
  verbose   = 0
)

5.1 Creating an explainer

hydro_explainer() is the single entry point for all XAI analyses. It wraps [DALEX::explain()] (Biecek and Burzykowski 2021) and attaches hydrology-specific metadata used for axis labels, titles, and interpretation messages throughout the package.

exp_flow <- hydro_explainer(
  model    = rf_model,
  data     = X_flow,
  y        = y_flow,
  variable = "streamflow",
  units    = "m3/s",
  label    = "RandomForest"
)

print(exp_flow)
# == xaiHydro Explainer ==
#   Target variable : streamflow [m3/s]
#   Model label     : RandomForest
#   Observations    : 365
#   Predictors      : 8
#   Model class     : randomForest

The explainer object is compatible with all xaiHydro functions and can also be passed directly to any DALEX function for additional analyses.

5.2 Explainer for XGBoost (custom predict wrapper)

For models that require a non-standard prediction call (XGBoost, Keras, PyTorch via torch), supply a predict_function:

pfun_xgb <- function(model, newdata) {
  predict(model, xgboost::xgb.DMatrix(as.matrix(newdata)))
}

exp_drought <- hydro_explainer(
  model            = xgb_model,
  data             = X_drought,
  y                = y_spei,
  variable         = "drought_index",
  units            = "SPEI",
  label            = "XGBoost",
  predict_function = pfun_xgb
)

6.1 Computing SHAP values

shap_vals <- hydro_shap(
  explainer = exp_flow,
  nsim      = 50,     # Monte Carlo replicates — increase for stability
  seed      = 42
)

# Returns a tidy data frame:
# id | feature | shap_value | feature_value | mean_abs_shap
head(shap_vals, 8)

Tip: nsim = 50 is adequate for publication figures. Use nsim = 100–200 for final manuscript outputs.

6.2 Global SHAP summary plot

plot_shap_summary(shap_vals, top_n = 8)

The beeswarm summary plot shows:

• Y-axis — features ranked by mean |SHAP| (most important at top).
• X-axis — SHAP value (positive = pushes prediction up; negative = pulls it down).
• Colour — raw feature value (red = high, blue = low).

Features where red points cluster on the right indicate that high values of that predictor increase predicted streamflow — physically consistent with rainfall–runoff theory.

6.3 Single-observation waterfall

# Explain a high-flow day (e.g., observation 180 after a storm event)
plot_shap_waterfall(shap_vals, obs_id = 180, top_n = 7)

The waterfall plot shows exactly how prediction for day 180 was built up from the global baseline (mean predicted streamflow) through additive feature contributions. Red bars push the prediction above baseline; blue bars pull it below.

8.1 PDP — marginal effects

pdp_result <- hydro_pdp(
  explainer   = exp_flow,
  variable    = c("precipitation", "soil_moisture",
                   "lag1_precip",  "temperature"),
  type        = "partial",
  grid_points = 60
)

plot_pdp_hydro(pdp_result, ncol = 2)

8.2 ALE — preferred for correlated predictors

ale_result <- hydro_pdp(
  explainer = exp_flow,
  variable  = c("precipitation", "soil_moisture"),
  type      = "accumulated"
)

plot_pdp_hydro(ale_result, ncol = 2)

In hydrology, precipitation and soil moisture are typically positively correlated (wet soils follow wet periods). ALE correctly handles this by computing effects within observed data slices rather than averaging over impossible combinations.

14.1 Physical interpretation of SHAP

When precipitation and lag1_precip appear in the top-ranked SHAP features with positive effects (red points on the right), this confirms the model has correctly learnt that higher antecedent and current rainfall drives higher streamflow — consistent with rainfall–runoff theory. Conversely, high evapotranspiration having negative SHAP (blue points on left) indicates that high ET demand reduces net water available for runoff.

14.2 ALE vs PDP selection guide

Use ALE when: - Predictors are correlated (precipitation–soil moisture, temperature–ET). - The training data does not uniformly cover the predictor space. - You suspect extrapolation artefacts in PDP.

Use PDP when: - Predictors are approximately independent. - You want a simpler figure for a non-technical audience.

14.3 LIME local R² threshold

An R² < 0.5 for the LIME surrogate means the local linear approximation is unreliable for that observation. Increase kernel_width (to use a broader neighbourhood) or reduce n_features (for a sparser model).