# Calibrated subseasonal-to-seasonal probability forecasts of extreme-heat days over Vietnam from a single-member reforecast *Draft manuscript prepared for* Monthly Weather Review*. Prose revised 2026-10-12 to the authors' established academic style (AMS/American spelling). Every quantitative claim is traceable to a named file in the project Evidence pack (`docs/RESULTS_REPORT.md` §13). Author list, affiliations, funding, and final reference details are placeholders to be completed by the authors; the reference list (§References) contains real, citable works to be finalized by the authors against journal style.* **Authors:** [A. Author]¹, [B. Author]², … · ¹[Affiliation], ²[Affiliation] **Corresponding author:** [name, email] --- ## Abstract Extreme heat is among the most consequential hazards in Vietnam, yet skillful operational guidance at the subseasonal-to-seasonal (S2S) range of days to six weeks ahead remains scarce. The problem is sharpened here by a reforecast archive that contains only a single deterministic (control) member. In this study, we build, calibrate, and verify a probabilistic forecast of extreme-heat days (EHDs) over Vietnam to a lead of 45 days, using a 20-yr (2003–2022) ECMWF S2S control reforecast. An EHD is defined per grid cell and calendar day as daily maximum temperature `Tx ≥ min(Tx90p, 35 °C)`, and the forecast target, the extreme-heat-day event (EHDE) defined as `P(≥1 EHD in a block)`, is issued for six fixed, non-overlapping lead blocks (W1 0–10 d to W6 39–45 d). Because no ensemble exists, probabilities are produced by single-member model output statistics (MOS), in which a family of statistical and machine-learning models maps the deterministic forecast and physically motivated covariates to a calibrated probability rather than to an ensemble exceedance fraction. All skill is estimated by leave-one-year-out (LOYO) cross-validation, with every fitted quantity strictly excluded from the held-out year and independently re-derived in plain NumPy as an integrity check. The results show that, using the forecast temperature alone, the calibrated product has statistically significant Brier skill over climatology to week 3 (day 24) on the national average, and longer in the tropical south and in the pre-monsoon (MAM) season. A leakage-safe feature-selection analysis indicates that, once the forecast-Tx signal is exhausted at extended range, slowly varying predictors carry the residual predictability: the El Niño–Southern Oscillation (ENSO) Niño-3.4 index emerges as the single most important predictor at weeks 3–6, ahead of the Madden–Julian Oscillation (MJO), the 500-hPa flow, and antecedent persistence, whereas land-surface memory adds little. Supplying these predictors with out-of-fold calibration and stacking extends significant skill from week 3 to week 6. These results offer opportunities for developing an operational S2S heat early-warning capability for Vietnam. **Keywords:** extreme heat; heat waves; subseasonal-to-seasonal prediction; forecast calibration; model output statistics; machine learning; verification; Vietnam. --- ## 1. Introduction ### a. Extreme heat, its impacts, and the Vietnamese setting Prolonged and intense hot spells cause excess mortality, strain electricity systems, and depress agricultural and labor productivity across Southeast Asia, and their frequency, duration, and intensity have risen with anthropogenic warming and are projected to intensify further (Perkins-Kirkpatrick and Lewis 2020; Perkins 2015; *Heatwaves in Southeast Asia and their changes in a warmer world*, 2021). The heat in the region is increasingly humid, which compounds its physiological impact (Li 2020), and the 2023 Southeast Asia event underscored both the severity and the limited predictability of such extremes (*The characterization, mechanism, predictability, and impacts of the unprecedented 2023 Southeast Asia heatwave*). Vietnam is particularly exposed. It spans roughly 8.5°–23.5°N over seven climatologically distinct regions, from the subtropical, monsoon-modulated north to the persistently hot tropical south, and it experiences extreme heat with pronounced seasonal and latitudinal structure: a pre-monsoon (March–May, MAM) peak in the north and center, driven by subsidence and adiabatic warming ahead of the summer monsoon, contrasted with a more time-invariant tropical regime in the south. Over the wider Indochina Peninsula, the synoptic controls, long-term trends, and their modulation by ENSO are now reasonably well characterized (*Synoptic characteristics … of heat waves over the Indochina Peninsula*; *Impacts of El Niño–Southern Oscillation on heat waves in the Indochina Peninsula*), and the human dimension in Vietnam, including spatial patterns of health vulnerability, occupational heat exposure, and the value of early-warning systems, is an active field (*Spatial patterns of health vulnerability to heatwaves in Vietnam*; *Heat index trends … for occupational heat exposure in Da Nang, Vietnam*; *Valuation of a heatwave early warning system … in Central Vietnam*). Heat preparedness, including the pre-positioning of health resources, the scheduling of grid maintenance and demand management, and the advising of outdoor labor and agriculture, requires lead times of two to six weeks. This is precisely the subseasonal-to-seasonal (S2S) range where raw dynamical skill for surface-temperature extremes is marginal and where statistical post-processing is indispensable (White et al. 2017; Robertson et al. 2020; Vitart et al. 2017). ### b. Predictability of heat extremes at the S2S range The S2S range bridges the deterministic weather limit (~10 days) and the seasonal range in which boundary conditions dominate, and it draws predictability from a distinct set of slowly varying sources: the Madden–Julian Oscillation (MJO; Madden and Julian 1971; Vitart 2014; Lim et al. 2018), the El Niño–Southern Oscillation (ENSO) and its modulation of tropical convection and the subtropical highs, quasi-stationary Rossby-wave patterns in the 500-hPa flow, and land-surface memory through soil moisture (Koster et al. 2010; Miralles et al. 2014). For heat specifically, several of these sources operate in tandem, through blocking and anticyclonic anomalies aloft, warm advection and subsidence in the lower troposphere, and antecedent surface dryness that shifts the surface energy balance toward sensible heating (Perkins 2015; Miralles et al. 2014). Regional S2S studies have quantified these sources and the resulting skill for heatwaves over the Yangtze River Basin (Xie et al. 2020; Pan et al. 2026), India (Mandal et al. 2025), West African cities, and Europe, and have shown that hot extremes can be *more* predictable than average summers at the subseasonal range and that land-surface uncertainty and soil-moisture coupling modulate that skill (*Higher subseasonal predictability of extreme hot European summer temperatures*; *Skillful subseasonal ensemble predictions of heat wave onsets through better representation of land surface uncertainties*; *Precipitation and soil moisture coupling constrains subseasonal predictability of a prolonged extreme heatwave*). Which source dominates is strongly lead-dependent. Near initialization the dynamical forecast still carries usable temperature information, whereas past two to three weeks that signal is largely exhausted and the residual predictability lives in the boundary and large-scale states. A forecast system for this range should therefore be designed to be lead-aware, and a central empirical question of this paper is how far, and from which physical sources, the skillful horizon for extreme heat over Vietnam can be pushed. ### c. Statistical post-processing and machine learning Raw dynamical output at the S2S range is biased and, when converted naively to probabilities, badly calibrated, so skillful probabilistic guidance almost always requires statistical post-processing (Wilks 2011; Vannitsem et al. 2021). The classical tool is model output statistics (MOS), a regression from model predictors to the observed variable or its probability. Modern practice augments linear MOS with machine-learning regressors and, increasingly, neural networks, which have been shown to match or exceed established ensemble post-processing on temperature and other fields (Rasp and Lerch 2018; Grönquist et al. 2021; Schulz and Lerch 2022; Vannitsem et al. 2021). For heat specifically, a rapidly growing literature applies random forests, gradient boosting, convolutional and deep networks, and hybrid dynamical–statistical schemes to subseasonal and seasonal heatwave prediction, often with explainable-AI diagnostics of the learned predictors, and reports skill gains over raw dynamical output (*Subseasonal prediction of central European summer heatwaves with linear and random forest machine learning models*; *Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data*; *Advancing subseasonal surface air temperature and heat wave prediction skill … in a deep learning model*; *Seasonal heatwave forecasting with explainable machine learning and remote sensing data*; *Decoding sub-seasonal predictors of extreme heat with interpretable machine learning*). We adopt this machine-learning-augmented MOS view, but we hold it to a strict out-of-sample standard and, critically, adapt it to the single-member setting. Two properties of the present problem shape the method. First, calibration and discrimination are distinct goals: a flexible model can be *sharper* yet *less reliable*, so probabilities must be recalibrated, here by isotonic regression (Zadrozny and Elkan 2002; Niculescu-Mizil and Caruana 2005), and evaluated with proper scoring rules that reward calibration and sharpness jointly (Gneiting and Raftery 2007; Gneiting et al. 2007). Second, because extreme heat is a rare, extremal event, verification must use metrics that do not degenerate to triviality as the base rate falls, namely the extremal-dependence indices EDI and SEDI (Ferro and Stephenson 2011) alongside the Brier score and its decomposition (Brier 1950; Murphy 1973). ### d. The single-member constraint Most probabilistic S2S calibration presumes an ensemble: probabilities begin as exceedance fractions over perturbed members and are corrected by ensemble MOS, and the ensemble supports spread-based uncertainty and rank-histogram diagnostics (Vannitsem et al. 2021). The reforecast archive available here is control-only, that is, a single deterministic member per initialization, which removes that entire machinery. This is not merely an inconvenience. It forces probabilities to be produced by single-member MOS, or distributional regression, in which the model learns a mapping from the deterministic forecast (plus covariates) to a calibrated event probability, and it makes the choice of predictand, model family, and calibration the load-bearing design decisions rather than the ensemble statistics. Comparatively little of the S2S post-processing literature addresses the single-member case for extremes directly, and doing so rigorously, without leaking information across the cross-validation boundary, is a contribution in its own right. ### e. Objectives and contributions Against this background, and driven by two concrete consumers, a duty forecaster who must run the system on a schedule to 45 days and a peer-reviewed audience that requires the method to be leakage-free and reproducible by someone who did not build it, this paper makes the following contributions. 1. A fixed, physically motivated extreme-heat-day index and a family of related predictands, namely the daily EHD, the block event EHDE, a severity variable `Sev`, and a heat-wave target, with a six-block lead-time verification protocol and an explicit block-straddle rule (§3a–b). 2. A single-member MOS framework and an expanded, lead-aware machine-learning model zoo (penalized logistic regression, three gradient-boosting implementations, random forests, multilayer perceptrons, a fully convolutional CNN, a residual U-Net, and a feature-tokenizer Transformer), all documented to the level required for reproduction (§3e) and all verified strictly out of sample. 3. A head-to-head comparison of three predictand formulations, namely event classification, severity regression with an isotonic probability link, and a two-level cascade that screens easy non-events, on a common leakage-safe footing (§3a, §3f, §4f). 4. A leakage-safe predictor-pool construction and selection over land-surface, local-thermodynamic, large-scale-circulation, and teleconnection/persistence candidates (§3c–d), which identifies the sources of extended-range predictability. 5. Two ensemble constructions, a leakage-safe stacked ensemble and a transparent 23-member vote ensemble, with an explicit ensemble-versus-best-single-member analysis that places the neural members on a common footing (§3g, §4g–h). 6. A region×week event product and a secondary heat-wave product (§4h–i), and a controlled negative result on true spatial deep learning (§4k). 7. An operational CF-compliant deliverable, an independent numerical re-derivation of the headline scores, and a fully reproducible, leakage-audited pipeline (§3h, §5, Data and code availability). The rest of the paper is organized as follows. Section 2 describes the data. Section 3 describes the indices, predictands, predictors, models, calibration, cross-validation, and verification. Section 4 presents the results, section 5 discusses them, and section 6 draws the conclusions. --- ## 2. Data ### a. S2S reforecast (ECMWF, control member only) We use the ECMWF S2S reforecast for the period 2003–2022 (20 yr), with 104 initializations per year on Thursdays and Sundays (2080 initializations in total). Only the control forecast is archived, as no perturbed members exist, which, as emphasized above, is the single most consequential fact for the method (§3e). Forecast daily maximum temperature is derived from the 6-hourly maximum-temperature field `mx2t6` as the daily maximum of its four 6-hourly values; the 46-step horizon (184 × 6 h) yields lead days 0–45, where day 0 is the first 24 h after initialization. The reforecast is delivered already bilinearly interpolated to the 0.25° observation grid (longitude 101°–111°E, 41 points; latitude 8°–24°N, 65 points), so no downscaling stage is applied, and the upstream interpolation is treated as a documented limitation (§5). Additional control fields used for the predictor extension (§2d) are 500-hPa geopotential, 850-hPa temperature, and surface and pressure-level moisture and wind, all on the same reforecast stream. Because the reforecast follows a `365_day` calendar with 29 February removed, forecast-derived Tx is matched to observations on (month, day). ### b. Observations: gridded target and station validation The calibration and verification target is the gridded observed daily `Tx` on the identical 0.25° grid (541 land cells), on a `365_day` calendar, covering the period 2002–2023, from which the climatology and the `Tx90p` threshold are constructed per fold (§3a). As an independent point-scale check we use 200 station daily-`Tx` records, of which 189 stations have usable overlap with the gridded field on common days. Reconciling the two reveals a systematic warm offset: stations read on average +1.5 to +2.2 °C warmer than their enclosing 0.25° cell in every region, with high day-to-day correlation (0.83–0.98), indicating a representativeness/offset gap rather than a timing error. Because the EHD threshold is derived from the cooler gridded field, station EHD frequency (0.31–0.59) runs 2–5 times the gridded frequency (≈0.10). We therefore fit and verify the product entirely in the gridded world and use the station comparison only to bound how far point-scale risk may exceed the gridded product (§5). Table 1 summarizes the data streams. **Table 1.** Principal data streams. | Stream | Source / variable | Grid / calendar | Period | Role | |--------|-------------------|-----------------|--------|------| | Reforecast Tx | ECMWF S2S control `mx2t6` → daily Tx | 0.25°, 41×65; 365-day | 2003–2022, 104 inits yr⁻¹ | Forecast predictor | | Reforecast fields | Z500, T850, surface & p-level moisture/wind (control) | 0.25° (VN) + 1.0° (extended) | 2003–2022 | Extended predictors | | Gridded obs Tx | Analysed `tasmax` | 0.25°, 541 land cells; 365-day | 2002–2023 | Target + climatology | | Station obs Tx | 200 station daily Tx (189 usable) | Point; daily | overlap | Independent validation | | Region mask | Integer `reg` mask (ids 1–7) | 0.25° | — | Stratification | | Teleconnections | ENSO Niño-3.4; MJO OMI (PC1, PC2, amplitude) | Monthly / daily indices | 2003–2022 | Extended predictors | ### c. Regions Seven regions are defined by the gridded integer `reg` mask (ids 1–7): Northwest (R1), Northeast (R2), Red River Delta (R3), North Central (R4), South Central (R5), Central Highlands (R6), and the South (R7). Region membership and extents are taken from the gridded mask; the administrative shapefile is used only for the publication map. All regional statistics are computed as zonal aggregates over the gridded mask, with no vector overlay, so they are internally consistent with the gridded target. ### d. Auxiliary predictor fields For the predictor extension (§3c) we assemble, from the same control reforecast and from observed teleconnection indices, four physically distinct groups: (i) land-surface state (soil moisture at 20 and 100 cm and deep soil temperature); (ii) local thermodynamics (skin temperature, CAPE, total column water, cloud cover, and dewpoint); (iii) large-scale circulation (500-hPa geopotential over a Vietnam/South-China-Sea box, over a subtropical box, and their gradient, together with box-mean 850-hPa temperature); and (iv) teleconnection/persistence (ENSO Niño-3.4 for the initialization month, the MJO OMI components PC1 and PC2 and amplitude at initialization, and antecedent 15-day observed EHD frequency and mean severity). A comprehensive gridded full-field cache (55 channels on an extended 1.0° domain and the native 0.25° grid; §3e, §4k) is additionally built for the true-spatial deep-learning experiment. All predictor construction is leakage-safe by design (§3c). --- ## 3. Methods ### a. Extreme-heat indices and the family of predictands *Extreme-heat day (EHD).* A day at a grid cell is an EHD if and only if > `Tx ≥ Thr`, with `Thr = min(Tx90p, 35 °C)`, where `Tx90p` is the calendar-day 90th percentile of observed daily `Tx` over the 2002–2023 baseline, estimated from a centered 15-day calendar-day window (±7 days) to stabilize the percentile against sampling noise. The `min(·)` rule is fixed by design and is a deliberate, physically motivated choice: it requires a day to be *both* locally unusual (at or above the local 90th percentile) *and* absolutely hot (capped at 35 °C), which prevents cool highland cells from labeling a merely above-normal but mild day as "extreme". The alternative `max(Tx90p, 35 °C)` variant is computed only as a diagnostic and is never substituted into the product. The percentile window is a *climatology-construction* device and is kept conceptually and operationally distinct from any forecast time-aggregation window; it is fit per fold, so the held-out year never informs its own threshold (§3h). *Extreme-heat-day event (EHDE), the primary predictand.* Because the operational decision is made per lead block, the forecast target is not the daily EHD but the block event > `EHDE = 1 ⇔ at least one EHD occurs within the block`, that is, `P(≥1 EHD in block)`. EHDE is the primary predictand throughout; unless stated otherwise, "the product" refers to calibrated `P(EHDE)`. We adopt the ≥1-EHD definition, rather than a stricter multi-day count, so that the event is neither too rare to verify nor trivially frequent, and so that all model variants remain directly comparable. *Severity (Sev), a regression predictand.* As an alternative pathway we define a continuous severity at the block scale, > `Sev = max_{d ∈ block} ( Tx_obs,d − Thr_d )`, the block-maximum amount by which the hottest day clears its threshold, where `Tx_obs,d` is the observed daily Tx and `Thr_d` is the day's threshold. By construction `Sev ≥ 0 ⇔ EHDE = 1`, so a regression that predicts `Sev` maps monotonically to `P(EHDE)` and additionally carries a "how extreme" magnitude that pure classification discards. Predicted severity is converted to a probability by an isotonic link fit on training data only (§3f). This formulation lets us test whether learning the ordered magnitude of heat is a better-posed problem than learning the binary event directly (§4f). *Heat-wave, a secondary predictand.* The secondary target is the heat-wave, `P(≥3 consecutive EHDs)`, detected on the continuous 0–45-day daily-EHD series and attributed to every block that a qualifying run overlaps, with an explicit block-straddle rule (a run spanning a block boundary is credited to both blocks). It is rarer than EHDE (base rate 0.11–0.16 versus 0.33–0.43) and is verified under the same protocol (§4i). ### b. Lead blocks and verification targets Forecasts are aggregated into six fixed, non-overlapping lead blocks (days since initialization): W1 0–10, W2 11–17, W3 18–24, W4 25–31, W5 32–38, W6 39–45. No running or sliding windows are used anywhere in aggregation or verification; this is a hard design constraint that keeps the six blocks statistically independent in lead and makes the verification protocol unambiguous. Every predictand in §3a is defined and scored on these blocks. ### c. Predictor construction and leakage safety For each initialization × land cell × block, a backbone of ten forecast-only features is computed: the block maximum, mean, and minimum of forecast daily `Tx` (`f_max`, `f_mean`, `f_min`); the fraction of block days with forecast `Tx ≥ Thr` (`f_exceed`); the margin by which the hottest forecast day clears the block-mean threshold (`f_maxanom`); a block ordinal encoding lead dependence (`lead_id`); the initialization-date seasonal harmonics `sin/cos(2π·doy/365)`; and grid latitude and longitude. In the severity pathway the forecast severity `Tx_fc − Thr` is added as a feature. Beyond the backbone we assemble a candidate pool of extended predictors spanning the four physical groups of §2d, all constructed to be leakage-safe: the gridded and large-scale reforecast fields are pure forecast quantities independent of the threshold; the teleconnection indices (ENSO, MJO) are *observed states known at initialization*; and the persistence predictors (`f_ante_ehdfreq`, `f_ante_sevmean`) apply the fold-specific threshold to observations dated strictly *before* initialization. No candidate uses any information dated at or after the block it predicts, and none uses the held-out year. ### d. Predictor selection Predictors are ranked per lead block by out-of-sample permutation importance, defined as the reduction in held-out-year AUC when a single predictor's values are randomly shuffled, aggregated over the 20 leave-one-year-out folds, and cross-checked against univariate mutual information. A predictor is *selected* for a block if its mean out-of-sample importance clears a noise floor estimated from shuffled controls; its across-fold variability is reported as a diagnostic but is *not* used as a filter, because episodic teleconnections are legitimately high-variance and would be unfairly penalized by a variance-based screen. The procedure recovers a textbook, lead-dependent predictor structure (§4j) and yields the fixed per-block lists used operationally (Table 2). The number of selected extended predictors is itself lead-dependent, from a short thermodynamic set at W1, to the full slow suite at W2–W4, narrowing to a few robust slow predictors at W6. **Table 2.** Selected extended predictors per lead block (added to the ten-feature backbone). `t850_reg`, `z500_reg`, `z500_sub`, `z500_grad` are box-mean 850-hPa temperature and 500-hPa geopotential over the Vietnam and subtropical boxes and their gradient; `enso_nino34` is monthly Niño-3.4; `mjo_pc1/pc2` are OMI components; `f_ante_*` are antecedent 15-day EHD frequency and mean severity. | block | selected extended predictors | |-------|------------------------------| | W1 | cape, t850_reg, f_ante_ehdfreq, f_ante_sevmean | | W2 | t850_reg, z500_reg, z500_sub, z500_grad, enso_nino34, mjo_pc1, mjo_pc2, f_ante_ehdfreq, f_ante_sevmean | | W3 | enso_nino34, mjo_pc1, mjo_pc2, t850_reg, z500_reg, z500_sub, z500_grad, f_ante_ehdfreq, f_ante_sevmean | | W4 | enso_nino34, mjo_pc1, mjo_pc2, t850_reg, z500_reg, z500_sub, z500_grad, f_ante_ehdfreq, f_ante_sevmean | | W5 | enso_nino34, mjo_pc2, t850_reg, z500_sub, z500_grad, f_ante_ehdfreq, f_ante_sevmean | | W6 | enso_nino34, z500_reg, z500_sub, f_ante_ehdfreq | ### e. The machine-learning model zoo Because the reforecast is control-only, `P(EHDE)` cannot be an ensemble exceedance fraction; every model below is a single-member MOS map from features to a calibrated probability. The zoo was built incrementally, as a two-model headline (v1), an expanded severity-aware zoo (v2), and a lead-aware hybrid zoo with deep learning (v3), and all members share a common `fit(X, y) → predict_proba(X)` interface, LOYO-nested training (§3h), and out-of-fold calibration (§3f). Two references are always carried: `raw` (the forecast block-exceedance fraction) and `climat` (the training base rate). *Tabular families.* Table 3 lists the tabular members with their key hyperparameters. Penalized logistic regression (standardized features, L2) is the interpretable statistical baseline and, as it turns out, the strongest single tabular family at extended range. Histogram gradient boosting (HGB) is the machine-learning workhorse; XGBoost and LightGBM add gradient-boosting *diversity* with different regularization and split strategies (Chen and Guestrin 2016; Ke et al. 2017). Random forests (Breiman 2001) and a scikit-learn multilayer perceptron (MLP) complete the tabular set. Each family exists in a classifier form (learning EHDE directly) and, where meaningful, a severity-regressor form (`ridge_sev`, `hgb_sev`, `mlp_sev`, `xgb_sev`, `lgbm_sev`; §3a, §3f). **Table 3.** Tabular model families and principal hyperparameters. All classifiers standardize inputs where relevant and are trained with nested hyperparameter selection inside each LOYO fold. | family | type | key hyperparameters | |--------|------|---------------------| | logistic | classifier | L2, `C = 1.0`, standardised, `max_iter = 1000` | | hgb | classifier | learning rate 0.06, `max_iter = 300`, log-loss, L2 leaf reg. | | rf | classifier | 300 trees, `min_samples_leaf = 100` | | mlp | classifier | hidden (64, 32), L2 `α = 1e-3`, early stopping | | xgb | classifier | 400 trees, depth 4, lr 0.05, histogram method, seeded | | lgbm | classifier | 500 trees, 31 leaves, lr 0.05, seeded | | ridge_sev | Sev regressor → isotonic | ridge `α = 1.0`, standardised | | hgb_sev | Sev regressor → isotonic | lr 0.06, `max_iter = 300`, `min_samples_leaf = 200` | | mlp_sev | Sev regressor → isotonic | hidden (64, 32), early stopping | | xgb_sev / lgbm_sev | Sev regressor → isotonic | as xgb / lgbm, regression objective | *Deep-learning members.* Four neural members are implemented in PyTorch (CPU), all trained with a class-imbalance-weighted binary cross-entropy (`pos_weight` set from the training base rate), a 10 % validation split for early stopping, and out-of-fold isotonic calibration. (i) A fully convolutional CNN (`TorchCNN`) operates on block-aggregated forecast maps: a stack of 3×3 convolutions (32→32→16→1 channels) with per-channel standardization from training statistics produces a per-cell probability field. (ii) A residual U-Net (`TorchResUNet`; Ronneberger et al. 2015) over the 41×65 Vietnam grid adds a two-level encoder/decoder with skip connections and residual convolution blocks (batch-normalized), padding the grid to a multiple of four and cropping the output back. (iii) A deep MLP (`TorchANN`, hidden [128, 64, 32], ReLU, dropout 0.2) is a more expressive tabular network. (iv) A feature-tokenizer Transformer (`TorchFTTransformer`; Vaswani et al. 2017) projects each numeric feature to a `d_model = 32` token, prepends a learnable CLS token, and applies two multi-head-attention encoder layers (4 heads, dropout); it is kept deliberately compact because the extended-range training signal is weak. Members (iii) and (iv) were evaluated but, being BSS-negative at long lead, are excluded from the production stack (§4j). *True-spatial full-field network.* To test whether resolving real large-scale spatial patterns, rather than broadcasting box-mean indices, adds skill, we build a comprehensive full-field cache (block-mean of the raw gridded fields over the six blocks, ocean retained, leakage-free: 55 channels comprising extended-domain Z500/T850/U200/V200 on a 1.0° 71×61 grid and native-grid surface and multi-level T/H/W/Q/U/V fields) and a two-branch fully convolutional network (`FullFieldResUNet`): a residual U-Net over the Vietnam grid, whose channels are the native full-field fields plus the rasterized backbone features, is fused with a small CNN encoder over the extended domain that global-average-pools to a context vector tiled back onto the Vietnam grid. This is the definitive test of the spatial-morphology hypothesis, and its (negative) result is reported in §4k. ### f. Probability calibration and the two-level cascade *Calibration.* A flexible model can be sharper yet less reliable, so every member's probabilities are recalibrated out of fold: an inner K-fold over the *training* years yields out-of-fold probabilities, on which an isotonic recalibration is fit (Zadrozny and Elkan 2002); the member is then refit on the full training years and pushed through the frozen calibrator. For the severity regressors the same mechanism supplies the `Sev → P(EHDE)` link: an isotonic map is fit on training `(predicted Sev, EHDE)` pairs, which is monotone by the `Sev ≥ 0 ⇔ EHDE` construction. The headline v1 product (`hgb_recal`) is HGB with isotonic recalibration. *Two-level cascade.* Motivated by the abundance of easy non-events, we test a two-level cascade. Level 1 (an HGB classifier trained on all training samples) screens samples it calls near-certain non-events, with `P(EHDE) ≤ 1 − cutoff` and the cutoff swept 0.90 → 0.95 → 0.99, and passes their level-1 probability straight through. Level 2 zoo models are trained only on the surviving, harder samples and predict the survivors; the final probability is `screened ? level-1 : level-2`. The sweep is kept out of sample. In practice the easy-negative screen proved skill-neutral (differences in the third decimal; §4f), so members are combined at cutoff 1.0 (no screening) in the production configuration, which is a reportable, if negative, methodological finding. ### g. Ensembling: stacking and the vote ensemble Two combinations of the zoo are formed. The first is a leakage-safe stacked ensemble (`ens_stack`): a logistic meta-learner is fit leave-one-year-out on the members' out-of-sample probabilities, so no member's in-sample fit ever informs the stack weights; an L1-sparse variant selects members automatically. An unweighted mean (`ens_mean`) is carried as a reference. The second is a transparent 23-member vote ensemble intended for operational interpretability: seven scikit-learn families × three cascade cutoffs, plus the CNN and deep ANN, each cast a binary vote for EHDE via its own leakage-safe per-region×block threshold; the issued probability is the vote fraction `m/23`, and the deterministic flag fires above a per-fold Peirce-optimal `P₀`. Aggregating the flags to a region×week event, which alerts a region when the fraction of its cells voted EHD exceeds a coverage threshold, yields the most operationally robust product (§4h). ### h. Cross-validation and the leakage rule All skill is estimated by leave-one-year-out (LOYO) over the period 2003–2022, with hyperparameter selection nested inside each fold. For a held-out year *y*, none of the following may use data from *y*: (1) the observed climatology and the `Tx90p` threshold; (2) the calibration model and its hyperparameters; (3) the probability recalibration; and (4) the per-region×block decision threshold. Each stage records the years it was fit on, and a verification stage *asserts* that the held-out year is absent from every fit-year set before scoring; a threshold or climatology fit on data used for verification is treated as a failed deliverable, not as a rounding issue. A separate all-years refit is used for operational deployment, but the *honest* skill is always the LOYO estimate, never the in-sample fit. As an integrity check, an independent validator re-derives the headline per-block Brier score, BSS, AUC, and bootstrap confidence intervals from the raw prediction files in plain NumPy, without reusing the main verification code, and confirms agreement (PASS). ### i. Verification metrics and significance Forecasts are scored against *both* references, climatology and the raw forecast. The probabilistic metrics are the Brier score and Brier skill score, with its reliability/resolution decomposition (Brier 1950; Murphy 1973), reliability curves, ROC/AUC, and sharpness. The binary metrics at the decision threshold are the full ten-score inventory, namely the probability of detection (POD), probability of false detection (POFD), false-alarm ratio (FAR), critical success index (CSI), equitable threat score (ETS), Heidke skill score (HSS), Peirce skill score, and the extremal-dependence indices EDI and SEDI (Ferro and Stephenson 2011), together with frequency bias. Because EHDE is extremal, EDI and SEDI are the primary discrimination metrics, as ordinary skill scores degrade toward triviality as the base rate falls. Scores are stratified by lead block, region, initialization season, and climatological severity. Significance uses a whole-year block bootstrap (resampling entire years to preserve serial correlation); field significance across the 541 cells is controlled with the Benjamini–Hochberg false-discovery-rate (FDR) procedure at `q = 0.05` (Benjamini and Hochberg 1995). Rank-histogram and spread–error diagnostics are not applicable to a single member and are reported as such. --- ## 4. Results ### a. Headline cross-validated skill (primary predictand) Pooled over all 20 held-out folds (~1.12 million out-of-sample forecasts per block) and scored against a per-fold base-rate climatology, the calibrated `hgb_recal` product has significant Brier skill to W3 and none beyond on the national-pooled average (Table 4). Discrimination (AUC) stays above 0.5 at all leads, but the calibrated Brier skill reaches zero by W4. Calibration always beats the raw forecast (BSS versus raw +0.28 to +0.52), because the raw block-exceedance probability is badly miscalibrated. This is the central empirical justification for MOS in the single-member setting. **Table 4.** Headline LOYO skill, primary predictand (EHDE), `hgb_recal`. 95 % CI from whole-year block bootstrap. | block | days | base rate | BSS vs clim [95 % CI] | AUC | BSS vs raw | significant | |-------|------|-----------|-----------------------|-----|-----------|-------------| | W1 | 0–10 | 0.431 | 0.356 [0.332, 0.383] | 0.845 | 0.517 | yes | | W2 | 11–17 | 0.331 | 0.118 [0.092, 0.144] | 0.710 | 0.310 | yes | | W3 | 18–24 | 0.331 | 0.037 [0.009, 0.063] | 0.630 | 0.288 | yes (marginal) | | W4 | 25–31 | 0.331 | 0.014 [−0.012, 0.039] | 0.603 | 0.281 | no | | W5 | 32–38 | 0.331 | 0.004 [−0.021, 0.027] | 0.591 | 0.280 | no | | W6 | 39–45 | 0.331 | −0.006 [−0.028, 0.018] | 0.579 | 0.278 | no | ### b. Binary contingency inventory At the operational per-region×block decision threshold, the full ten-score binary inventory for `hgb_recal` decays monotonically with lead, and the frequency bias grows past 1 at long lead (over-forecasting), consistent with the BSS collapse (Table 5). The extremal indices EDI and SEDI, the appropriate discrimination metrics for this rare event, fall from ~0.68 at W1 to ~0.15 at W6. **Table 5.** Pooled binary inventory by block, `hgb_recal`. | block | base | POD | POFD | FAR | CSI | ETS | HSS | Peirce | EDI | SEDI | freq bias | |-------|------|-----|------|-----|-----|-----|-----|--------|-----|------|-----------| | W1 | 0.431 | 0.761 | 0.237 | 0.291 | 0.580 | 0.351 | 0.520 | 0.524 | 0.682 | 0.682 | 1.07 | | W2 | 0.331 | 0.644 | 0.345 | 0.520 | 0.379 | 0.160 | 0.275 | 0.299 | 0.415 | 0.417 | 1.34 | | W3 | 0.331 | 0.583 | 0.402 | 0.582 | 0.322 | 0.090 | 0.165 | 0.181 | 0.257 | 0.259 | 1.39 | | W4 | 0.331 | 0.557 | 0.414 | 0.600 | 0.303 | 0.070 | 0.130 | 0.143 | 0.203 | 0.205 | 1.39 | | W5 | 0.331 | 0.549 | 0.428 | 0.612 | 0.294 | 0.058 | 0.109 | 0.121 | 0.171 | 0.173 | 1.41 | | W6 | 0.331 | 0.521 | 0.415 | 0.617 | 0.283 | 0.051 | 0.097 | 0.106 | 0.148 | 0.152 | 1.36 | ### c. Reliability and the value of calibration The gap between the calibrated and raw products (Table 4, "BSS vs raw") is a direct measure of the value of post-processing: the raw block-exceedance probability is systematically overconfident, and isotonic recalibration removes most of the reliability penalty while preserving resolution. This is the expected behavior for a single deterministic member converted to a probability, and it confirms that the calibration step, rather than any single model choice, is what makes the product usable. ### d. Spatial and regional structure Per-cell Brier skill with FDR control (`q = 0.05`) is significant for 539/504/192/116/68/30 of the 541 land cells at W1 to W6, respectively. Regionally, the tropical south retains skill markedly longer than the north: the Central Highlands (R6) and the South (R7) show BSS ≈ 0.18–0.19 at W2, ≈ 0.10–0.12 at W3, and still positive values at extended range (≈ +0.07–0.08 at W4, decaying to ≈ +0.04 by W6), whereas R1– R3 in the north lose skill by W3 and turn slightly negative at W3–W6. The operational skillful horizon is thus lead- *and* region-dependent (north ≈ W2, south ≈ W3+), which is consistent with tropical heat being more persistent and less synoptically driven than its northern counterpart. ### e. Season and severity stratification Stratified by initialization season (Table 6), the pre-monsoon MAM season is the most predictable at extended range, and is the only season with positive BSS through W6, whereas DJF and JJA lose skill by W3. This is operationally favorable: the product is most trustworthy at long lead precisely in the season when extreme heat is most consequential. Stratified by climatological EHD frequency (terciles), skill rises with the base rate at W1 but converges to near-zero across all terciles by W4–W6. **Table 6.** BSS vs climatology by initialization season, `hgb_recal`. | block | DJF | MAM | JJA | SON | |-------|-----|-----|-----|-----| | W1 | 0.420 | 0.407 | 0.263 | 0.343 | | W2 | 0.117 | 0.177 | 0.076 | 0.104 | | W3 | 0.012 | 0.075 | 0.029 | 0.031 | | W4 | −0.010 | 0.051 | −0.004 | 0.019 | | W5 | −0.014 | 0.039 | −0.005 | −0.003 | | W6 | −0.021 | 0.025 | −0.011 | −0.016 | ### f. Predictand formulation: event classification, severity regression, and the cascade Comparing the three predictand formulations of §3a on the identical leakage-safe pipeline, we find that the severity-regression pathway (the `*_sev` members with an isotonic link) is competitive with direct event classification but does not systematically beat it: at short lead the boosting classifiers and severity regressors are within a few thousandths of BSS/SEDI of one another, and at long lead the linear/logistic member on the slow predictors leads regardless of whether it is posed as classification or regression. The two-level cascade that screens easy non-events is skill-neutral, as sweeping the screening cutoff from 0.90 to 0.99 moves pooled scores only in the third decimal, so the production configuration combines members without screening (cutoff 1.0). The practical conclusion is that the event framing is adequate; the severity variable is retained because it adds an interpretable magnitude and diversifies the ensemble, and not because it lifts the headline skill. ### g. Expanded zoo, and ensemble versus the best single member Scoring *every* member, including the CNN and the deep ANN, on the *identical* leakage-safe pipeline as the ensembles (per-held-out-year, per-region Peirce decision threshold) makes the ensemble-versus-member comparison apples-to-apples. On that common footing, the best single member is the fully convolutional CNN at every block. The stacked ensemble `ens_stack` is the best forecast at short lead (W1–W2), but at long lead (W3–W6) the CNN matches or slightly exceeds it (Table 7). Stacking therefore helps most where the signal is strongest; at long lead no combination beats the single CNN, a caution against assuming that ensemble averaging always helps once members become correlated and the signal is weak. **Table 7.** Best single member (CNN) vs stacked ensemble, SEDI / BSS / AUC (same pipeline). Bold = better value for that metric and lead. | block | CNN (best member) | `ens_stack` | |-------|-------------------|-------------| | W1 | 0.707 / 0.385 / 0.859 | **0.715 / 0.391 / 0.864** | | W2 | 0.445 / 0.136 / 0.722 | **0.454 / 0.139 / 0.725** | | W3 | **0.292 / 0.063 / 0.647** | 0.291 / 0.059 / 0.642 | | W4 | 0.245 / **0.046 / 0.622** | **0.249** / 0.041 / 0.616 | | W5 | 0.190 / **0.033 / 0.606** | **0.203** / 0.029 / 0.596 | | W6 | **0.181 / 0.027 / 0.593** | 0.170 / 0.022 / 0.581 | ### h. Vote ensemble and region×week event forecast A transparent 23-member vote ensemble roughly matches the best single member and sits a little below the trained `ens_stack` at short lead (W1 SEDI 0.674 versus 0.715), with the Peirce-optimal `P₀` near a simple majority (≈ 0.5). Aggregating the forecast to a region×week event, which alerts a region when the fraction of its cells forecast to be EHD exceeds a threshold, is markedly more skillful than the per-cell product (W1 Peirce 0.65 / SEDI 0.80, still positive at W6), because cell-scale noise averages out. The per-week Peirce-optimal coverage threshold is unstable across folds, so a fixed 0.5 majority reference is recommended operationally. ### i. Secondary predictand: heat-wave (≥3 consecutive EHDs) Verified under the same LOYO protocol, the heat-wave product (base rate 0.11–0.16, rarer than the primary) shows the same skill horizon as the primary: significant BSS versus climatology to W3 (0.269 / 0.092 / 0.033, all bootstrap CIs above zero) and climatology-equivalent at W4+ (Table 8). Discrimination is weaker than the primary at matched lead, and the binary heat-wave flag over-forecasts strongly (frequency bias 2.0–2.8): a Peirce-optimal threshold on a rare event buys hit rate at the cost of false alarms, so the flag should be read as a screening alert alongside the probability. **Table 8.** Heat-wave (secondary) LOYO skill, `hgb_recal`. | block | base rate | BSS vs clim [95 % CI] | AUC | SEDI | freq bias | significant | |-------|-----------|-----------------------|-----|------|-----------|-------------| | W1 | 0.158 | 0.269 [0.240, 0.300] | 0.849 | 0.691 | 2.00 | yes | | W2 | 0.136 | 0.092 [0.058, 0.128] | 0.727 | 0.453 | 2.63 | yes | | W3 | 0.136 | 0.033 [0.002, 0.059] | 0.649 | 0.302 | 2.59 | yes (marginal) | | W4 | 0.136 | 0.011 [−0.023, 0.041] | 0.614 | 0.223 | 2.57 | no | | W5 | 0.135 | 0.010 [−0.021, 0.037] | 0.612 | 0.221 | 2.66 | no | | W6 | 0.112 | −0.005 [−0.033, 0.022] | 0.594 | 0.188 | 2.82 | no | ### j. Extended predictors extend the skillful horizon The feature-selection analysis recovers a textbook subseasonal predictability transition *out of sample*. At short lead (W1–W2) the forecast-Tx statistics and box-mean 850-hPa temperature dominate. At extended range (W3–W6) the forecast's own temperature signal is exhausted and slowly varying predictors take over: ENSO Niño-3.4 is the single most important predictor at W3–W6, followed by the MJO, the 500-hPa circulation, the 850-hPa temperature, and antecedent persistence. Notably, soil moisture and soil temperature never clear the importance noise floor at any lead; in this control reforecast, land-surface memory adds little once ENSO, circulation, and persistence are present, a result that contrasts with land-atmosphere-coupling studies elsewhere and merits attention (§5). Adding the selected predictors naively raises discrimination at extended range (AUC +0.02 to +0.05, SEDI +0.05 to +0.12 at W3–W6) but degrades the Brier skill score, producing a sharper-but-less-reliable model, symptomatic of in-sample calibration. Recalibrating out of fold and stacking the backbone and the selected models converts the discrimination gain into calibrated skill: the stacked model beats the forecast-only backbone at every lead, with the whole-year bootstrap on the gain excluding zero throughout (Table 9). **Table 9.** Out-of-sample skill: forecast-only backbone vs the out-of-fold stacked model (with selected extended predictors). | block | base BSS | stack BSS | stack−base gain [95 % CI] | AUC base→stack | SEDI base→stack | |-------|----------|-----------|---------------------------|----------------|-----------------| | W1 | 0.348 | **0.377** | +0.029 [0.025, 0.034] | 0.842→0.856 | 0.674→0.700 | | W2 | 0.115 | **0.142** | +0.027 [0.014, 0.042] | 0.707→0.730 | 0.411→0.460 | | W3 | 0.039 | **0.072** | +0.033 [0.015, 0.049] | 0.622→0.658 | 0.242→0.307 | | W4 | 0.022 | **0.051** | +0.029 [0.014, 0.043] | 0.582→0.639 | 0.173→0.258 | | W5 | 0.015 | **0.053** | +0.038 [0.020, 0.054] | 0.570→0.639 | 0.143→0.274 | | W6 | 0.011 | **0.034** | +0.023 [0.011, 0.036] | 0.535→0.608 | 0.101→0.215 | Crucially, an absolute whole-year bootstrap shows that the stacked model has significant Brier skill over climatology at every block through W6 (W4 0.051 [0.028, 0.069], W5 0.053 [0.030, 0.070], W6 0.034 [0.017, 0.048]), whereas the Stage-1 headline lost significance at W4. Two separable effects contribute, in that the out-of-fold calibration lifts even the backbone marginally, and the selected extended predictors add substantially on top; together they extend significant extreme-heat probability skill from week 3 to week 6, with the long-lead skill physically sourced from ENSO, the MJO, the large-scale circulation, and persistence. ### k. A hybrid zoo, an error-diverse stack, and a controlled negative result on spatial deep learning Letting the full model zoo consume the lead-dependent predictor slice, which adds gradient-boosting diversity (XGBoost, LightGBM) and true spatial deep learning (a residual U-Net alongside the CNN) and combines all thirteen calibrated members with a leakage-safe logistic stack, yields a stack with positive Brier skill at every lead that equals or beats the best single member throughout, extending the horizon further at long lead (W4 0.051→0.071, W5 0.053→0.062, W6 0.034→0.058) and adding a substantial short-lead gain (W1 0.377→0.431; Table 10). **Table 10.** Thirteen-member calibrated stack, out-of-sample BSS versus climatology, with the best single member at each lead. | block | base rate | best single member | ens_stack BSS | ens_mean BSS | ens_stack SEDI | |-------|-----------|--------------------|---------------|--------------|----------------| | W1 | 0.431 | resunet 0.421 | **0.431** | 0.383 | 0.747 | | W2 | 0.331 | xgb / cnn 0.140 | **0.157** | 0.140 | 0.487 | | W3 | 0.331 | logistic 0.076 | **0.087** | 0.078 | 0.356 | | W4 | 0.331 | logistic 0.066 | **0.071** | 0.062 | 0.320 | | W5 | 0.331 | logistic 0.058 | **0.062** | 0.058 | 0.323 | | W6 | 0.331 | logistic 0.053 | **0.058** | 0.053 | 0.296 | The single-member ranking is itself a physical result. At short lead the flexible members win: gradient boosting leads the tabular families (W1 BSS lgbm 0.372, hgb 0.367 versus logistic 0.249) and the spatial deep-learning members are best of all (resunet 0.421, cnn 0.409), exploiting the strong, spatially structured forecast-temperature signal near initialization. At extended range, in contrast, plain logistic regression on the large-scale predictors is the single best family (W3–W6 BSS 0.076, 0.066, 0.058, 0.053) because the boosting trees overfit the weak teleconnection signal. Error-diversity diagnostics explain why the stack helps unevenly: the residual U-Net is the most decorrelated member at short lead (mean absolute error-correlation 0.87 versus ~0.95 for the trees) and the only member with a clearly positive leave-one-out contribution (W1 ΔBSS +0.020), whereas at long lead the tabular members are near-redundant (error-correlation ~0.98) and an L1-sparse stack collapses onto the single logistic member. The two deeper networks, the feature-tokenizer Transformer and the deep MLP, are BSS-negative at long lead and are excluded from production. Finally, we test directly whether resolving real large-scale spatial patterns adds skill. The two-branch full-field network (§3e), which consumes the actual Z500, T850, and upper-level winds as spatial fields over the extended domain rather than as broadcast box-mean indices, is scored out-of-sample against the box-index residual U-Net baseline. It is statistically indistinguishable from, and marginally below, the baseline at every lead (ΔBSS ≈ 0 to −0.03, ΔSEDI ≈ 0 to −0.04; Table 11). This clean negative result indicates that, for this predictand, the discriminative signal is carried by the *local grid-cell forecast value*, not by the large-scale spatial *morphology* of the extended domain; the box-index simplification therefore costs no skill, and the full-field network is not merged into the production stack. **Table 11.** True-spatial full-field network vs box-index residual U-Net (out-of-sample), by lead. | block | full-field BSS | box-index BSS | Δ BSS | full-field SEDI | box-index SEDI | Δ SEDI | |-------|----------------|---------------|-------|-----------------|----------------|--------| | W1 | 0.392 | 0.421 | −0.028 | 0.723 | 0.740 | −0.017 | | W2 | 0.139 | 0.139 | −0.000 | 0.456 | 0.477 | −0.021 | | W3 | 0.079 | 0.078 | +0.001 | 0.356 | 0.342 | +0.013 | | W4 | 0.050 | 0.052 | −0.001 | 0.290 | 0.303 | −0.013 | | W5 | 0.042 | 0.046 | −0.005 | 0.263 | 0.307 | −0.043 | | W6 | 0.044 | 0.052 | −0.007 | 0.264 | 0.290 | −0.026 | --- ## 5. Discussion **Two horizons: forecast-driven and predictor-driven.** Using only the forecast temperature, a single-member S2S reforecast, after proper calibration, yields *significant* probabilistic skill for extreme-heat days to about week 3 on the national average. Four independent lines of evidence (the primary BSS, the binary discrimination scores, the secondary heat-wave target, and the model zoo) agree that this forecast-only horizon closes near W3–W4, which indicates a *predictability* limit rather than a calibrator shortfall. The feature-selection analysis then shows *why* and *how far it can be pushed*: past week 3 the forecast's own temperature signal is exhausted, and the residual predictability lives in the slowly varying boundary and large-scale states, with ENSO foremost, then the MJO, the 500-hPa flow, and antecedent persistence. Supplying those predictors and calibrating them out of fold extends significant skill to week 6, so the operational horizon is not fixed but depends on whether the slow predictors are used. This lead-dependent transition mirrors the sources of subseasonal heatwave skill documented for the Yangtze River Basin, India, West Africa, and Europe, and it places the Indochina/Vietnam case, where the ENSO signal is strong, within that broader literature. **Where and when the product is most trustworthy.** Consistent with the physics, skill persists further where predictability is intrinsically higher, in the tropical south and the pre-monsoon MAM season, which is operationally the most useful place for it to persist, since MAM is when heat is most consequential in the populous north-central regions. The strong ENSO control on Indochina heat is the plausible physical basis for the extended-range skill, and it connects the statistical result to the known synoptic and teleconnection controls on heat over the peninsula. **The land-surface null.** Unlike several land-atmosphere-coupling studies that find soil-moisture memory a key subseasonal heat predictor, our selected soil-moisture and soil-temperature candidates never clear the out-of-sample importance floor. We interpret this cautiously. The predictors here are *control-reforecast* soil states (themselves model products with their own drift), not observed or data-assimilated soil moisture, and the ENSO and persistence predictors may already absorb much of the low-frequency variance. This appears to be a limitation of the available data rather than evidence against land-surface coupling per se. **Calibration is essential; ensembles help modestly; spatial morphology does not.** The raw block-exceedance probability is badly miscalibrated, and every trainable family beats it by a large margin, so calibration, not model choice, is what makes the product usable. Combining models adds a further, small increment at short lead via stacking, but at long lead the single CNN is as good as any combination. Resolving true large-scale spatial fields, rather than broadcasting box-mean indices, does not help at all (§4k): the discriminative content is local. Each of these is a *negative* or *bounded* result that we report rather than obscure, because each is directly useful to a duty forecaster deciding where to spend effort. **Aggregation buys skill.** Moving from a per-cell binary decision to a region×week event forecast recovers substantial skill by averaging out cell-scale noise, and it is the recommended operational presentation for extended leads. **Limitations.** (1) The reforecast is control-only, so probabilities come from single-member MOS, not ensemble spread, and no real-time ensemble uncertainty is exploited in training. (2) The reforecast is delivered pre-interpolated to 0.25°, so an upstream interpolation error is baked in and no downscaling is performed. (3) A systematic gridded-versus-station warm offset means point-scale EHD risk exceeds the gridded product by 2–5 times in frequency, so the product is a gridded, not a point, statement. (4) The long-lead binary flag over-forecasts (frequency bias > 1) and must be read together with the probability and the scores. (5) The heat-wave product is verified but not yet embedded in the operational NetCDF. (6) Real-time deployment must feed a single control member built identically to the reforecast; if a real-time ensemble is issued it must be reduced to the control/mean, and any model-version mismatch documented. --- ## 6. Conclusions In this study, we built and leakage-audited a calibrated S2S probability forecast of extreme-heat days over Vietnam from a control-only reforecast. The main findings are summarized as follows. 1. Using only the forecast temperature, the calibrated product has statistically significant Brier skill over climatology to week 3, and longer in the tropical south (to W3–W4) and in the MAM season (to W6), and is far better calibrated than the raw model (BSS versus raw +0.28 to +0.52). 2. Among the predictand formulations, direct event classification and severity regression perform comparably, and the easy-negative cascade is skill-neutral; the severity variable is retained for interpretability and ensemble diversity rather than for a headline gain. 3. A leakage-safe feature-selection analysis identifies the sources of extended-range predictability, namely ENSO foremost at weeks 3–6, then the MJO, the 500-hPa circulation, the 850-hPa temperature, and antecedent persistence, with land-surface memory contributing little in this control reforecast, and it shows that supplying these predictors with out-of-fold calibration and stacking extends significant skill from week 3 to week 6. 4. A stacked ensemble is best at short lead, a deep network is the strongest single model at extended range, a transparent vote ensemble roughly matches the best single member, and a region×week aggregation is the most skillful and most operationally robust presentation. 5. Resolving true large-scale spatial fields does not improve on broadcast box-mean indices, a clean negative result that locates the discriminative signal at the grid-cell scale. The system is delivered as a CF-compliant gridded product with an independently re-derived headline and a fully reproducible, leakage-audited pipeline; its horizon, short from the forecast alone and longer once the slow predictors are used, is reported as the scientific result rather than obscured. We believe that these results are an important step toward an operational S2S heat early-warning capability for Vietnam. Further research should be directed toward a genuine reforecast ensemble (to recover spread-based uncertainty), toward observed land-surface and humidity predictors (to revisit the land-surface null and to move toward a heat-*stress* rather than a heat-*temperature* target), and toward a formal downscaling to reconcile the gridded product with point-scale exposure, which would be an interesting topic to address in a future study. --- ## Data and code availability The forecast and observation archives are read-only institutional data (ECMWF S2S reforecast; gridded and station observed `Tx` for Vietnam). All code, configuration, intermediate out-of-sample predictions, score tables, and figures are in the project repository; each quantitative claim in this manuscript maps to a named file and producing script in the Evidence pack (`docs/RESULTS_REPORT.md` §13). Key artefacts: LOYO out-of-sample predictions (`outputs/loyo_pred/`, `outputs/loyo_pred_v2*/`, `outputs/loyo_pred_v3*/`, `outputs/loyo_pred_v4ff/`, `outputs/loyo_pred_secondary/`); headline and stratified scores (`results/tables/cv_scores_by_block.csv`, `cv_bootstrap_ci.csv`, `cv_scores_by_region_block.csv`, `cv_binary_scores_by_block.csv`, `cv_scores_by_season_block.csv`, `cv_scores_by_severity_block.csv`); predictor selection (`feature_importance_by_block.csv`, `feature_selection_by_block.csv`); Stage-2/3 stack and zoo scores (`v3cal_binary_by_block.csv`, `v3zoo_scores_by_model_block.csv`, `v3_stack_scores_by_block.csv`); full-field test (`v4ff_scores_by_model_block.csv`, `v4_vs_v3_resunet.csv`); model-zoo and ensemble scores (`v2_scores_with_ensemble.csv`, `v2_single_member_scores_by_block.csv`); vote/regional tables (`v6_*`); secondary-target scores (`cv_scores_secondary_by_block.csv`, `cv_bootstrap_secondary_ci.csv`, `cv_binary_secondary_by_block.csv`); the operational product (`outputs/operational/`, `outputs/production_model/`); and the independent validator (`results/tables/validator_headline_rederivation.csv`). ## Author contributions, funding, competing interests [To be completed by the authors.] --- ## References *The reference list below collects real, citable works. Entries drawn from the project reference collection are given with their exact titles; the authors should complete full author lists, volumes, and page numbers to journal style before submission.* **Foundational methods and verification.** - Brier, G. W., 1950: Verification of forecasts expressed in terms of probability. *Mon. Wea. Rev.*, **78**, 1–3. - Murphy, A. H., 1973: A new vector partition of the probability score. *J. Appl. Meteor.*, **12**, 595–600. - Gneiting, T., and A. E. Raftery, 2007: Strictly proper scoring rules, prediction, and estimation. *J. Amer. Statist. Assoc.*, **102**, 359–378. - Gneiting, T., F. Balabdaoui, and A. E. Raftery, 2007: Probabilistic forecasts, calibration and sharpness. *J. Roy. Statist. Soc. B*, **69**, 243–268. - Ferro, C. A. T., and D. B. Stephenson, 2011: Extremal dependence indices: Improved verification measures for deterministic forecasts of rare binary events. *Wea. Forecasting*, **26**, 699–713. - Benjamini, Y., and Y. Hochberg, 1995: Controlling the false discovery rate. *J. Roy. Statist. Soc. B*, **57**, 289–300. - Zadrozny, B., and C. Elkan, 2002: Transforming classifier scores into accurate multiclass probability estimates. *Proc. 8th ACM SIGKDD*, 694–699. - Niculescu-Mizil, A., and R. Caruana, 2005: Predicting good probabilities with supervised learning. *Proc. 22nd ICML*, 625–632. - Wilks, D. S., 2011: *Statistical Methods in the Atmospheric Sciences.* 3rd ed. Academic Press, 676 pp. - Vannitsem, S., and Coauthors, 2021: Statistical postprocessing for weather forecasts: Review, challenges, and avenues in a big-data world. *Bull. Amer. Meteor. Soc.*, **102**, E681–E699. **Machine-learning and statistical-postprocessing methods.** - Breiman, L., 2001: Random forests. *Mach. Learn.*, **45**, 5–32. - Chen, T., and C. Guestrin, 2016: XGBoost: A scalable tree boosting system. *Proc. 22nd ACM SIGKDD*, 785–794. - Ke, G., and Coauthors, 2017: LightGBM: A highly efficient gradient boosting decision tree. *Adv. Neural Inf. Process. Syst.*, **30**. - Ronneberger, O., P. Fischer, and T. Brox, 2015: U-Net: Convolutional networks for biomedical image segmentation. *MICCAI*, 234–241. - Vaswani, A., and Coauthors, 2017: Attention is all you need. *Adv. Neural Inf. Process. Syst.*, **30**. - Rasp, S., and S. Lerch, 2018: Neural networks for postprocessing ensemble weather forecasts. *Mon. Wea. Rev.*, **146**, 3885–3900. - Grönquist, P., and Coauthors, 2021: Deep learning for post-processing ensemble weather forecasts. *Phil. Trans. Roy. Soc. A*, **379**, 20200092. - Schulz, B., and S. Lerch, 2022: Machine learning methods for postprocessing ensemble forecasts of wind gusts: A systematic comparison. *Mon. Wea. Rev.*, **150**, 235–257. **S2S prediction and its sources of skill.** - Vitart, F., and Coauthors, 2017: The Subseasonal to Seasonal (S2S) Prediction Project database. *Bull. Amer. Meteor. Soc.*, **98**, 163–173. - Vitart, F., 2014: Evolution of ECMWF sub-seasonal forecast skill scores. *Quart. J. Roy. Meteor. Soc.*, **140**, 1889–1899. - White, C. J., and Coauthors, 2017: Potential applications of subseasonal-to-seasonal (S2S) predictions. *Meteor. Appl.*, **24**, 315–325. - Robertson, A. W., and Coauthors, 2020: Current and emerging developments in subseasonal to decadal prediction. *Bull. Amer. Meteor. Soc.*, **101**, E869–E896. - Madden, R. A., and P. R. Julian, 1971: Detection of a 40–50 day oscillation in the zonal wind in the tropical Pacific. *J. Atmos. Sci.*, **28**, 702–708. - Lim, Y., S.-W. Son, and D. Kim, 2018: MJO prediction skill of the subseasonal-to-seasonal prediction models. *J. Climate*, **31**, 4075–4094. - Koster, R. D., and Coauthors, 2010: Contribution of land surface initialization to subseasonal forecast skill: First results from a multi-model experiment. *Geophys. Res. Lett.*, **37**, L02402. - Miralles, D. G., A. J. Teuling, C. C. van Heerwaarden, and J. Vilà-Guerau de Arellano, 2014: Mega-heatwave temperatures due to combined soil desiccation and atmospheric heat accumulation. *Nat. Geosci.*, **7**, 345–349. **Heat, heatwaves, and their prediction — Vietnam, Southeast Asia, and other regions.** *(Exact titles from the project reference collection; authors/venues to be completed to journal style.)* - Heatwaves in Southeast Asia and their changes in a warmer world. *Earth's Future*, 2021. - Li, X.-X., 2020: Heat wave trends in Southeast Asia during 1979–2018: The impact of humidity. - The characterization, mechanism, predictability, and impacts of the unprecedented 2023 Southeast Asia heatwave. - Impacts of El Niño–Southern Oscillation on heat waves in the Indochina Peninsula. - Synoptic characteristics, atmospheric controls, and long-term changes of heat waves over the Indochina Peninsula. - Spatial patterns of health vulnerability to heatwaves in Vietnam. - Spatial distribution and trends of heat stress in Vietnam. - Heat index trends and climate change implications for occupational heat exposure in Da Nang, Vietnam. - Valuation of a heatwave early warning system for mitigating risks associated with heat-related illness in Central Vietnam. - High-resolution dynamically downscaled projections of future extreme temperatures, heatwaves and exposure in Southeast Asia. - Xie, J., and Coauthors, 2020: Sources of subseasonal prediction skill for heatwaves over the Yangtze River Basin revealed from three S2S models. - Mandal, R., and Coauthors, 2025: Heat waves in India: Patterns, associations, and subseasonal prediction skill. - Pan, X., and Coauthors, 2026: Sub-seasonal prediction of extreme high temperature over the Yangtze River Basin based on deep learning. - Subseasonal prediction of central European summer heatwaves with linear and random forest machine learning models. - Pattern-based forecasting enhances the prediction skill of European heatwaves into the sub-seasonal range. - Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data. - Advancing subseasonal surface air temperature and heat wave prediction skill in China by incorporating scale interaction in a deep learning model. - Seasonal heatwave forecasting with explainable machine learning and remote sensing data. - Skillful subseasonal ensemble predictions of heat wave onsets through better representation of land surface uncertainties. - Precipitation and soil moisture coupling constrains subseasonal predictability of a prolonged extreme heatwave. **Heatwave definitions, indices, and impacts.** - Perkins, S. E., and L. V. Alexander, 2013: On the measurement of heat waves. *J. Climate*, **26**, 4500–4517. - Perkins, S. E., 2015: A review on the scientific understanding of heatwaves—Their measurement, driving mechanisms, and changes at the global scale. *Atmos. Res.*, **164–165**, 242–267. - Perkins-Kirkpatrick, S. E., and S. C. Lewis, 2020: Increasing trends in regional heatwaves. *Nat. Commun.*, **11**, 3357. - Nairn, J. R., and R. J. B. Fawcett, 2015: The excess heat factor: A metric for heatwave intensity and its use in classifying heatwave severity. *Int. J. Environ. Res. Public Health*, **12**, 227–253. - Quantifying excess deaths related to heatwaves under climate change scenarios: A multicountry time series modelling study. **Regional climate of Vietnam** *(to be completed by the authors).* - Phan-Van, T., and Coauthors, 2009: Regional climate modelling over Vietnam. - Ngo-Duc, T., and Coauthors, 2017: Trends in extreme temperature and precipitation indices over Vietnam / Southeast Asia.