Evaluating Machine Learning Models for Post-Wildfire Debris-Flow Prediction

# Evaluating Machine Learning Models for Post-Wildfire Debris-Flow Prediction
Source: [https://arxiv.org/html/2608.05265](https://arxiv.org/html/2608.05265)
1\]organization=Department of Geomatics Engineering, University of Calgary, city=Calgary, state=AB, country=Canada

\\cortext

\[1\]Corresponding author

Zhengsen XuYimin ZhuZack DewisMabel HeffringSaeid TaleghanidoozdoozanMotasem AlkayidMegan GreenwoodLincoln Linlin Xulincoln\.xu@ucalgary\.ca\[

###### Abstract

Prediction of post\-wildfire debris flows is critical for mitigating hazards to communities, infrastructure, and resources during intense rainfall in recently burned areas\. However, identifying reliable machine learning models for this task is complicated by several factors, including the overlap of debris\-flow and no\-debris\-flow events in feature space, the lack of physical interpretability, and limited training data\. Therefore, the goal of this paper is to address these challenges through a systematic evaluation of a broad set of machine learning models in terms of model performance, feature importance, and response to synthetic data augmentation\. Using basin\-scale observations of post\-wildfire debris\-flow events across the western United States, we compare 15 models, including a new foundation model, the Tabular Prior\-Data Fitted Network \(TabPFN\)\. The results of repeated stratified cross\-validation indicate that TabPFN achieves the strongest unaugmented performance with a threat score of 0\.637, closely followed by the leading tree\-based models\. We use SHapley Additive exPlanations \(SHAP\) for the feature importance evaluation to identify which features the top\-performing models rely on most when predicting debris flows\. This evaluation finds that short\-duration rainfall intensity and storm accumulation receive the highest SHAP rankings across the explained models, with burn severity and terrain features ranked lower\. Finally, we quantify the utility of synthetic data augmentation to mitigate the scarcity of debris\-flow event observations\. Using TabPFN to generate synthetic training data, we improve the performance of all models except CNN, with the largest mean gain in threat score of\+0\.041\+0\.041among the deep\-learning models\. By integrating rigorous model benchmarking, transparent feature evaluation, and data augmentation, this work provides a comprehensive framework for enhancing the accuracy and reliability of post\-wildfire debris\-flow prediction\.

###### keywords:

Post\-wildfire debris\-flow prediction\\sepMachine learning\\sepTabular foundation model\\sepTabPFN\\sepSHAP\\sepSynthetic data augmentation

## 1Introduction

![Refer to caption](https://arxiv.org/html/2608.05265v1/x1.png)Figure 1:Study\-area overview of the 34 fire events in the USGS post\-wildfire debris\-flow dataset \(2000–2013, western United States\)\. Panel \(a\) shows the regional extent with labeled fire centroids on a digital elevation model \(DEM\), panel \(b\) zooms to the southern California cluster, and panel \(c\) shows a United States locator\. Panels \(d\)–\(h\) show delineated basins for the Station, Grand Prix–Old, Schultz, Little Bear, and Waldo Canyon fires, with polygons filled to show the response \(no\-debris\-flow or debris flow\)\.Post\-wildfire debris flows threaten lives, infrastructure, and resources in burned watersheds across fire\-prone regions worldwide, and operational warning depends on models that estimate the probability of initiation from basin and storm features\(staley\_prediction\_2017\)\. The problem of predicting this hazard is difficult for three reasons\. First, selecting the best model for this task is challenging\. Many candidate models are available, and debris\-flow and non\-event basins overlap in their rainfall, burn\-severity, and terrain signatures \(Table[2](https://arxiv.org/html/2608.05265#S2.T2); Figure[8](https://arxiv.org/html/2608.05265#S3.F8)\), so no single threshold cleanly separates the two classes of outcomes\(nikolopoulos\_evaluation\_2018\), making the choice of model consequential\. Second, when classes overlap in feature space, different model families can place decision boundaries differently, so identifying which features a high\-performing model actually relies on is non\-trivial\. This feature importance is needed to support model\-selection and trust judgments, to diagnose whether predictions reflect robust signal rather than spurious patterns confined to a particular subset of observations, and to assess whether a model would generalize\(ribeiro\_why\_2016\)\. Third, the available data for this problem are typically small and class\-imbalanced, with debris\-flow events the minority outcome, which constrains how well any model can learn the boundary between the two outcomes and disproportionately limits identification of the rare positive class\(he\_learning\_2009\)\. Together, these difficulties motivate a rigorous, multi\-model evaluation under controlled conditions, similar to the strategies used in other specialized remote\-sensing prediction studies\(xu\_comparative\_2014\)\.

A typical post\-wildfire debris\-flow hazard assessment workflow proceeds from delineation of burned basins to assembly of rainfall and terrain features as a tabular collection of features to prediction of debris flows\(staley\_updated\_2016\)\. Within this workflow, the prediction stage is critical because missed events leave exposed communities without warning, whereas false alarms increase operational burden\(sattele\_reliability\_2015\)\. This paper, therefore, focuses on the prediction stage, where input features are translated into warning decisions\. Accordingly, we evaluate models using event\-detection metrics, such as threat score, that balance missed events and false alarms under class imbalance, rather than overall accuracy alone\(schaefer\_critical\_1990\)\.

Several machine learning models have been applied to post\-wildfire debris\-flow prediction\. Logistic regression models that relate rainfall intensity and duration to initiation probability underlie current operational approaches by the United States Geological Survey \(USGS\)\(staley\_objective\_2013;staley\_updated\_2016;staley\_prediction\_2017\)\. Based on the USGS approach, improved threat scores have been reported for naïve Bayes, mixture discriminant analysis, and classification trees\(kern\_machine\_2017\), decision trees\(addison\_assessment\_2019\), eXtreme Gradient Boosting \(XGBoost\) trained on satellite\-based inputs\(orland\_scalable\_2022\), and Random Forest and neural network models\(nikolopoulos\_evaluation\_2018;roten\_machine\_2022\)\. Beyond the United States, logistic regression models have been applied to post\-wildfire debris\-flow prediction in the Mediterranean and southwest China\(diakakis\_exploring\_2023;jin\_susceptibility\_2022\)\. Three gaps remain in this body of work\. First, prior studies typically compare only a handful of models and rely on a single train\-test split\. This makes performance estimates highly variable and heavily dependent on the specific split, particularly when performance gaps are small\(cawley\_over\-fitting\_2010;krstajic\_cross\-validation\_2014;dietterich\_approximate\_1998\)\. In addition, several high\-performing recently developed models for tabular data remain unevaluated for this problem\. Tabular Prior\-Data Fitted Network \(TabPFN\) is particularly relevant, as it performs approximate Bayesian prediction via in\-context learning to excel on small\-to\-medium datasets\(hollmann\_tabpfn\_2023;hollmann\_accurate\_2025\)\. Beyond the debris\-flow setting, TabPFN has outperformed conventional models in thaw\-hazard mapping\(zhu\_method\_2026\)and has shown strong performance for landslide susceptibility mapping under data\-limited conditions\(zhou\_landslide\_2025;yang\_knowledge\-data\_2026\)\. Similarly, the tree\-structured Recursive Feature Machine \(xRFM\) provides a scalable approach to feature learning that has shown strong performance on tabular benchmarks\(beaglehole\_xrfm\_2025\)\. Second, the features that high\-performing models for post\-wildfire debris\-flow prediction rely on are not well characterized\. A meta\-analysis of debris\-flow machine learning studies confirmed that interpretation methods are rarely applied across this literature\(yang\_machine\-learning\-based\_2024\)\. Relatedly, whether different model families attend to the same features or to divergent signals remains an open question\. Third, while synthetic data augmentation has improved related hazard models under similar class imbalance\(wang\_optimizing\_2019;kumar\_addressing\_2024\), whether such augmentation transfers to the post\-wildfire debris\-flow setting, and which models benefit when it is added, is not well understood\.

Therefore, this paper presents a rigorous evaluation framework for post\-wildfire debris\-flow prediction\. We conduct a systematic comparison of fifteen models across five model families: logistic, distance\- and kernel\-based, deep\-learning, tree\-based, and tabular foundation models\. Each model is trained under a shared preprocessing pipeline and a common hyperparameter\-tuning framework where applicable, then evaluated with repeated stratifiedkk\-fold cross\-validation\. This isolates model architecture as the primary variable and averages performance over many resamplings, so observed gaps reflect model choice rather than uneven configuration effort or sampling noise from a single split\. Beyond model ranking, we use SHapley Additive exPlanations \(SHAP\) feature importances\(lundberg\_unified\_2017\)both to characterize the features that individual high\-performing models rely on and to contrast importance patterns across model families on shared splits, asking whether high\-ranked models converge on the same features or rely on divergent signals\. Finally, we generate synthetic training observations with TabPFN and quantify, under the same evaluation protocol, whether augmentation improves prediction in this class\-imbalanced setting and which of the fifteen models possibly benefit from it\. Together, these results clarify which model families best handle the overlap between debris\-flow and non\-event basins, identify which features high\-performing models rely on and where families agree or diverge, and quantify whether synthetic augmentation can mitigate the small, class\-imbalanced training data typical of this task\.

## 2Methodology

Data and fold\-wise pipeline[Data Sources](https://arxiv.org/html/2608.05265#S2.F4)– USGS post\-wildfire debris\-flow dataset – 1,550 observations; 33 fires; 248 storms – Seven western U\.S\. states \(2000–2013\)[Features and Response](https://arxiv.org/html/2608.05265#S2.SS1)– Ten features \(Table[1](https://arxiv.org/html/2608.05265#S2.T1)\) – Six rainfall, four basin features – Binary response \(debris flow: 0/1\)[Cross\-Validation](https://arxiv.org/html/2608.05265#S2.SS2)– Repeated stratified cross\-validation – 10×\\times5 splits; identical folds across models – Outer test sets contain real observations only[Preprocessing of Features](https://arxiv.org/html/2608.05265#S2.SS3)– Fold\-fitted imputation and standardization – Monotonic rainfall constraints; median gap fill[Hyperparameter Tuning](https://arxiv.org/html/2608.05265#S2.SS4)– Threat\-score\-optimized search \(100 trials\) – Separate tuning for unaugmented and augmented\(1\) Model Comparison[Models Compared](https://arxiv.org/html/2608.05265#S2.SS5)– 15 models, 5 families:[Logistic Regression](https://arxiv.org/html/2608.05265#S2.SS6),[Distance\- and Kernel\-based](https://arxiv.org/html/2608.05265#S2.SS7),[Deep\-learning](https://arxiv.org/html/2608.05265#S2.SS8),[Tree\-based](https://arxiv.org/html/2608.05265#S2.SS9), and[TabPFN](https://arxiv.org/html/2608.05265#S2.SS10) – Unaugmented vs\. Augmented training [Evaluation Metrics](https://arxiv.org/html/2608.05265#S2.SS11) – Threat score \(primary\), AUC, F1[\(2\) SHAP Feature Importance Evaluation](https://arxiv.org/html/2608.05265#S2.SS12)SHAP Computation: – Permutation SHAP \(TabPFN\); TreeSHAP \(Random Forest, XGBoost, CatBoost\) – Per fold: refit on train, explain held\-out rows Importance Evaluation: – Beeswarm, fold\-stability, cross\-model ranking, feature\-dependence plots – Basin\-level SHAP importance maps\(3\) Synthetic Data Augmentation Evaluation[Synthetic Data Generation:](https://arxiv.org/html/2608.05265#S2.SS13)– Class\-conditional TabPFN generator per fold – Fidelity \(KS, Wasserstein\-1, energy distance\) [Augmentation Evaluation:](https://arxiv.org/html/2608.05265#S2.SS14) – All 15 models share identical synthetic data – Real rows duplicated then synthetics appended – Primary Metric: Per\-model threat scoreΔ\\Delta

Figure 2:Conceptual overview of the experimental pipeline\. The left column lists the sequential steps applied within each outer training fold: data sources, features and response, cross\-validation splitting, preprocessing, and hyperparameter tuning\. The right column shows the three analysis modules: the 15\-model benchmark, SHAP feature importance evaluation, and TabPFN\-based synthetic data augmentation evaluation\.Figure[2](https://arxiv.org/html/2608.05265#S2.F2)summarizes the overall methodology\. The left column covers data compilation, features, cross\-validation, preprocessing, and hyperparameter tuning\. The right column maps to the three primary contributions: the 15\-model benchmark, SHAP feature importance, and TabPFN\-driven synthetic data augmentation\. This section follows the same order\. We first describe the dataset and feature set, then the repeated stratified cross\-validation design, within\-fold preprocessing, and hyperparameter tuning\. We then present the model families, evaluation metrics, paired model comparisons, SHAP feature importance, and the TabPFN\-based synthetic data augmentation procedure\.

### 2\.1Dataset

Table 1:Features included in the model comparison\. Features are grouped by category, with feature names, units, and brief descriptions\.CategoryFeatureUnitDescriptionFire & TerraindNBR\-Average differenced normalized burn ratio of watershedPropHM23\-Proportion of watershed burned at high/mod severity, gradient\>23∘\>23^\{\\circ\}ContribAreak​m2km^\{2\}Contributing area of observation locationSoilKF\-Average KF\-Factor \(erodibility index of the fine fragments of the soil\) of the watershedStormStormAccummmTotal rainfall accumulation of stormStormDurhTotal duration of stormStormAvgImm/hAverage storm intensityRainfallPeakI15mm/hPeak 15\-minute rainfall intensityPeakI30mm/hPeak 30\-minute rainfall intensityPeakI60mm/hPeak 60\-minute rainfall intensityTable 2:Summary statistics of post\-wildfire debris\-flow predictors by response class\. For each feature, the table reports minimum, mean, median, maximum, and interquartile range \(IQR\)\.FeatureNo debris flowDebris flowMinMeanMedMaxIQRMinMeanMedMaxIQRStormDur0\.0017\.8412\.4765\.0023\.350\.0020\.6423\.1763\.5029\.75StormAccum0\.0026\.5817\.00214\.1225\.652\.7965\.5352\.83222\.2568\.11StormAvgI0\.003\.791\.4058\.672\.650\.008\.843\.2058\.674\.12PeakI151\.6019\.4714\.22112\.0012\.546\.1036\.6827\.43112\.0028\.60PeakI301\.2013\.9210\.0080\.009\.203\.5626\.3819\.6080\.0018\.17PeakI600\.609\.037\.1150\.805\.372\.0016\.7513\.7250\.8012\.39ContribArea0\.020\.990\.417\.891\.050\.031\.250\.577\.791\.42PropHM230\.000\.460\.490\.990\.490\.010\.560\.590\.990\.30dNBR0\.010\.330\.301\.000\.260\.040\.370\.380\.870\.26KF0\.000\.210\.240\.980\.110\.000\.300\.2411\.360\.06![Refer to caption](https://arxiv.org/html/2608.05265v1/x2.png)Figure 3:Upper\-triangle Spearman rank correlation matrix for the ten benchmark features and the binary debris\-flow response\. Each cell gives a rank correlation value\. Dashed lines separate theResponserow and column from the feature–feature entries\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x3.png)Figure 4:Grand Prix–Old fire: spatial distribution of nine features across delineated basins for a representative storm\. The top row shows burn\- and soil\-related features, and the lower rows show rainfall accumulation, duration, mean intensity, and peak 15\-, 30\-, and 60\-minute intensities\.For this evaluation, we used the publicly available basin\-level dataset compiled by the USGS for post\-wildfire debris\-flow prediction\(staley\_updated\_2016\), which integrates response labels with rainfall, burn\-severity, terrain, and soil features assembled under a consistent operational workflow\. The dataset contains 1,550 basin\-storm observations from 34 fires and 244 storms collected between 2000 and 2013 across seven states in the western United States \(California, Arizona, Colorado, Idaho, Montana, New Mexico, and Utah\)\. Each row corresponds to one burned basin under one storm event, and the response is binary: debris flow is positive, and no debris flow is negative, with 334 debris\-flow and 1,216 no\-debris\-flow observations\. We chose this dataset because it underpins the operational USGS model\(staley\_prediction\_2017\)and is a standard public benchmark for this prediction task, making our model rankings comparable to prior machine learning evaluations on the same problem and USGS data lineage\(kern\_machine\_2017;nikolopoulos\_evaluation\_2018;addison\_assessment\_2019;orland\_scalable\_2022;roten\_machine\_2022\)\. The basin\-level labels and features are also relatively high quality, having been compiled under a consistent USGS operational protocol, and the sample size and class imbalance are representative of post\-wildfire debris\-flow datasets more broadly, where confirmed events are scarce, and non\-events dominate\. This makes findings on small\-sample model behavior and class\-imbalance handling more likely to transfer to other studies of this hazard\. The dataset is geographically concentrated, with 61% of records from southern California fires and 39% from the remaining states \(Figure[1](https://arxiv.org/html/2608.05265#S1.F1)\), reflecting the historical focus of the USGS operational debris\-flow warning program on the western United States, where post\-wildfire debris flows have frequently been observed\. To delineate the basin boundaries for visualization, we used the USGS postfire debris\-flow Python package\(jonathan\_m\_king\_pfdf\_2025\)\.

Table[1](https://arxiv.org/html/2608.05265#S2.T1)shows the ten numeric features used in the comparison, spanning rainfall, burn severity, basin properties, and soil erodibility\. Rainfall features areStormAccum,StormDur,StormAvgI,PeakI15,PeakI30, andPeakI60\(staley\_updated\_2016\)\. Watershed and burn\-severity features areContribArea,PropHM23, anddNBR\(key\_landscape\_2006\)\. The soil featureKFis the USLE soil erodibility factor \(KK\) for the fine \(FF\) fraction of the soil, derived from STATSGO\(schwarz\_state\_1995\)\. The pairwise Spearman rank correlations among retained features and the response \(Figure[3](https://arxiv.org/html/2608.05265#S2.F3)\) confirm that the three peak\-intensity windows are strongly intercorrelated \(ρ\>0\.90\\rho\>0\.90\) while burn\-severity and soil features are largely independent of rainfall\. This characterizes the dependence structure of the feature space, with a strongly intercorrelated rainfall\-intensity block alongside largely independent burn\-severity and soil features, and establishes the marginal response associations against which the feature importance findings are later interpreted\.

Figure[4](https://arxiv.org/html/2608.05265#S2.F4)complements this global view of feature correlations by mapping each feature across the basins of a single fire, showing the within\-fire spatial structure that the global correlations average over\. Many of these features exhibit clear spatial autocorrelation within the fire, with neighboring basins carrying similar values\. This confirms that the features vary meaningfully at the basin level rather than being uniform within a fire, supporting the use of the basin as the modeling unit for this dataset\. The within\-fire similarity among neighboring basins also reflects genuine physical structure in the data rather than label noise, and it informs the cross\-validation design and the spatial interpretation of feature importances\.

Three characteristics of the dataset make this a difficult prediction problem\. First, the two classes overlap heavily in feature space\. The per\-feature summary statistics in Table[2](https://arxiv.org/html/2608.05265#S2.T2)show that debris\-flow and no\-debris\-flow basins differ in the expected directions on rainfall intensity, burn severity, and soil erodibility, but the distributions overlap substantially, and several features are nearly indistinguishable between the two classes\. The t\-SNE projection of the feature space in Figure[8](https://arxiv.org/html/2608.05265#S3.F8)a reinforces this, with the two classes overlapping heavily in the low\-dimensional embedding and no clear cluster structure separating them\. Second, the dataset is class\-imbalanced, with only 334 of 1,550 observations labeled as debris flow, a 21\.5% positive rate, so a model can reach high accuracy by defaulting to the majority no\-debris\-flow class while failing to identify the debris\-flow events of interest\. Third, the dataset is small, with 1,550 observations across ten features, which limits the training data available to each model and increases the variance of the performance estimates\. Together, class overlap, class imbalance, and small sample size imply that no single feature or low\-dimensional combination cleanly separates the classes, motivating the comparison of model families with differing inductive biases\.

### 2\.2Cross\-Validation

Given that the dataset is small and class\-imbalanced, reliance on a single train\-test split would yield unstable estimates that are sensitive to split choice\. For that reason, we evaluated the performance of the models with repeated stratified cross\-validation\(kohavi\_study\_1995;arlot\_survey\_2010;krstajic\_cross\-validation\_2014\)\. In particular, we used stratified 5\-fold cross\-validation repeated 10 times, producing 50 outer evaluation splits with an equal percentage of no\-debris\-flow and debris\-flow observations in each split\. Splits were precomputed once with a fixed random seed and reused unchanged across all models so that every model was evaluated on identical train\-test splits\. Because 5\-fold cross\-validation was used, 20% of observations were held out for testing, and 80% were used for training\. All preprocessing and hyperparameter tuning were performed using only outer\-training data\. Within each outer training fold, a stratified 80/20 split produced an inner training set used to fit candidate model configurations and a validation set used to score them during hyperparameter tuning, keeping the held\-out test fold unseen during tuning to avoid selection bias\(cawley\_over\-fitting\_2010\)\. After hyperparameter selection, each model was refit on the full outer training fold using the chosen configuration so that the final fit used all available training data, while reported performance was measured on the held\-out test fold that was never seen during tuning or fitting\.

### 2\.3Preprocessing of Features

Preprocessing followed a consistent protocol intended to make comparisons across models methodologically fair\. Missing values were imputed within each outer training fold using a multi\-step procedure that exploited known physical relationships among the features before falling back on summary statistics\. Missing entries among the 15\-, 30\-, and 60\-minute peak rainfall intensity windows, and analogously the corresponding accumulation windows, were filled by fitting a log\-log linear regression across durations, which corresponds to the empirical power\-law intensity\-duration relationship characteristic of rainfall\(koutsoyiannis\_mathematical\_1998\)\. Decreasing peak intensity and increasing accumulation with longer durations were enforced for all rows\. MissingdNBRvalues were imputed fromPropHM23using a fold\-specific polynomial regression\. Any residual missing values were filled with the per\-feature training median\.

The imputer was fit on the outer training fold only and applied unchanged to the corresponding validation and held\-out test splits to prevent information leakage\. Features were then provided to each model in the form most appropriate to its learning assumptions\. Distance\-, kernel\-, and gradient\-based models received features standardized by a per\-fold z\-score transformation fit on the outer training data, while tree\-based models, which are invariant to monotone feature transformations, operated directly on the features\. This per\-model treatment of inputs was chosen to give each model its best chance of performing well on this dataset rather than to impose identical preprocessing across model families with different learning requirements\. Because the same imputation protocol was applied within each fold for every model, the same train\-test splits were shared across the comparison, and scaling was applied only where required by the learning algorithm, observed performance differences can be attributed to model behavior rather than to differences in data handling\.

### 2\.4Hyperparameter Tuning

Hyperparameter selection followed a single, consistent protocol designed to keep model comparison fair across model families\. For each model and training condition, we used Bayesian optimization with 100 trials and Hyperband pruning to maximize validation threat score under fold\-aware training\(akiba\_optuna\_2019;li\_hyperband\_2018\)\. The selected configuration for each model was then fixed for final evaluation so that performance differences reflect model behavior rather than repeated retuning noise\. Tuning was conducted separately for unaugmented training and for augmented training with TabPFN\-generated synthetic observations \(Section[2\.13](https://arxiv.org/html/2608.05265#S2.SS13)\), so each configuration matched the training distribution under that condition\. Not all models were included in this shared tuning framework\. Some have no applicable hyperparameters under our setup, and others use their own native tuning procedures\. These differences are described alongside each model in the sections that follow\. As with the preprocessing protocol, these per\-model differences reflect a deliberate choice to give each model its best chance of performing well rather than to enforce an artificially uniform tuning procedure across model families that would favor some models over others\.

Two model\-fitting choices that affect performance under class imbalance were held fixed rather than tuned per fold\. Class\-weight rebalancing was applied through each library’s standard balanced\-weighting mode where supported, since this is the appropriate default for imbalanced binary classification\(he\_learning\_2009\)and is not data\-dependent in a way that would benefit from per\-fold tuning\. Decision\-threshold optimization was likewise not performed to avoid conflating hyperparameter selection with the metric used for model comparison\.

### 2\.5Models Compared

The 15 models span five model families, chosen so that the benchmark compares substantively different ways of learning from the same dataset rather than minor variants of one model architecture\.

[*Logistic Regression Models*](https://arxiv.org/html/2608.05265#S2.SS6)are represented by two logistic regression variants: the operational Staley17 model\(staley\_prediction\_2017\), restricted to four features, and a standard logistic regression\(cox\_regression\_1958\)trained on the full ten\-feature set so that the link function can be isolated from feature\-set restriction\.

[*Distance\- and kernel\-based Models*](https://arxiv.org/html/2608.05265#S2.SS7)include K\-Nearest Neighbors \(KNN\)\(cover\_nearest\_1967\)and the Support Vector Classifier \(SVC\)\(cortes\_support\-vector\_1995\)\.

[*Deep\-learning Models*](https://arxiv.org/html/2608.05265#S2.SS8)include a fully connected multi\-layer perceptron \(MLP\)\(rumelhart\_learning\_1986\)together with four models originally designed for sequential or image\-like input: a one\-dimensional convolutional network \(CNN\)\(lecun\_gradient\-based\_1998\), a long short\-term memory \(LSTM\)\(hochreiter\_long\_1997\), a Transformer encoder applied to per\-feature token projections\(vaswani\_attention\_2017\), and a Mamba state\-space model\(gu\_mamba\_2024\)\.

[*Tree\-based Models*](https://arxiv.org/html/2608.05265#S2.SS9)include Random Forest\(breiman\_random\_2001\), Extremely Randomized Trees \(ExtraTrees\)\(geurts\_extremely\_2006\), eXtreme Gradient Boosting \(XGBoost\)\(chen\_xgboost\_2016\), Categorical Boosting \(CatBoost\)\(prokhorenkova\_catboost\_2018\), and tree\-structured Recursive Feature Machines \(xRFM\)\(beaglehole\_xrfm\_2025\)\.

[*Tabular Prior\-Data Fitted Network \(TabPFN\)*](https://arxiv.org/html/2608.05265#S2.SS10)is in a family of its own; it represents a model class that performs prediction via in\-context learning, pre\-trained on millions of synthetic tasks\(hollmann\_tabpfn\_2023\)\.

### 2\.6Logistic Regression Models

TheStaley17\(staley\_prediction\_2017\)logistic regression model employed in the USGS operational system defines the link function as:

χ=β\+C1​T​R\+C2​F​R\+C3​S​R\\chi=\\beta\+C\_\{1\}TR\+C\_\{2\}FR\+C\_\{3\}SR\(1\)whereC1C\_\{1\},C2C\_\{2\}, andC3C\_\{3\}are learned weights,RRdenotes peak 15\-minute rainfall accumulation in mm,TTisPropHM23,FFisdNBR, andSSisKF\. This formulation ensures that whenR=0R=0mm, the predicted probability asymptotically approaches zero for sufficiently negativeβ\\betavalues\. Among the many candidate features examined by Staley et al\., these four yielded the strongest predictive performance and were retained in the published model formulation\(staley\_prediction\_2017\)\. Because the other 14 models consumePeakI15directly as a 15\-minute peak intensity, we derived the accumulationRRrequired by the Staley17 formulation asR=PeakI15×0\.25R=\\textit\{PeakI15\}\\times 0\.25h, which preserves the published equation while keeping all models conditioned on the same rainfall information\. In addition to Staley17, we included a standardLogistic Regression\(cox\_regression\_1958\)trained on the full set of available features to isolate the effect of feature restriction from the choice of model family\. Staley17 was excluded from the shared hyperparameter tuning because it faithfully reproduces the original formulation\. Logistic Regression was included\.

### 2\.7Distance\- and Kernel\-based Models

Distance\- and kernel\-based models form class boundaries by finding local neighborhoods or maximum\-margin separators in transformed feature space\(hastie\_elements\_2009\)\.K\-Nearest Neighbors \(KNN\)\(cover\_nearest\_1967\)predicts class membership from the labels of nearby training observations under a chosen distance metric, with decision smoothness controlled by the neighbor count and distance weighting\.Support Vector Classifier \(SVC\)\(cortes\_support\-vector\_1995\)learns a separating hyperplane with maximum margin and, through nonlinear kernels, can represent curved class boundaries in the ten\-dimensional feature space\. Both of these models used the shared hyperparameter tuning framework, with the tuned hyperparameters including neighborhood size and distance weighting for KNN and kernel, regularization, and kernel\-scale settings for SVC\.

### 2\.8Deep\-learning Models

As a feed\-forward deep\-learning baseline, we included anMLP\(rumelhart\_learning\_1986\)with fully connected layers and nonlinear activations, which imposes no structural prior over the input features\. We then evaluated four models whose inductive biases instead target sequential, spatial, or token\-structured features: a one\-dimensionalCNN\(lecun\_gradient\-based\_1998\), anLSTM\(hochreiter\_long\_1997\), aTransformer\(vaswani\_attention\_2017\)encoder over feature tokens, and aMamba\(gu\_mamba\_2024\)state\-space model\. These are not obvious choices for a small set of tabular, non\-sequential features, but we included them to test whether convolution, recurrence, attention, or state\-space structure can nevertheless extract useful signal from the feature ordering, a question that remains open in the tabular deep\-learning literature\(grinsztajn\_why\_2022;borisov\_deep\_2024\)\. All models in this family used the shared hyperparameter tuning framework, with tuned hyperparameters covering architectural width and depth\. Optimizer and training\-schedule parameters were held fixed across all deep\-learning models\.

### 2\.9Tree\-based Models

Tree\-based models aggregate many decision trees to reduce variance and capture nonlinear interactions among features\(breiman\_random\_2001\), and tree\-based approaches have demonstrated strong performance on post\-wildfire debris\-flow prediction\(roten\_machine\_2022\)\. We included five tree\-based models\.Random Forest\(breiman\_random\_2001\)averages predictions across bootstrap\-resampled trees with random feature subsampling at each split\.Extremely Randomized Trees \(ExtraTrees\)\(geurts\_extremely\_2006\)extends this by selecting split thresholds at random, further reducing variance at a slight bias cost\.eXtreme Gradient Boosting \(XGBoost\)\(chen\_xgboost\_2016\)andCategorical Boosting \(CatBoost\)\(prokhorenkova\_catboost\_2018\)are gradient\-boosted decision tree \(GBDT\) implementations that iteratively correct residual errors with shallow trees\. XGBoost applies regularized shrinkage and depth control while CatBoost uses ordered boosting to reduce overfitting\.tree\-structured Recursive Feature Machines \(xRFM\)\(beaglehole\_xrfm\_2025\)combines feature\-learning kernel machines with an adaptive tree structure, allowing the kernel to adapt to local data geometry at each split of the input space\. All models except for xRFM used the shared hyperparameter tuning framework\. xRFM used an internal tuning framework\.

### 2\.10Tabular Prior\-Data Fitted Network \(TabPFN\)

![Refer to caption](https://arxiv.org/html/2608.05265v1/x4.png)Figure 5:Overview of TabPFN architecture, training, prediction, and synthetic data generation\. Panel \(a\) shows the TabPFN layer block, repeated twelve times, combining feature attention, sample attention, and a multi\-layer perceptron\. Panel \(b\) depicts TabPFN pre\-training, where millions of synthetic datasets are sampled from a Structural Causal Model prior and used to train the Transformer weights via a supervised loss\. Panel \(c\) applies the pre\-trained TabPFN to the USGS dataset: labeled training basins and unlabeled test basins enter the model together, and class probabilities \(DF / No DF\) emerge in one forward pass with no parameter updates\. Panel \(d\) illustrates TabPFN\-based synthetic data generation, in which class\-conditional distributions are modeled autoregressively via conditional feature modeling \(p​\(xi∣x−i,y\)p\(x\_\{i\}\\mid x\_\{\-i\},y\)\), the chain rule yields the jointp​\(X∣y\)p\(X\\mid y\), and sampling produces synthetic debris\-flow and no\-debris\-flow rows that augment the USGS training set\.Tabular Prior\-Data Fitted Network \(TabPFN\) departs from the standard fit\-from\-scratch paradigm for machine learning\. Whereas conventional models such as gradient\-boosted trees are trained anew on every dataset, TabPFN is a pre\-trained Transformer that predicts on an unseen dataset in a single forward pass by treating the labeled training examples as context\(hollmann\_tabpfn\_2023\)\. This capability is grounded in the Prior\-Data Fitted Network framework, which reformulates supervised learning as a meta\-learning task\. The model is trained offline on millions of synthetic datasets sampled from a prior combining Structural Causal Models and Bayesian Neural Networks \(Figure[5](https://arxiv.org/html/2608.05265#S2.F5)b\), learning to approximate the posterior predictive distribution without any dataset\-specific parameter updates\(muller\_transformers\_2024\)\.‘ At inference time \(Figure[5](https://arxiv.org/html/2608.05265#S2.F5)c\), each training observation is encoded as a token, and the Transformer’s self\-attention layers jointly process the training context and the unlabeled test data to produce class probabilities, analogous to few\-shot learning in large language models\(brown\_language\_2020\)but specialized to tabular inputs\. We used TabPFN\-2\.5\(grinsztajn\_tabpfn\-25\_2025\), which extends the original architecture by incorporating large\-scale open\-source tabular data alongside the synthetic prior to improve generalization across diverse data distributions\.

The internal architecture \(Figure[5](https://arxiv.org/html/2608.05265#S2.F5)a\) consists of twelve stacked layer blocks\. Each block combines feature\-wise attention, which models dependencies among features, with sample\-wise attention, which conditions predictions on the labeled training context, followed by a position\-wise multi\-layer perceptron\. This dual\-attention design allows the model to simultaneously relate features to one another and relate test points to the labeled training distribution, a form of implicit non\-parametric reasoning that requires no user\-specified model of the data\-generating process\. Three properties of this design are particularly relevant to the present comparison\. First, this deep architecture gives the model substantial capacity to learn complex, nonlinear interactions among features, well beyond the shallow functional forms of logistic regression\. Second, because the model is pre\-trained on millions of synthetic tabular datasets, it brings transferable tabular inductive biases to small tabular datasets without requiring large domain\-specific training data, a natural fit for the 1,550\-observation USGS dataset\. Third, TabPFN requires no dataset\-specific hyperparameter tuning, all model weights remain fixed after pre\-training, and the model is applied with static settings, in contrast to the Bayesian search applied to the other tuned models\.

Beyond prediction, the TabPFN software ecosystem includes an unsupervised extension that turns the same pretrained tabular Transformer into a class\-conditional synthetic data generator \(Figure[5](https://arxiv.org/html/2608.05265#S2.F5)d\)\(grinsztajn\_priorlabstabpfn\-extensions\_2026\)\. Using this model as the synthetic data generator aligns with its architectural emphasis on modeling dependencies among tabular features and with offline meta\-training across many synthetic tabular datasets, both of which support coherent conditional modeling when data is limited\.

Thus, TabPFN serves two roles in this study: as one of the fifteen evaluated models, and as the generator of the synthetic training observations used for augmentation\.

### 2\.11Evaluation Metrics

Model performance was evaluated using a suite of metrics derived from the confusion matrix\. The dataset for this comparison is class\-imbalanced, with no\-debris\-flow events outnumbering debris\-flow events, so selecting metrics that remain informative under such skew is essential\. The four outcomes of binary classification are: true positives \(T​PTP, debris\-flow events correctly predicted\), false positives \(F​PFP, debris flow predicted but not observed\), false negatives \(F​NFN, debris flow observed but not predicted\), and true negatives \(T​NTN, non\-events correctly identified\)\.

The primary metric is thethreat score\(TS\), also known as the Jaccard Score or Critical Success Index, defined as:

T​S=T​PT​P\+F​P\+F​NTS=\\frac\{TP\}\{TP\+FP\+FN\}\(2\)
In all metrics, debris flow is the positive class\. The threat score is particularly suited to hazard assessment because it disregards true negatives, which dominate imbalanced datasets, and measures how well the model recovers the event of interest relative to the total number of times that event was either predicted or observed\(schaefer\_critical\_1990\)\. Supplementary metrics reported for completeness include accuracy, precision, recall, F1\-score, and area under the receiver\-operating\-characteristic \(ROC\) curve \(AUC\)\. Accuracy shows the proportion of correctly predicted events\. Precision and recall characterize the false\-alarm and missed\-event trade\-off, respectively\. F1 is their harmonic mean\. ROC\-AUC provides a threshold\-independent measure of discrimination\(hanley\_meaning\_1982\)\. For each metric, we reported the mean and standard deviation across all 50 outer test folds to characterize both central performance and fold\-to\-fold stability\(kohavi\_study\_1995;arlot\_survey\_2010\)\. On a single inventory, fold\-to\-fold dispersion can be comparable to mean gaps among top models, so model comparisons must account for cross\-validation variability rather than table rank alone\(krstajic\_cross\-validation\_2014;nadeau\_inference\_2003\)\.

To assess whether differences between models are large relative to this variability, we compared pre\-specified pairs using paired threat scores from the same outer splits, following single\-dataset comparison guidance for models evaluated on identical resampled partitions\(dietterich\_approximate\_1998;nadeau\_inference\_2003\)\. For each pair, we reported the mean paired differenceΔ\\DeltaTS, the fraction of folds in which one model achieved a higher threat score, and repeat\-level summaries obtained by averaging threat score over the five folds within each of the ten CV repeats\. We summarized repeat\-levelΔ\\DeltaTS with its mean, standard deviation, and 95% confidence interval, and applied a two\-sided Wilcoxon signed\-rank test to the ten repeat\-level differences as a nonparametric supportive check only\(demsar\_statistical\_2006\)\. Because all folds draw from the same inventory, cross\-validation scores are not independent observations, and repeat\-level tests were not treated as full replication\(nadeau\_inference\_2003\)\. Primary emphasis was on effect size, fold win rate, and overlap of repeat\-level confidence intervals rather than on exhaustive pairwise hypothesis testing\(benavoli\_time\_2017\)\.

### 2\.12SHAP Feature Importance Evaluation

We used SHAP feature importances to interpret which features the different models rely on when predicting debris\-flow probability on the held\-out basins\(lundberg\_unified\_2017\)\. SHAP assigns each feature a value representing its marginal contribution to the model output relative to a baseline, satisfying local\-accuracy and consistency guarantees\. We chose SHAP over permutation\-based importance because its per\-observation additive scores support both global ranking and the basin\-level spatial maps used here, and because Shapley values provide a principled additive decomposition relative to an explicit baseline\.

Importance scores were computed on the same 50 outer test splits as the benchmark, so that SHAP estimates inherited the fold structure used for performance evaluation\. SHAP values are produced by an*explainer*, an algorithm that estimates each feature’s Shapley contribution for a given model class\. The choice of explainer depended on the model\. TabPFN used the library’s built\-in permutation explainer with a fixed evaluation budget\(grinsztajn\_priorlabstabpfn\-extensions\_2026\)\. Meanwhile, Random Forest, XGBoost, and CatBoost used TreeSHAP, an exact polynomial\-time algorithm for tree\-based models\(lundberg\_explainable\_2019\)\. For each fold, the model was refit on the training split using the tuned hyperparameters from the main benchmark run, and SHAP values were then computed with the corresponding explainer on a uniformly subsampled subset of the held\-out test rows\.

For basin\-level spatial analyses, we used a leave\-one\-fire\-out procedure rather than the per\-fold splits described above\. A TabPFN model was fit on all data except the target fire, and SHAP values were computed on the withheld basins using the same permutation explainer\. This ensured that all basins within a fire were explained by a single trained model, yielding spatially coherent importance polygons rather than a mosaic of per\-fold explanations\. Basin\-level maps included only basin observations from one storm event per fire, so rainfall\-related features corresponded to a single meteorological event rather than a mixture of storms on the same terrain\. We focused this analysis on the Station and Grand Prix–Old fires because they have a large amount of observations for a single storm\.

We then evaluated the SHAP feature importances along five complementary analyses\. First, we compared pooled\-fold global rankings between TabPFN and the leading tree\-based models \(Random Forest, XGBoost, CatBoost\) to identify points of agreement and divergence across model families\. Second, we examined the fold\-to\-fold stability of TabPFN’s feature importance by tracking mean\|SHAP\|\|\\mathrm\{SHAP\}\|per feature across the 50 stratified outer folds\. Third, we mapped the resulting basin\-level TabPFN SHAP values onto delineated watersheds, which addresses the limitation that global rankings cannot distinguish whether TabPFN conditions on basin\-level feature variation within a fire\. Fourth, we related basin\-level feature values to their SHAP contributions through Spearman rank correlation within each fire, which characterizes the dependence between feature value and importance independent of map geometry\. Fifth, we tested whether the sign of the SHAP–feature dependence aligns with established physical expectations by computing the Spearman correlation between each feature’s values and its corresponding SHAP contributions pooled across all 50 folds for the four leading models, treating\|r\|<0\.15\|r\|<0\.15as a weak or negligible association\. We defined those expectations from previous post\-wildfire debris\-flow studies\(cannon\_predicting\_2010;staley\_updated\_2016;staley\_prediction\_2017\)as a positive association between each feature and predicted debris\-flow probability for all features exceptStormDurandContribArea, for which no fixed directional match between classes was assumed\.

Throughout this study, SHAP values are read as model\-level feature reliance applicable to this study area and comparable burned\-watershed settings, not as universal physical drivers of debris flows in all burned landscapes\.

### 2\.13Synthetic Data Generation

Synthetic observations were generated using the TabPFN unsupervised extension\(grinsztajn\_priorlabstabpfn\-extensions\_2026\)\. Within each outer fold, a class\-conditional generator was fit per class on that class’s real training rows, and observations were drawn autoregressively one feature at a time from the learned conditional distributions, producing coherent multivariate observations without a fixed parametric form\.

Generator fidelity was assessed per fold using the two\-sample Kolmogorov–Smirnov statistic and Wasserstein\-1 distance per feature, together with the multivariate energy distance between real and synthetic observations within each class\(massey\_kolmogorov\-smirnov\_1951;villani\_wasserstein\_2009;szekely\_energy\_2013\)\. These three metrics capture, respectively, the largest pointwise gap between univariate empirical distributions, the overall transport distance between those distributions, and any multivariate dependence mismatch across features\.

### 2\.14Synthetic Data Augmentation Evaluation

Augmentation was treated as a fixed study component rather than a benchmark dimension, so a single synthetic generator was used throughout\. TabPFN’s autoregressive conditional modeling was chosen because it is designed for small tabular datasets and is available within the same modeling framework\. The objective was to measure how this one consistent mechanism affects downstream model performance, not to compare generators\. Alternative tabular synthesis approaches, such as SMOTE\-style interpolation\(chawla\_smote\_2002\)and conditional generative adversarial networks\(xu\_modeling\_2019\), are therefore outside the scope of this study\.

For each fold, synthetic observations were generated such that the count per class matched the real count of the opposing class, aiming for a balanced distribution while preserving the empirical per\-class feature distributions\. However, to prevent synthetic data from overwhelming the training signal, the real observations were duplicated once before appending the synthetic rows, resulting in a final training set that was partially rebalanced rather than perfectly uniform\. Crucially, test sets remained entirely real, and all synthetic observations for a given fold were generated once and shared identically across all models, ensuring that benchmarking differences reflected model behavior rather than generator variance\.

All 15 models were then compared again under the augmented training condition using the same outer cross\-validation, inner tuning, and evaluation metrics as the unaugmented comparison, isolating the effect of augmentation on each model\.

## 3Results

We first compared models under unaugmented training\. Then we interpreted model behavior with SHAP\-based feature importances\. Finally, we evaluated the quality of the synthetic data generation and the impact of augmentation on model performance\.

### 3\.1Model Comparison

![Refer to caption](https://arxiv.org/html/2608.05265v1/x5.png)Figure 6:Mean ROC curves under unaugmented training, with one curve per model averaged over the same 50 stratified outer folds \(10 repeats of 5 folds\) and a dashed diagonal at random chance\. The legend lists models in descending unaugmented threat score order, and the bottom\-right inset magnifies the upper\-left region\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x6.png)Figure 7:Per\-model ROC traces under unaugmented training, arranged in a three\-row by five\-column grid sorted by descending mean AUC, with Staley17 in the bottom\-right panel\. Within each panel, thin green lines are the 50 individual outer\-fold ROC curves, the navy line is the fold mean, the gray dashed lines are the±1\\pm 1standard deviation envelope, and the red dashed diagonal marks random chance\.Table 3:Fifteen models under unaugmented training on the same 50 stratified outer folds\. For each metric,μ\\muis the mean across folds andσ\\sigmais the standard deviation\. Rows are sorted by descending mean threat score \(TS\)\. Within each metric, the three highest means are tinted by rank \(1st: sky blue; 2nd: orange/gold; 3rd: bluish green\)\.ModelAccF1RecallPrecTSμ\\muσ\\sigmaμ\\muσ\\sigmaμ\\muσ\\sigmaμ\\muσ\\sigmaμ\\muσ\\sigmaStaley170\.7830\.0130\.2620\.0490\.1800\.0420\.5040\.0920\.1520\.032LogisticRegression0\.7730\.0260\.5820\.0320\.7320\.0460\.4850\.0380\.4120\.032Mamba0\.7660\.0300\.5970\.0330\.8010\.0610\.4780\.0420\.4260\.034Transformer0\.7690\.0330\.6040\.0320\.8110\.0600\.4840\.0440\.4330\.034LSTM0\.7770\.0420\.6060\.0510\.7910\.0670\.4950\.0600\.4370\.053CNN0\.7890\.0240\.6270\.0290\.8220\.0470\.5080\.0340\.4570\.031MLP0\.7940\.0320\.6280\.0450\.8020\.0630\.5180\.0500\.4590\.048SVC0\.8340\.0180\.6660\.0320\.7690\.0510\.5890\.0360\.5000\.036xRFM0\.8660\.0150\.6730\.0420\.6410\.0670\.7140\.0500\.5090\.048KNN0\.8990\.0150\.7550\.0390\.7250\.0600\.7920\.0450\.6080\.051RandomForest0\.8980\.0180\.7640\.0400\.7690\.0570\.7630\.0540\.6200\.052XGBoost0\.8980\.0180\.7650\.0410\.7700\.0590\.7630\.0510\.6210\.054CatBoost0\.8940\.0160\.7660\.0340\.8090\.0570\.7300\.0420\.6220\.045ExtraTrees0\.8940\.0190\.7770\.0350\.8500\.0460\.7170\.0490\.6360\.047TabPFN0\.9030\.0160\.7770\.0360\.7840\.0570\.7750\.0500\.6370\.048![Refer to caption](https://arxiv.org/html/2608.05265v1/x7.png)Figure 8:Two\-panel t\-SNE\(maaten\_visualizing\_2008\)projection \(perplexity 30\) of the full dataset, with class\-wise marginal histograms along each axis\. Panel \(a\) uses preprocessed raw features and panel \(b\) uses out\-of\-fold TabPFN internal representations\.Table[3](https://arxiv.org/html/2608.05265#S3.T3)summarizes the comparison under unaugmented training\. TabPFN ranked first with a mean threat score of 0\.637, followed by a narrow top tier \(span 0\.017\) comprising ExtraTrees \(0\.636\), CatBoost \(0\.622\), XGBoost \(0\.621\), and Random Forest \(0\.620\), with fold dispersion \(σ\\sigma\) between 0\.045 and 0\.054 across these five models\. TabPFN and ExtraTrees tied for the highest F1 \(both 0\.777\), and ExtraTrees achieved the highest recall \(0\.850\)\. Relative to this fold\-to\-fold variability, the leading models are closely matched\. The mean gap between TabPFN and ExtraTrees is onlyΔ\\DeltaTS=\+0\.001=\+0\.001, TabPFN exceeded ExtraTrees on 30 of 50 paired folds, and the repeat\-level 95% confidence interval forΔ\\DeltaTS \(−0\.007\-0\.007to\+0\.009\+0\.009\) includes zero\. Among the top five models, repeat\-level rank 1 was split across ExtraTrees \(five repeats\), TabPFN \(four\), and CatBoost \(one\), so table order should not be read as a stable ordinal ranking within the top tier\. By contrast, TabPFN exceeded the operational Staley17 baseline byΔ\\DeltaTS=\+0\.486=\+0\.486on every fold with repeat\-level 95% CI\[\+0\.478,\+0\.493\]\[\+0\.478,\+0\.493\], and exceeded CatBoost byΔ\\DeltaTS=\+0\.016=\+0\.016\(36 of 50 folds; repeat\-level CI\[\+0\.005,\+0\.026\]\[\+0\.005,\+0\.026\]\)\. Below the top tier, SVC \(0\.500\) and xRFM \(0\.509\) occupied a middle band, the four sequential or image\-oriented deep\-learning models clustered at 0\.426–0\.459 \(Mamba, Transformer, LSTM, CNN\) with high recall \(0\.79–0\.82\) but low precision \(0\.48–0\.51\), the ten\-feature Logistic Regression reached 0\.412, and the four\-feature Staley17 model trailed at 0\.152 with a recall of 0\.180\.

Figure[6](https://arxiv.org/html/2608.05265#S3.F6)shows the mean ROC curves averaged across the 50 outer folds for each model\. The threshold\-independent ordering reproduces the threat score ranking: TabPFN \(AUC 0\.953\), ExtraTrees \(0\.948\), Random Forest \(0\.945\), XGBoost \(0\.944\), and CatBoost \(0\.939\) trace the upper\-left envelope, while the Staley17 model is relegated to a markedly lower curve at AUC 0\.805\.

Figure[7](https://arxiv.org/html/2608.05265#S3.F7)resolves the mean curves into their 50 per\-fold traces and±1\\pm 1standard deviation envelopes\. TabPFN and the leading tree\-based models occupy the most favorable ROC region across folds with tightly bundled traces, whereas the deep\-learning and distance\-based models exhibit visibly wider standard deviation envelopes\.

### 3\.2Feature Importance Evaluation

Table 4:SHAP\-based feature rankings by model \(1 = most important\)\. Rankings from mean\|SHAP\|\|\\mathrm\{SHAP\}\|per feature\. Feature numbers show overall rank averaged across models\.FeatureRandomForestXGBoostCatBoostTabPFN1\. StormAccum22122\. PeakI1541213\. PeakI3016834\. PropHM2355385\. PeakI60310556\. KF103477\. StormAvgI74968\. StormDur691049\. dNBR887910\. ContribArea97610![Refer to caption](https://arxiv.org/html/2608.05265v1/x8.png)Figure 9:Combined SHAP summary for TabPFN, pooling held\-out basin explanations across all 50 stratified outer folds\. The left panel shows mean\|SHAP\|\|\\mathrm\{SHAP\}\|per feature as a horizontal bar chart, and the right panel shows the corresponding beeswarm in which each dot is one basin\-level explanation plotted by its SHAP contribution to log\-odds and colored by feature value from low to high\. Features are ordered by mean\|SHAP\|\|\\mathrm\{SHAP\}\|so both panels share the same feature axis\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x9.png)Figure 10:Fold\-wise stability of TabPFN SHAP feature importance as a heatmap, with rows for the ten features and columns for the 50 stratified outer folds\. Each cell gives the mean\|SHAP\|\|\\mathrm\{SHAP\}\|on the fold’s held\-out basins\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x10.png)Figure 11:Agreement between SHAP\-derived feature directions and physical expectations\. Orange = positive direction \(toward debris flow\); teal = negative; gray = weak \(\|r\|<0\.15\|r\|<0\.15\)\. Y/N indicates agreement or contradiction with physics;\+\+/−\-is shown directly for features with ambiguous physical sign\. All four models correctly capture the expected positive direction for burn severity, terrain steepness, soil erodibility, and rainfall intensity\.The SHAP feature importance evaluation comprises four analyses: global rankings, fold\-to\-fold stability, basin\-level spatial structure within fires, and alignment with physical expectations\. Rainfall intensity dominated global importance, while burn severity and soil erodibility showed the most consistent within\-fire structure, and the four leading models agreed on physically expected directions\.

The beeswarm summary \(Figure[9](https://arxiv.org/html/2608.05265#S3.F9)\) ranks peak 15\-minute intensity, storm accumulation, and peak 30\-minute intensity as the top three features by mean\|SHAP\|\|\\mathrm\{SHAP\}\|, together carrying the bulk of the feature importance\. Table[4](https://arxiv.org/html/2608.05265#S3.T4)compares SHAP\-derived rankings for TabPFN with high\-performing tree\-based models\. Rainfall metrics ranked highly across all four models, but burn severity and soil erodibility varied substantially\. Soil erodibility \(KF\) ranked 3–4 for XGBoost and CatBoost but 10 for Random Forest, andPropHM23ranked 3–5 for the tree\-based models but 8 for TabPFN\.

The fold\-stability heatmap \(Figure[10](https://arxiv.org/html/2608.05265#S3.F10)\) shows that the three rainfall features \(PeakI15,StormAccum, andPeakI30\) maintain high mean\|SHAP\|\|\\mathrm\{SHAP\}\|across nearly all folds, whereasdNBRandContribAreaexhibit lower and more variable importance from fold to fold\.

Figures[13](https://arxiv.org/html/2608.05265#S4.F13)and[14](https://arxiv.org/html/2608.05265#S4.F14)show SHAP maps and dependence scatters for the Station and Grand Prix–Old fires across four features \(PeakI15,PropHM23,dNBR, andKF\)\. The spatial maps \(top and middle rows\) reveal contrasting behavior between the two fires\. On Grand Prix–Old, all four SHAP fields visually co\-locate with their corresponding feature fields, with high\-intensity, high\-severity, and high\-erodibility basins carrying positive SHAP\. On Station, only thedNBRandKFSHAP fields show the same visual co\-location with their feature fields, while thePeakI15andPropHM23SHAP fields appear more uniform across basins than their feature fields\. The basin\-level dependence scatters \(bottom rows\) quantify these patterns independent of map geometry\. On Grand Prix–Old, all four features show strong positive Spearman associations between feature value and SHAP contribution \(ρ=0\.67\\rho=0\.67–0\.920\.92\)\. On Station,dNBRandKFshow similarly strong positive dependence \(ρ=0\.71\\rho=0\.71and0\.750\.75\), whereasPeakI15andPropHM23show only weak associations \(ρ=0\.27\\rho=0\.27–0\.450\.45\)\. Burn severity and soil erodibility therefore exhibit more consistent spatial structure and SHAP–feature coupling across both fires than peak 15\-minute intensity and terrain steepness at these events\.

Figure[11](https://arxiv.org/html/2608.05265#S3.F11)summarizes the SHAP–feature sign check for the four leading models against the consensus physical expectation\. All four correctly recovered a positive direction for burn severity, terrain steepness, soil erodibility, and every rainfall metric, with no contradictions observed for features with a clear physical sign\.

### 3\.3Synthetic Data Augmentation Evaluation

Table 5:Fifteen models under augmented training on the same 50 stratified outer folds and split assignments as Table[3](https://arxiv.org/html/2608.05265#S3.T3)\. For each metric,μ\\muis the mean across folds with synthetic observations added during training, andΔ\\Deltais the change relative to the unaugmented run\. Rows are sorted by descending mean threat score \(TS\)\. Within each metric, the three largest improvements are tinted by rank \(1st: sky blue; 2nd: orange/gold; 3rd: bluish green\)\.ModelAccF1RecallPrecTSμ\\muΔ\\Deltaμ\\muΔ\\Deltaμ\\muΔ\\Deltaμ\\muΔ\\Deltaμ\\muΔ\\DeltaStaley170\.772\-0\.0110\.443\+0\.1810\.421\+0\.2410\.472\-0\.0310\.285\+0\.133LogisticRegression0\.773\+0\.0000\.587\+0\.0040\.745\+0\.0130\.486\+0\.0000\.416\+0\.004CNN0\.787\-0\.0020\.624\-0\.0030\.819\-0\.0030\.506\-0\.0030\.454\-0\.003Mamba0\.790\+0\.0250\.626\+0\.0300\.812\+0\.0110\.512\+0\.0330\.457\+0\.031Transformer0\.786\+0\.0170\.627\+0\.0230\.830\+0\.0190\.506\+0\.0210\.457\+0\.024LSTM0\.829\+0\.0520\.677\+0\.0700\.826\+0\.0340\.576\+0\.0810\.513\+0\.076SVC0\.839\+0\.0050\.688\+0\.0220\.825\+0\.0570\.592\+0\.0030\.526\+0\.025MLP0\.841\+0\.0470\.698\+0\.0700\.853\+0\.0510\.593\+0\.0750\.538\+0\.078xRFM0\.864\-0\.0030\.722\+0\.0490\.821\+0\.1800\.647\-0\.0670\.566\+0\.058KNN0\.901\+0\.0020\.774\+0\.0190\.789\+0\.0640\.763\-0\.0290\.632\+0\.024ExtraTrees0\.899\+0\.0050\.778\+0\.0020\.819\-0\.0310\.744\+0\.0260\.638\+0\.002XGBoost0\.899\+0\.0010\.780\+0\.0150\.826\+0\.0560\.740\-0\.0230\.640\+0\.019RandomForest0\.900\+0\.0020\.780\+0\.0160\.826\+0\.0570\.742\-0\.0220\.641\+0\.020TabPFN0\.897\-0\.0060\.781\+0\.0030\.846\+0\.0620\.727\-0\.0490\.641\+0\.004CatBoost0\.899\+0\.0060\.782\+0\.0160\.840\+0\.0310\.734\+0\.0040\.643\+0\.022![Refer to caption](https://arxiv.org/html/2608.05265v1/x11.png)Figure 12:Two\-panel t\-SNE\(maaten\_visualizing\_2008\)comparison of TabPFN\-generated synthetic training rows against real rows for a single representative outer fold \(fold 0\), with class\-wise marginal histograms along each axis\. All four groups are projected under one global t\-SNE fit \(perplexity 30\), so positions are comparable across panels\. Panel \(a\) shows real debris\-flow basins and synthetic debris\-flow observations, and panel \(b\) shows real no\-debris\-flow basins and synthetic no\-debris\-flow observations\.Before evaluating the impact of augmentation using synthetic observations, we assessed the realism of TabPFN\-generated synthetic observations against the real training data\. Figure[15](https://arxiv.org/html/2608.05265#S4.F15)shows class\-conditional kernel density estimates for real and synthetic observations across the benchmark features, with per\-feature Kolmogorov–Smirnov and Wasserstein summaries\. The synthetic\-versus\-real comparison indicates close agreement between classes\. The two\-sample Kolmogorov–SmirnovDDranges from 0\.024 to 0\.154 across features, with soil erodibility \(KF\) the largest single\-feature discrepancy and rainfall features mostly belowD=0\.10D=0\.10, while the multivariate energy distance is 0\.020 for no\-debris\-flow observations and 0\.013 for debris\-flow observations\.

Figure[12](https://arxiv.org/html/2608.05265#S3.F12)shows a two\-panel t\-SNE for one representative outer fold with class\-wise marginal histograms, where panel \(a\) compares real versus synthetic debris\-flow observations and panel \(b\) compares real versus synthetic no\-debris\-flow observations under a single t\-SNE fit computed jointly across all four groups\. The t\-SNE embedding confirms the agreement between synthetic and real observations visually, the synthetic points overlay the real point cloud in both classes without forming isolated sub\-clusters\.

Introducing the synthetic observations into training reshuffled the top tier within a narrow band \(Table[5](https://arxiv.org/html/2608.05265#S3.T5)\)\. The five highest\-threat\-score models, CatBoost \(0\.643\), TabPFN \(0\.641\), Random Forest \(0\.641\), XGBoost \(0\.640\), and ExtraTrees \(0\.638\), tightened into a span of only 0\.005 in mean threat score, smaller than their typical foldσ\\sigma\(0\.040–0\.042\)\. TabPFN exceeded ExtraTrees on 27 of 50 paired folds with a meanΔ\\DeltaTS=\+0\.003=\+0\.003and a repeat\-level 95% confidence interval that includes zero \(−0\.006\-0\.006to\+0\.012\+0\.012\), and CatBoost led TabPFN byΔ\\DeltaTS=−0\.003=\-0\.003\(25 of 50 folds in favor of TabPFN\)\. Repeat\-level rank 1 among these five models was distributed across CatBoost \(three repeats\), ExtraTrees \(three\), Random Forest \(two\), and XGBoost \(two\), with no single model dominating\. CatBoost moved from third under unaugmented training to first under augmentation and also achieved the highest augmented F1 \(0\.782\) among these top\-performing models\.

The largest gains appeared for the most restricted models\. Staley17 threat score rose by 13\.3 percentage points and its recall by 24\.1 percentage points, while the ten\-feature Logistic Regression gained only 0\.4 percentage points\. Among the deep\-learning models, LSTM gained 7\.6 percentage points in threat score and the MLP 7\.8 percentage points, while CNN declined slightly by 0\.3 percentage points\. Among middle\-tier models, xRFM improved by 5\.8 percentage points, SVC by 2\.5 percentage points, and KNN by 2\.4 percentage points\. The leading tree\-based models and TabPFN improved by at most 2\.2 percentage points in threat score\. Figure[16](https://arxiv.org/html/2608.05265#S4.F16)summarizes the per\-model change in threat score\.

## 4Discussion

### 4\.1Model Comparison

TabPFN ranked first by threat score \(0\.637\) and accuracy \(0\.903\) but tied ExtraTrees on F1 \(both 0\.777\), and ExtraTrees achieved the highest recall \(0\.850\)\. The mean threat\-score margin between TabPFN and ExtraTrees \(Δ\\DeltaTS=\+0\.001=\+0\.001\) is far smaller than fold\-to\-fold dispersion \(σ≈0\.05\\sigma\\approx 0\.05\), and repeat\-level comparisons do not support a durable ordering within the top tier, even though TabPFN clearly exceeds the Staley17 baseline and mid\-tier models by large, stable margins\. The leading five nonlinear models, therefore, form a practical top tier rather than a single universally dominant model\. The mean ROC curves in Figure[6](https://arxiv.org/html/2608.05265#S3.F6)reproduce this ordering under a threshold\-independent criterion, with TabPFN and the four leading tree\-based models tracing the upper\-left envelope and the Staley17 model relegated to a markedly lower curve, so the threat score ranking is not an artifact of the operating threshold\. This top\-tier interpretation also explains why repeated cross\-validation is necessary for model selection in this dataset\(kohavi\_study\_1995;krstajic\_cross\-validation\_2014\)\. Because the fold\-level variation of each leading model exceeds the mean separation among those models, repeated evaluation supports a stable top\-tier conclusion while discouraging over\-interpretation of the exact within\-tier order from any single split\. With top\-tier performance this tightly grouped, the cost of obtaining it becomes the practical differentiator\. Tree\-based models such as ExtraTrees, CatBoost, and XGBoost required Bayesian hyperparameter search, whereas TabPFN was run with static settings with no dataset\-specific tuning and still matched or exceeded them\. For users with limited computational resources or rapid\-response requirements, this tuning\-free property could be as operationally relevant as the marginal threat score difference\.

The tightness of the top tier is more consistent with a dataset\-imposed ceiling than with the models being functionally equivalent\. The 1,550 basin\-storm records and the southern California concentration of training observations jointly bound the achievable threat score regardless of model family\. The two\-panel t\-SNE in Figure[8](https://arxiv.org/html/2608.05265#S3.F8)corroborates this directly: panel \(a\) shows the two classes overlapping heavily in the raw feature space, while panel \(b\) shows substantial separation under TabPFN’s out\-of\-fold internal representations\. Even in panel \(b\), residual overlap persists at the class boundary, indicating that TabPFN’s learned embedding reduces but does not eliminate the class overlap in the available features\. Under these constraints, the tight clustering of top\-tier threat scores reflects the limits of the available data as much as the capabilities of the models themselves\.

The substantial gap between the logistic models and this top tier indicates that the feature\-to\-response mapping is not adequately represented by a fixed linear decision surface\. A linear model cannot express interaction terms or threshold effects among features, so both the four\-feature Staley17 model \(threat score 0\.152\) and the ten\-feature Logistic Regression \(0\.412\) trail the top nonlinear tier by 0\.2 threat score or more, consistent with prior work showing that tree\-based and neural models improve prediction over logistic regression\(addison\_assessment\_2019;roten\_machine\_2022\)\. The intermediate value of the ten\-feature Logistic Regression decomposes the Staley17\-to\-top\-tier gap into two roughly equal contributions\. First, moving from Staley17 to the ten\-feature Logistic Regression recovers about 0\.26 threat score by relaxing the four\-feature restriction alone, and second, moving from the ten\-feature Logistic Regression to the top tier recovers a further 0\.22 by relaxing the linear decision surface\. Feature breadth and decision\-surface flexibility, therefore, both contribute to improved performance\.

The four sequential or image\-oriented deep\-learning models cluster at threat score 0\.43–0\.46 with uniformly high recall \(0\.79–0\.82\) and low precision \(0\.48–0\.51\), indicating systematic over\-prediction of debris\-flow occurrence rather than training instability\. The per\-fold ROC traces in Figure[7](https://arxiv.org/html/2608.05265#S3.F7)corroborate this reading\. The deep\-learning panels show visibly wider±1\\pm 1SD envelopes than the top tree\-based models and TabPFN, so their lower mean skill is accompanied by larger fold\-to\-fold variability rather than a tightly clustered but biased decision surface\. A plausible mechanism is a mismatch between model capacity and effective sample size\. Convolutional, recurrent, attentional, and state\-space architectures are designed for image, sequence, or language domains and impose no inductive prior suited to tabular data, leaving them prone to overfitting on a dataset of this size\. This pattern is directionally consistent with tabular benchmarks, in which tree\-based models often outperform high\-capacity deep learning when no structure matched to tabular data is imposed\(grinsztajn\_why\_2022\)\. The comparison is only indicative here because the inventory has roughly 1,550 observations and because the deep models apply convolutional, recurrent, attentional, and state\-space architectures to a fixed feature ordering rather than standard tabular pipelines\.

Several constraints should remain explicit, and they are directly connected to how well these operational considerations transfer beyond the fire settings represented in the inventory\. The dataset is geographically concentrated, with 61% of records from southern California fires, which may favor features and model structures that generalize well within that region but perform differently in the drier, more continental environments represented by the remaining fires\. Because all observations are from the first two post\-wildfire years, the models are not evaluated under later recovery conditions, when vegetation re\-establishment and soil stabilization likely reduce debris\-flow susceptibility\(degraff\_timing\_2015\)\. Moderate class imbalance, with a positive\-event rate of 21\.5%, is a further limitation\. Threat score remains informative at this rate, but operational deployment would still require recalibrating decision thresholds to the false\-alarm and missed\-event tolerances of the intended warning context\. Future work extending evaluation to multi\-region datasets beyond the western United States and incorporating longer post\-wildfire temporal windows would directly test whether these geographic and temporal limits affect transferability and robustness\.

### 4\.2Feature Importance Evaluation

![Refer to caption](https://arxiv.org/html/2608.05265v1/x12.png)Figure 13:Basin\-level SHAP feature importance for the Station fire in a three\-row, four\-column layout forPeakI15,PropHM23,dNBR, andKF\. The top row shows feature\-value maps and the middle row shows SHAP\-value maps under a shared diverging colormap\. The bottom row shows dependence scatter plots of SHAP value versus feature value, with teal dots for basins, a solid orange LOWESS curve, a dashed dark\-gray OLS fit, a dotted blue rank\-based fit, a horizontal gray reference at zero, and a text box reporting Spearmanρ\\rho\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x13.png)Figure 14:Basin\-level SHAP feature importance for the Grand Prix–Old fire under the same three\-row, four\-column layout as Figure[13](https://arxiv.org/html/2608.05265#S4.F13)forPeakI15,PropHM23,dNBR, andKF\. Row encoding follows Figure[13](https://arxiv.org/html/2608.05265#S4.F13): feature\-value maps on per\-column sequential colormaps, SHAP\-value maps on a shared diverging colormap, and dependence scatter plots with a solid orange LOWESS curve, a dashed dark\-gray OLS line, a dotted blue rank\-based line, a horizontal gray reference at zero, and a Spearmanρ\\rhotext box\.SHAP analysis clarifies which features the evaluated models use when predicting debris\-flow probability and supports interpretability for practitioners\(rundel\_interpretable\_2024;lundberg\_unified\_2017\)\. For the four models we explained with SHAP, rainfall features occupy the top tier in every ranking \(Table[4](https://arxiv.org/html/2608.05265#S3.T4)\), indicating cross\-model agreement on which features these models emphasize when predicting debris\-flow probability\. Short\-window intensity and storm accumulation receive the largest contributions, in line with the intensity–duration framing of the operational Staley17 model developed for the same regional warning context\(staley\_prediction\_2017;staley\_updated\_2016\)\. Within the rainfall tier itself, the three peak\-intensity windows are strongly intercorrelated \(ρ\>0\.90\\rho\>0\.90in Figure[3](https://arxiv.org/html/2608.05265#S2.F3)\), so the specific rainfall feature each model promotes can shift while agreement among the four models on rainfall as the dominant predictive signal is preserved\. Below the rainfall tier, the rankings diverge substantially\. Burn severity \(PropHM23,dNBR\), soil erodibility \(KF\), and contributing area receive different positions across the four models, indicating that agreement on the dominant predictive signal does not extend to which secondary feature carries residual predictive content\. Contributing area \(ContribArea\) ranks consistently low across all four models, suggesting that basin size contributes little marginal signal once rainfall and burn\-severity features are present\.

The most striking divergences below the rainfall tier admit plausible mechanistic readings\.KFranks 3–4 for XGBoost and CatBoost but 10 for Random Forest\. Gradient\-boosted trees can exploit a well\-placedKFthreshold as a sequential corrective split once a rainfall feature has captured the dominant predictive signal, whereas Random Forest’s split\-level random feature subsampling means rainfall features are usually available at each split andKFrarely surfaces as the chosen feature\(louppe\_understanding\_2013\)\.PropHM23ranks 3–5 for the tree\-based models but 8 for TabPFN\. A plausible explanation is that TabPFN’s in\-context attention routes burn\-severity signal through interaction terms with rainfall context that do not register cleanly in first\-order marginal importance\. These divergent secondary\-feature rankings across the four explained models reinforce that no single model’s SHAP ranking should be read as a definitive statement of physical importance\.

The four features Staley17 retained \(PropHM23,dNBR,KF, and 15\-min rainfall accumulation\) map only partially onto the SHAP rankings of the modern models\. Rainfall andKFrank highly, butdNBRranks last or near\-last across all four models, andPropHM23is the lowest\-ranked Staley17 feature in TabPFN \(rank 8 of 10\)\. This is not a direct contradiction of the Staley17 model\. It was selected to maximize logistic\-regression skill under a four\-feature constraint, and the features it retains may exploit different structural niches than those on which modern nonlinear models most rely\. The constructive reading is that rainfall,KF, andPropHM23receive consistently high SHAP ranks across the four explained models, whiledNBRcontributes little marginal SHAP signal once those features are present\.

TabPFN’s mean absolute SHAP values \(Figure[9](https://arxiv.org/html/2608.05265#S3.F9), left panel\) place short\-window peak intensity \(15 min\), storm accumulation, and 30\-minute peak intensity as the three largest contributors to predicted debris\-flow probability, followed by storm duration and 60\-minute peak intensity\. The beeswarm panel of Figure[9](https://arxiv.org/html/2608.05265#S3.F9)resolves the sign and per\-basin spread behind these magnitudes\. For the top rainfall features, high values sit on the positive side and low values on the negative side, confirming a monotonic mapping from rainfall intensity and accumulation to predicted debris\-flow probability at the basin level rather than an averaging artifact\. Lower\-ranked features cluster more tightly near zero with mixed color, indicating that their small mean\|SHAP\|\|\\mathrm\{SHAP\}\|reflects genuinely weak per\-observation effects rather than cancellation between large positive and negative contributions\.

The fold\-stability heatmap \(Figure[10](https://arxiv.org/html/2608.05265#S3.F10)\) shows that rainfall features maintain high mean\|SHAP\|\|\\mathrm\{SHAP\}\|across nearly all 50 outer folds, whereas burn ratio \(dNBR\) and contributing area exhibit lower and more variable importance\. This asymmetry indicates that the rainfall\-driven importances are reproducible across training folds, while the burn\-severity and terrain importances are more sensitive to which folds the model is fit on\. The rainfall ranks can therefore be read as a stable property of TabPFN feature importance, while thedNBRand contributing\-area ranks are less trustworthy as point estimates and are better reported with fold\-level uncertainty\. This within\-model instability also lines up with the cross\-model picture in Table[4](https://arxiv.org/html/2608.05265#S3.T4)\. The features whose TabPFN importance is least fold\-stable are the same features whose rank disagrees most across CatBoost, Random Forest, and XGBoost, so low\-signal features are unstable both within TabPFN across folds and across model families, while rainfall is stable on both axes\.

The basin\-level maps in Figures[13](https://arxiv.org/html/2608.05265#S4.F13)and[14](https://arxiv.org/html/2608.05265#S4.F14)let us read the spatial structure of the SHAP field directly against the spatial structure of each feature field\. Because the features themselves are spatially autocorrelated within a fire \(Section[2\.1](https://arxiv.org/html/2608.05265#S2.SS1)\), this comparison is a non\-trivial test of whether the model is conditioning on the feature rather than on a global bias\. If the model were ignoring a feature, the SHAP map for that feature would be spatially flat or unrelated to the feature map\. Instead, on the Grand Prix–Old fire each of the four Staley17\-style features produces a SHAP map whose high and low patches track the corresponding feature map \(Figure[4](https://arxiv.org/html/2608.05265#S2.F4)\), with basins of high storm intensity, highPropHM23, highdNBR, and highKFcarrying positive SHAP and the converse basins carrying negative SHAP\. On Station,dNBRandKFSHAP fields mirror the spatial variation of those feature fields, while the SHAP maps for peak 15\-minute intensity andPropHM23show weaker spatial correspondence\. This co\-variation between the feature map and the SHAP map is the spatial signature of the model conditioning its prediction on the physical quantity rather than on a global bias or on incidental basin identity\.

The dependence scatter plots in the bottom row of each figure make this relationship direct\. SHAP value is plotted against feature value, so a monotone trend indicates the model is mapping more of the feature into a larger contribution to debris\-flow probability\. On Grand Prix–Old all four features show clear monotonic SHAP–feature dependence, with positive Spearmanρ\\rhofor the rainfall, burn\-severity, and erodibility features, indicating that the model has internalized the expected sign for each feature with a clear physical direction\. On Station,dNBRandKFretain strong monotonic dependence \(ρ=0\.71\\rho=0\.71and0\.750\.75\), whilePeakI15andPropHM23flatten \(ρ=0\.27\\rho=0\.27–0\.450\.45\), showing that the contribution of those two features to the prediction varies little with their basin\-level values at this event\. Read together for these two fires, the agreement between the spatial structure of the feature and SHAP maps and the monotonic SHAP–feature scatter relationships indicates that TabPFN’s predictions are guided by feature–response relationships with the expected sign rather than by spurious associations\.

Beyond feature reliance, a meaningful concern for any data\-driven hazard prediction model is whether the learned relationships point in the physically expected direction\. Figure[11](https://arxiv.org/html/2608.05265#S3.F11)addresses this directly\. All four models assign positive SHAP contributions to higher values of burn severity, terrain steepness, soil erodibility, and every rainfall metric, with no reversals observed for any feature with a clear physical sign\. This concordance between learned and expected directions supports treating the top\-tier models not merely as pattern\-matching devices but as learning feature–response relationships with the expected sign for rainfall, burn severity, and soil erodibility, strengthening the case for their use alongside the operational Staley17 framework in the regional warning context this dataset was assembled for\.

### 4\.3Synthetic Data Augmentation Evaluation

![Refer to caption](https://arxiv.org/html/2608.05265v1/x14.png)Figure 15:Class\-conditional kernel density estimates of the benchmark features in a nine\-panel grid, with panels \(a\)–\(i\) shown separately for each feature and a shared bottom legend\. Solid curves are imputed real training observations and dotted curves are fold 0 synthetic observations from the TabPFN unsupervised generator, with teal denoting the no\-debris\-flow \(no DF\) class and orange denoting the debris\-flow \(DF\) class\. Inset tables report the two\-sample Kolmogorov–Smirnov statisticDD\(bounded in\[0,1\]\[0,1\]\) and the Wasserstein\-1 distanceW1W\_\{1\}between real and synthetic observations for each class\.![Refer to caption](https://arxiv.org/html/2608.05265v1/x15.png)Figure 16:Per\-model threat score for the unaugmented and augmented evaluations as overlaid box plots \(left axis\) and the mean change in threat score from unaugmented to augmented as bars \(right axis\)\. Models are ordered along the horizontal axis by descending unaugmented threat score to match the order in previous figures\.The synthetic data augmentation evaluation shows that TabPFN\-generated synthetic rows can provide a measurable benefit for many models, but the effect is not uniform across models\. TabPFN serves two roles in this evaluation, as the generator of synthetic rows and as one of the evaluated models, which raises the possibility that augmentation could either favor TabPFN with synthetic rows that match its expected distribution or, conversely, leave it unchanged because TabPFN can already represent that distribution from the unaugmented training data alone\. The KDE and t\-SNE realism checks \(Figures[15](https://arxiv.org/html/2608.05265#S4.F15)and[12](https://arxiv.org/html/2608.05265#S3.F12)\) confirm the synthetic rows are distributionally close to real, which narrows the dual\-role question to whether TabPFN benefits more than other models from data it can already represent\.

The observed pattern is that TabPFN does not benefit\. TabPFN gained only\+0\.004\+0\.004threat score while several non\-TabPFN models gained substantially more \(Figure[16](https://arxiv.org/html/2608.05265#S4.F16)\), which argues against a generator\-model bias and is instead consistent with TabPFN already saturating the dataset’s information\. That said, per\-feature realism is uneven in a way that bears on interpretation: soil erodibility \(KF\) carries the largest realism gap while rainfall features remain close, andKFis also the feature whose SHAP ranking diverges most across models \(rank 3–4 for XGBoost and CatBoost versus 10 for Random Forest\), so any model that relies on a sharpKFthreshold may be the most exposed to residual generator error\. Such model\-specific sensitivity to augmentation quality is expected in imbalanced classification, where resampling can improve minority sensitivity yet still interact with model bias and decision\-boundary geometry\(chawla\_smote\_2002;he\_learning\_2009;branco\_survey\_2016\)\.

Gains in the augmented top tier \(Table[5](https://arxiv.org/html/2608.05265#S3.T5)\) are marginal and uneven\. CatBoost, Random Forest, and XGBoost each gained roughly two percentage points in threat score, while TabPFN and ExtraTrees moved by less than half a point, consistent with the unaugmented top tier already approaching the dataset’s predictive ceiling\.

Practically, synthetic data augmentation should be treated as a tunable component of the pipeline rather than a universal improvement step\. The strongest deployment strategy is to validate synthetic data augmentation choices on held\-out real folds, report delta metrics, and retain only settings that improve both event detection and error trade\-offs for the intended warning objective\.

## 5Conclusion

This study presented a unified evaluation of 15 models spanning five model families for post\-wildfire debris\-flow prediction across the western United States, to our knowledge, the first benchmark for this problem to include the new foundation model TabPFN, evaluate SHAP feature importance, and evaluate synthetic data augmentation\. TabPFN and the leading tree\-based models \(ExtraTrees, CatBoost, XGBoost, Random Forest\) clustered tightly at threat score 0\.62–0\.64 against 0\.15 for the Staley17 model, a roughly fourfold gap, with TabPFN matching the tuned models while requiring no dataset\-specific hyperparameter search\. SHAP analyses across the four explained models converged on the same three highest\-ranked and most fold\-stable rainfall features \(PeakI15,StormAccum, andPeakI30\), indicating that these model families emphasized the same predictive signals rather than divergent ones\. The four features Staley17 used to define its logistic regression model mapped only partially onto the SHAP rankings of the modern models\. Basin\-level SHAP maps for two illustrative fires further showed spatial SHAP fields that tracked the corresponding feature fields and monotonic SHAP–feature dependence scatter relationships, indicating that TabPFN conditioned its predictions on feature–response behavior with the expected sign at those events rather than on spurious associations\. TabPFN\-generated synthetic data augmentation, applied to mitigate the small, class\-imbalanced training data, yielded the largest gain for the most restricted model, Staley17, 13\.3 percentage points in threat score, while the top nonlinear tier improved by at most 2\.2 percentage points, indicating that the benefit of synthetic augmentation was governed primarily by feature\-set restriction and model capacity rather than acting as a universal enhancement\. These findings were most directly applicable to datasets with similar spatial and temporal scope, because model rankings were estimated under random stratified cross\-validation and the data were concentrated in southern California within the first two post\-wildfire years\. Repeated cross\-validation strengthened ranking reliability relative to a single split, but it did not replace the need for stricter spatial transfer tests when models are deployed across new fire regions\. Future work should extend evaluation to multi\-region datasets with strict spatial holdouts, longer post\-wildfire temporal windows, and explicit threshold\-optimization and cost\-sensitive objectives that reflect the asymmetric consequences of missed events versus false alarms\.

## References