the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 23 Feb 2026
Physically Coherent Machine Learning for Tropical Cyclone Storm Surge Emulation
Abstract. Climate change is projected to impact tropical cyclone magnitude and frequency, with high magnitude events becoming more common. The destructive nature of event derived storm surges and associated coastal flooding necessitates risk management. However, the historic record is too short and too sparse to assess risk effectively, resulting in incomplete probability distributions of surge heights, particularly for distribution tails. Hydrodynamic simulation can fill these gaps, but the number of simulations required, both spatially and under diverse climates, coupled with their high computational cost, is prohibitive. To address this, we present an Artificial Neural Network storm surge emulator, deployed in the northwest Gulf of Mexico. This is trained on a database of hydrodynamic simulations, and outputs spatially coherent time series of surge. Our model achieves an R2 of 0.91, with a RMSE of 13 cm when compared to an independent test set of hydrodynamic simulations, while exhibiting a computational gain factor of over 1500. Our approach is novel in its use of feature engineering to improve performance. Here variables which are physically relevant to surge are derived from commonly used features, such as wind and pressure, allowing us to maintain a simple model architecture, while steering the model towards physically coherent learning. Shapley Values are utilised for model interpretation and demonstrate that the model is making physically justified inference. Success is demonstrated by comparing our feature engineered model to a control, which engages in minimal feature engineering. The control achieves an R2 of 0.71 and a RMSE of 23 cm only.
- Preprint
(3007 KB) - Metadata XML
-
Supplement
(30732 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-781', Anonymous Referee #1, 24 Jun 2026
-
AC1: 'Reply on RC1', Hamish Wilkinson, 02 Aug 2026
Thank you for taking the time to review our paper, and for your helpful suggestions to improve it. For clarity, we will first write your comment, and then our reply.
- “The claim of broad generalisability should be moderated. Although the model performs well for the northwest Gulf of Mexico, the validation is still based on synthetic storms within a single geographical region. Additional validation using held-out storm-track families, historical events, or another coastal domain would strengthen the conclusions.”
We agree that the claims of generalisability are overly broad. To address this, we will moderate claims of generalisability to focus on within domain generalisation. In essence the headline point becomes the removal of the requirement of coordinate specific features to achieve good skill. We will reframe the rest of the paper with this context and make it clear that domain scale generalisability is the focus of subsequent work. - “The use of the term “physically coherent” needs more careful justification. SHAP analysis provides useful interpretability, but it does not by itself prove physical consistency or causality. This is particularly important because some feature responses, such as the depth-related behaviour, appear physically counterintuitive.”
We agree that SHAP does not in of itself ensure physical coherence. To address this we shall change the title to “Towards” physically coherent learning and adjust framing to highlight that SHAP values are valuable in that they flag anomalous behaviour, allowing for researcher led correction to prevent physically unjustifiable features, as well as providing confidence and guidance in feature selection. - “The model performance at extreme surge levels requires further assessment. Since storm-surge risk is mainly governed by high-end events, the authors should include additional metrics such as peak-surge error, high-quantile error, timing error of maximum surge, and nearshore-only performance.”
We agree that peak performance requires a more detailed analysis considering its importance. To do so we will:- Subset nodes to within the 95th percentile and within the 2m depth contours.
- Calculate error statistics including bias, MAE, RMSE, and Scatter Index (normalised RMSE) of peak surge.
- Calculate phase error for these same peaks with metrics on phase bias, phase MAE and phase RMSE, with these normalised to timestep.
- Provide Kling Gupta Efficiency scores on these peaks to further evaluate model performance for peak surge.
- Provide a figure (likely in SI) of a frequency distribution of surge heights for all models.
- “The computational speed-up should be clarified. The definition of the gain factor should be explicitly stated, and the reported acceleration should include preprocessing and feature-generation costs, not only inference time, to provide a realistic estimate of practical deployment efficiency.”
We agree that the computational speed up section lacks clarity, and as it is a clear driver of the use case for ML models requires further attention. To do so we shall:- Rerun the feature engineering scripts to provide timings for this stage of the process.
- Scale all timings (including ADCIRC runtimes) to the hardware used for this study to allow a like for like comparison.
- Clarify gain factor. For example, a gain factor of one hundred indicates that the ML model could run 100 simulations in the time ADCIRC runs one.
- Clarify the methods used to produce average times.
- The references and editorial presentation require careful checking. Some methodological claims, particularly regarding CNN limitations on unstructured meshes, SHAP-based physical interpretation, and computational comparison with ADCIRC, would benefit from more precise citations or more cautious wording. The manuscript also contains minor typographical issues such as “ADCRIC,” “final output later,” “ReLu,” and “Army Corp of Engineers,” which should be corrected.
As discussed in responses through comments 1 to 4 we shall address concerns regarding SHAP and computational comparisons. In regard to CNNs we will broaden our explanation to acknowledge that in this case we are referring to computer vision methods which required structured images (2D) as opposed to 1D convolutions. Finally, we shall reproof read and correct any typographical errors.
Thank you once again for your time and your expertise, it has been most helpful and will leave the paper much improved.
Citation: https://doi.org/10.5194/egusphere-2026-781-AC1 - “The claim of broad generalisability should be moderated. Although the model performs well for the northwest Gulf of Mexico, the validation is still based on synthetic storms within a single geographical region. Additional validation using held-out storm-track families, historical events, or another coastal domain would strengthen the conclusions.”
-
AC1: 'Reply on RC1', Hamish Wilkinson, 02 Aug 2026
-
RC2: 'Comment on egusphere-2026-781', Anonymous Referee #2, 29 Jun 2026
The paper presents a point-wise feedforward ANN that emulates ADCIRC storm surge over the NW Gulf of Mexico. SHAP is used for interpretation and feature selection, and a feature-engineered Fullset is benchmarked against a Naive control with about 688 to 1535 claimed speed-up. The manuscript is a clean, well-written study, and the coordinate-ablation experiment is a nice piece of design. My concerns are mostly about a gap between the framing and what the experiments actually establish. I'd see this as a solid contribution after a substantial revision of claims, not necessarily of method.
Major comments
1. The test set is not independent in the way the headline metrics imply. The 40 storms are split using a ration of 7/2/1, leaving only about 4 test storms. Note that the arithmetic doesn't match with 7/2/1 of 40, since 28/8/4 produces 36 train+val, but Table 2’s caption says 37. More importantly, all 40 are perturbations of just 5 primary synthetic tracks. So a held-out storm is almost certainly a near-neighbour in feature space of a training storm. R^2 computed over millions of spatially/temporally autocorrelated nodes from 3~4 storm realisations has a tiny effective sample size, and no uncertainty is attached to any headline number. At minimum I'd want leave-one-track-out cross-validation and error bars on R^2/RMSE before the 0.91-vs-0.71 contrast is presented as definitive.
2. Generalization is the central motivation but is never tested. The abstract, intro and conclusion are built around global deployability and extrapolation beyond the training domain, but every evaluation is in-domain on the same synthetic family. Removing coordinates and retaining performance shows the model isn't memorizing location. The paper concedes “extrapolating beyond the training domain is the next step,” which I agree with, but that concession is in tension with the title and framing. I guess adding a real out-of-domain test (another region) would be necessary.
3. “Physically coherent” is partly undercut by your own SHAP results. Depth shows an inverse-to-theory relationship in the Fullset model, which you attribute to unspecified “feature interaction” and frame as “encouraging.” That reading is charitable. Given that SHAP's additivity assumes feature independence, and you've engineered many features from shared roots, the SHAP-based physical-coherence argument is the shakiest load-bearing claim in the paper. I'd recommend softening “physically coherent” to something like “broadly consistent with surge physics,” and acknowledging the correlated-feature caveat at the point where SHAP physics conclusions are drawn, not only in Section 2.6.2.
4. The tool is weakest exactly where the stated application needs it most. The entire motivation is filling distribution tails for return periods. However the models under-predict high extremes and over-predict low extremes, and the Fullset over-predicts negative surge in bays. RMSE of 13 cm is dominated by the enormous mass of near-zero/moderate points. Tail skill is what matters for probabilistic approach (e.g., PTHA-style) use and is where performance is worst. This trade-off deserves to be foregrounded in the abstract/conclusion rather than emerging gradually through Figs 4 to 8.
5. The speed-up accounting is one-sided. 688 to 1535 times speed-up compares single-GPU inference to 1024-core ADCIRC wall-clock, but excludes the cost of the 446 ADCIRC runs needed to create the training set, and the admitted hidden per-storm feature-generation cost, and preprocessing, which you show scales with feature-set size. For the stated use case (10^5+ runs) amortization is favorable, but the headline number should be stated as inference-only and the training-data cost acknowledged.
Minor comments
1. Eq. 4 is missing a closing parenthesis.
2. Many typos were found. “ADCRIC” for ADCIRC (several places), “Fullest” for Fullset, “final output later” for “layer”, and “Long-Term-Short-Term Networks” should be “Long Short-Term Memory.”
3. The 1013 mb inclusion threshold is essentially standard SLP, so “minimum pressure never falls below 1013 mb” is a very weak filter. It is worth a sentence justifying it.
4. Table 1 leads with R^2=0.93 (with XY) while the abstract uses 0.91 (no XY). Fine given your no-coordinate focus, but state the convention once up front.
5. Tides and waves are correctly flagged as omitted, but given the “risk management” framing and your own citations on nonlinear tide-surge interaction, I'd soften “risk management” to “a step toward risk-relevant emulation.”
6. “Novel in its use of feature engineering” overstates. Coastal geometry/bathymetric predictors appear in prior surge-ML work. The contribution is better cast as a systematic feature-engineering blueprint with interpretability, which is defensible.
Citation: https://doi.org/10.5194/egusphere-2026-781-RC2 -
AC2: 'Reply on RC2', Hamish Wilkinson, 02 Aug 2026
Thank you for providing your time and your expertise to reviewing and improving our paper. It is very much appreciated. For clarity, we will first write your comment, and then our reply. In some cases, your comments align with reviewer one, and in such cases, we will refer to our replies there.
- “The test set is not independent in the way the headline metrics imply. The 40 storms are split using a ration of 7/2/1, leaving only about 4 test storms. Note that the arithmetic doesn't match with 7/2/1 of 40, since 28/8/4 produces 36 train+val, but Table 2’s caption says 37. More importantly, all 40 are perturbations of just 5 primary synthetic tracks. So a held-out storm is almost certainly a near-neighbour in feature space of a training storm. R^2 computed over millions of spatially/temporally autocorrelated nodes from 3~4 storm realisations has a tiny effective sample size, and no uncertainty is attached to any headline number. At minimum I'd want leave-one-track-out cross-validation and error bars on R^2/RMSE before the 0.91-vs-0.71 contrast is presented as definitive.”
We will split this comment out and answer them directly:
“The test set is not independent in the way the headline metrics imply. The 40 storms are split using a ration of 7/2/1, leaving only about 4 test storms.”
We agree that this is a small number of storms. We currently do sensitivity testing in SI against the 345 storms which were not used in the 7/2/1 split. However, we agree that this is inadequate and will expand on this and move it into the main text.
“Note that the arithmetic doesn't match with 7/2/1 of 40, since 28/8/4 produces 36 train+val, but Table 2’s caption says 37”
This is a typo, we will correct.
“More importantly, all 40 are perturbations of just 5 primary synthetic tracks. So a held-out storm is almost certainly a near-neighbour in feature space of a training storm.”
We agree on the limitations of the synthetic set, and we do directly address this in the text. We also take steps to reduce over fitting by only using 40 storms out of 385 candidates for the ML pipeline. Nevertheless, we agree that this cannot be overstated, so we will adjust the text to foreground this limitation. However, while the tracks do all come from 5 primary it should be noted that the issue here is primarily the shape of the track, as shown in Figure 1. By perturbing the track set characteristics used to generate wind fields, such as the max wind speed, the Holland B parameter, and indeed the radius of maximum winds, the wind field, and therefore surge, produced for each track (even when not perturbed in space) will be quite different from the original. When perturbed in space, and so with different landfall locations, this will generate unique surge profiles. Note that the original purpose of the synthetic set was for FEMA to assess flood insurance risk categories.
"R^2 computed over millions of spatially/temporally autocorrelated nodes from 3~4 storm realisations has a tiny effective sample size, and no uncertainty is attached to any headline number. At minimum I'd want leave-one-track-out cross-validation and error bars on R^2/RMSE before the 0.91-vs-0.71 contrast is presented as definitive.”
We agree that this is a flaw in the methodology. To address this, we will expand on sensitivity testing by taking the 345 storms which were not used in the ML pipeline and compute error metrics for each of these storms. We will use these to calculate standard deviations and 95% confidence intervals and report these in the text. - “Generalization is the central motivation but is never tested. The abstract, intro and conclusion are built around global deployability and extrapolation beyond the training domain, but every evaluation is in-domain on the same synthetic family. Removing coordinates and retaining performance shows the model isn't memorizing location. The paper concedes “extrapolating beyond the training domain is the next step,” which I agree with, but that concession is in tension with the title and framing. I guess adding a real out-of-domain test (another region) would be necessary.”
This closely follows reviewer 1s first comment, and I believe that our response there address this comment too. For clarity we shall repeat it here:
We agree that the claims of generalisability are overly broad. To address this, we will moderate claims of generalisability to focus on within domain generalisation. In essence the headline point becomes the removal of the requirement of coordinate specific features to achieve good skill. We will reframe the rest of the paper with this context and make it clear that domain scale generalisability is the focus of subsequent work. - “Physically coherent” is partly undercut by your own SHAP results. Depth shows an inverse-to-theory relationship in the Fullset model, which you attribute to unspecified “feature interaction” and frame as “encouraging.” That reading is charitable. Given that SHAP's additivity assumes feature independence, and you've engineered many features from shared roots, the SHAP-based physical-coherence argument is the shakiest load-bearing claim in the paper. I'd recommend softening “physically coherent” to something like “broadly consistent with surge physics,” and acknowledging the correlated-feature caveat at the point where SHAP physics conclusions are drawn, not only in Section 2.6.2.”
This makes the same point as the first reviewers second point. And I believe that our response there also broadly addresses this comment. To address some points here directly:- We agree that depths inversion should not be attributed to feature engineering where there is no clear proof of this. We will update the text and frame this as a clear limitation of ML models. Even with interpretability methods such as SHAP there remains an element of black box.
- We agree that SHAP assumes independence, and as such that means that it is almost certain that features will be sharing importance and indeed interacting. We also agree that this limits analysis of physical coherence. We will reframe the text in the manner recommended, as discussed in our response to reviewer one.
- We agree that depths inversion should not be attributed to feature engineering where there is no clear proof of this. We will update the text and frame this as a clear limitation of ML models. Even with interpretability methods such as SHAP there remains an element of black box.
- “The tool is weakest exactly where the stated application needs it most. The entire motivation is filling distribution tails for return periods. However the models under-predict high extremes and over-predict low extremes, and the Fullset over-predicts negative surge in bays. RMSE of 13 cm is dominated by the enormous mass of near-zero/moderate points. Tail skill is what matters for probabilistic approach (e.g., PTHA-style) use and is where performance is worst. This trade-off deserves to be foregrounded in the abstract/conclusion rather than emerging gradually through Figs 4 to 8.”
This is similar to reviewer one’s views in comment 3. Our response is repeated here for clarity:
We agree that peak performance requires a more detailed analysis considering its importance. To do so we will:
-
- Subset nodes to within the 95th percentile and within the 2m depth contours.
- Calculate error statistics including bias, MAE, RMSE, and Scatter Index (normalised RMSE) om peak surge.
- Calculate phase error for these same peaks with metrics on phase bias, phase MAE and phase RMSE, with these normalised to timestep.
- Provide Kling Gupta Efficiency scores on these peaks to further evaluate model performance for peak surge.
- Provide a figure (likely in SI) of a frequency distribution of surge heights for all models.
In direct response to reviewer two, we will foreground the results of this analysis in both the abstract and conclusion, alongside the updated error metrics from sensitivity testing.
- “The speed-up accounting is one-sided. 688 to 1535 times speed-up compares single-GPU inference to 1024-core ADCIRC wall-clock, but excludes the cost of the 446 ADCIRC runs needed to create the training set, and the admitted hidden per-storm feature-generation cost, and preprocessing, which you show scales with feature-set size. For the stated use case (10^5+ runs) amortization is favorable, but the headline number should be stated as inference-only and the training-data cost acknowledged.”
This aligns with reviewer one’s 4th comment, and we agree it is justified. For clarity we provide our response to reviewer one here:
We agree that the computational speed up section lacks clarity, and as it is a clear driver of the use case for ML models requires further attention. To do so we shall:
-
- Rerun the feature engineering scripts to provide timings for this stage of the process.
- Scale all timings (including ADCIRC runtimes) to the hardware used for this study to allow a like for like comparison.
- Clarify gain factor. For example, a gain factor of one hundred indicates that the ML model could run 100 simulations in the time ADCIRC runs one.
- Clarify the methods used to produce average times.
“Minor comments"
- Eq. 4 is missing a closing parenthesis.
- Many typos were found. “ADCRIC” for ADCIRC (several places), “Fullest” for Fullset, “final output later” for “layer”, and “Long-Term-Short-Term Networks” should be “Long Short-Term Memory.”
- The 1013 mb inclusion threshold is essentially standard SLP, so “minimum pressure never falls below 1013 mb” is a very weak filter. It is worth a sentence justifying it.
- Table 1 leads with R^2=0.93 (with XY) while the abstract uses 0.91 (no XY). Fine given your no-coordinate focus, but state the convention once up front.
- Tides and waves are correctly flagged as omitted, but given the “risk management” framing and your own citations on nonlinear tide-surge interaction, I'd soften “risk management” to “a step toward risk-relevant emulation.”
- “Novel in its use of feature engineering” overstates. Coastal geometry/bathymetric predictors appear in prior surge-ML work. The contribution is better cast as a systematic feature-engineering blueprint with interpretability, which is defensible.”
Authors response:
- We will add the missing parenthesis.
- Re-proofread and correct typos.
- We will remove this – it refers to subsetting the storm to arrival in the ML domain, not a filter on storm selection. This is better explained later in the text. Agree it is confusing.
- We will make the case for giving primacy for coordinate removed error metrics clear.
- Agreed, will update the text to reflect.
- Agreed, will update the text to reflect.
Citation: https://doi.org/10.5194/egusphere-2026-781-AC2 - “The test set is not independent in the way the headline metrics imply. The 40 storms are split using a ration of 7/2/1, leaving only about 4 test storms. Note that the arithmetic doesn't match with 7/2/1 of 40, since 28/8/4 produces 36 train+val, but Table 2’s caption says 37. More importantly, all 40 are perturbations of just 5 primary synthetic tracks. So a held-out storm is almost certainly a near-neighbour in feature space of a training storm. R^2 computed over millions of spatially/temporally autocorrelated nodes from 3~4 storm realisations has a tiny effective sample size, and no uncertainty is attached to any headline number. At minimum I'd want leave-one-track-out cross-validation and error bars on R^2/RMSE before the 0.91-vs-0.71 contrast is presented as definitive.”
-
AC2: 'Reply on RC2', Hamish Wilkinson, 02 Aug 2026
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,949 | 922 | 80 | 2,951 | 220 | 94 | 161 |
- HTML: 1,949
- PDF: 922
- XML: 80
- Total: 2,951
- Supplement: 220
- BibTeX: 94
- EndNote: 161
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The manuscript presents a timely and practically relevant study on storm-surge emulation using physically motivated feature engineering within an artificial neural network framework. The work addresses an important computational bottleneck in hydrodynamic storm-surge modelling and demonstrates clear improvement over a naive machine-learning baseline. The manuscript is generally well motivated and has publication potential; however, the following points should be addressed before publication.