Empty Shell Syndrome

Rapid surface greening conceals arrested woody biomass recovery and elevates secondary mortality Nature Plants (2025)

Multi-sensor satellite evidence that post-disturbance canopy greenness masks biomass stagnation. DTW clustering (K=3) identifies the Empty Shell syndrome in 183 of 601 training sites. XGBoost (AUC=0.909) + SHAP attribution identifies maximum temperature, elevation, and burn severity as dominant drivers. Cox PH survival analysis: HR=1.60 (P=1.61e-10). Global extrapolation: ~2.65 Pg C hidden carbon shortfall.

---

Repository Structure

empty-shell-syndrome/
├── README.md                          # This file (EN/CN)
├── README_CN.md                       # Chinese version
├── config.py                          # Auto-resolves data/ and temp/ paths
├── requirements.txt                   # Python dependencies
│
├── data/                              # Immutable input data (15 files, ~512MB)
│   │                                   #   GEE-extracted satellite + climate records
│   │                                   #   Do NOT edit — scripts read from here
│   └── XGBoost_Global_20F.json        #   Trained binary classifier (paper model)
│
├── temp/                              # Runtime outputs (auto-cleaned)
│   │                                   #   All pipeline/validation/figure outputs
│   │                                   #   Safe to delete — regenerated by scripts
│
├── pipeline/                          # Main experiments (Methods 5.1–5.5)
│   ├── clean_sites.py                 # 5.1: ESA+WRI+land-cover filter (1613→1393)
│   ├── assemble_panel.py              # 5.1: Align Landsat + merge 6 data layers
│   ├── pipeline.py                    # 5.2–5.4: [CORE] Full Yan pipeline: DTW→XGBoost→Cox
│   ├── XGBoost.py                     # 5.3: Binary XGBoost classifier (ES vs Others)
│   ├── SHAP.py                        # 5.3: SHAP TreeExplainer feature ranking
│   ├── Empty_Shell_probability.py     # 5.3: Predict ES probability for 1,393 sites
│   ├── Cox_PH.py                      # 5.4: Nested Cox PH + Schoenfeld + AIC
│   ├── survival.py                    # 5.4: Kaplan-Meier + log-rank test
│   └── global_carbon_sink.py          # 5.5: Global 5km projection + carbon deficit
│
├── validation/                        # Robustness checks (7 scripts)
│   ├── ablation_single_metrics.py     # Drop one satellite metric at a time
│   ├── ablation_biome_type.py         # Remove Biome_Type, retrain (ΔAUC)
│   ├── window_3-10_stability.py       # W=3–10 observation window sweep
│   ├── validate_dtw_2-8__clusters.py  # DTW cluster quality (Silhouette + CH)
│   ├── spatial_thinning.py            # Grid-based thinning + criterion ES
│   ├── tertile_mortality_gradient.py  # ES probability tertile mortality
│   └── compute_all_pvalues.py         # All p-values reported in paper
│
├── figures/                           # Paper figures (19 scripts)
│   ├── fig1a_es_probability_map.py    # Fig 1a: Global 1,393-site ES probability
│   ├── fig1b_recovery_trajectories.py # Fig 1b: 4-metric group trajectories (95% CI)
│   ├── fig1b_site118_zscore.py        # Fig 1b inset: Site 118 Z-score anomaly
│   ├── fig1b_site118_raw.py           # Fig 1b inset: Site 118 raw physical values
│   ├── fig1c_violin_endstate.py       # Fig 1c: Year-8 4-metric violin + Welch's t
│   ├── fig2_shap_drivers.py           # Fig 2a–h: SHAP importance + beeswarm + 6 dependence
│   ├── fig3_kaplan_meier.py           # Fig 2 (KM): Kaplan-Meier survival curves
│   ├── fig3_cox_forest.py             # Fig 2 (Cox): Nested Cox forest plot
│   ├── fig4_cartopy_render.py         # Fig 3a,b: Cartopy vector maps
│   ├── figS1a_k_trend.py              # Fig S1a: K=2–5 decoupling trajectories
│   ├── figS1b_k_auc_hr.py             # Fig S1b: K vs AUC & HR
│   ├── figS1c_k_spatial_map.py        # Fig S1c: K=2–5 spatial distribution
│   ├── figS2a_window_hr.py            # Fig S2a: Observation window HR stability
│   ├── figS2_roc_pr.py                # Fig S2a/b: ROC + PR curves
│   ├── figS2_regional_prevalence.py   # Fig S2a: Regional ES prevalence
│   ├── figS2_raw_data_check.py        # Fig S2b: Raw NDVI/VOD/SPEI reality check
│   ├── figS2_biome_mortality.py       # Fig S2c: Biome-level Fisher analysis
│   ├── figS2_feature_ablation.py      # Fig S2d: Feature ablation SHAP
│   └── figS2_training_distribution.py # Training set biome composition
│
└── gee/                               # Google Earth Engine extraction (8 scripts)
    ├── Hansen GFC+WRI forest loss driver+ESA WorldCover 2021.py
    ├── Landsat_NDVI_NDII.py
    ├── SPEI.py
    ├── GOSIF_SIF.py
    ├── VODCA_VOD.py
    ├── ABG_canopy_height.py
    ├── 20_environmental_features.py
    └── global 5km multi-band GeoTIFF.py

---

Quick Start

pip install -r requirements.txt
cd pipeline
python clean_sites.py                # Step 1: filter sites
python assemble_panel.py             # Step 2: build master panel
python pipeline.py                   # Step 3: run full pipeline

All scripts auto-detect data/ for input and temp/ for output via config.py.

---

Data Pipeline

1,613 raw sites → ESA+WRI+land-cover filter → 1,393 analysis sites
  → Yan pipeline (NDVI-minimum onset, Forest+No-repeat+≥5 obs/metric)
  → 601 DTW training sites, K=3 clustering identifies Empty Shell (max NDVI-VOD decoupling)
       ES = 183, Others = 418

601 DTW labels → XGBoost binary (20 features, AUC=0.909±0.026)
  → Predict 1,393 sites → Prob_ES > 0.5 = 511 (36.7%)
  → Cox PH: HR=1.60 [1.39,1.85], P=1.61e-10
  → Global extrapolation: 22.2% forest at risk, ~2.65 Pg C deficit

Key Results

Analysis	Method	Result
DTW training set	Yan pipeline K=3	601 sites, 183 ES
XGBoost	5-fold CV, max_depth=2	AUC = 0.909 ± 0.026
SHAP top features	TreeExplainer	tmmx (0.70), elevation (0.28), dNBR (0.21)
Cox PH	ES alone (Z-standardized)	HR = 1.60 [1.39, 1.85], P = 1.61e-10
Cox PH + Sand Content	Mediation test	HR = 1.10, P = 0.22 (ns)
Kaplan-Meier	Log-rank test	χ² = 28.52, P = 9.29e-08
Global pixel prediction	5km raster, 25.8M pixels	22.2% Prob_ES > 0.5
Phantom carbon deficit	Σ Prob×δ×AGB×Area	~2.65 Pg C

Requirements

Python 3.10+, RAM ≥ 16GB (32GB recommended for Fig 4). See requirements.txt.

References

• Hammond et al. (2022) Nature Communications 13, 1761.
• Yan et al. (2025) Nature Plants 11, 731-742.
• Forzieri et al. (2022) Nature 608, 534-539.