Coastal areas are increasingly vulnerable to compound 昀氀ood events, which occur when storm surges, river over昀氀ow, and heavy rainfall occur simultaneously or in close succession (Eilander et al., 2023;MuÞ noz et al., 2021;Xu et al., 2024). Although these events have historically been rare, their frequency and intensity have risen signi昀椀cantly in recent years (Taherkhani et al., 2020;Wing et al., 2022). This increase is primarily attributed to climate change, which drives sea level rise, marine heatwaves, torrential rainfall, intense winds from tropical cyclones, glacier melting, land subsidence, and altered ocean currents (Bloom昀椀eld et al., 2023;Levy et al., 2024;Ohenhen et al., 2023;Radfar et al., 2024;Santiago-Collazo et al., 2019;Thiéblemont et al., 2024;Volkov et al., 2023). Currently, approximately 700 million people and an estimated $13 trillion in assets located in areas less than 10 m above mean sea level worldwide are at risk of 昀氀ooding (Kirezci et al., 2023). In the United States alone, tropical cyclones intensify rainfall and river over昀氀ow, contributing to over $13.8 billion in 昀氀ood damages and 86 昀氀ood fatalities out of a total $280 billion in cyclone-related damages and 683 deaths over the past 昀椀ve years (NOAA-NCEI, 2024). Hence, the imperative to understand 昀氀ood dynamics caused by compound events is driving the exploration of advanced modeling approaches that can enhance current prediction methods and support 昀氀ood management practices.
Methods for predicting compound and coastal 昀氀ood dynamics rely on both physically-based and data-driven approaches (Chiang et al., 2024;Fraehr et al., 2022;Hu et al., 2019;Löwe et al., 2021;Marsooli and Wang, 2020;Sampurno et al., 2022;Shahabi and Tahvildari, 2024). Physically-based approaches, such as hydrodynamic and hydraulic models, are particularly effective at estimating the interaction of multiple 昀氀ood drivers such as rainfall, river discharge, and storm surges (Bates, 2022(Bates, , 2023;;Leijnse et al., 2021;Santiago-Collazo et al., 2024). These approaches demonstrate strong predictive skills across multiple 昀氀ood scenarios by complying with physical constraints such as the conservation of mass and momentum over time and space (Alipour et al., 2022;Camus et al., 2021;Jibhakate et al., 2023;PeÞ na et al., 2022;Zhong et al., 2024). However, they require suf昀椀cient computational resources and observational (forcing) data, even for short-term (hourly to daily) predictions over large-scale domains (Bilskie et al., 2021;Nezhad et al., 2023). Data-driven models are generally ef昀椀cient in this regard since they leverage statistical and machine learning techniques to learn from nonlinear associations and infer data patterns, thus rapidly generating 昀氀ood predictions without relying on complex physical simulations (Gomez et al., 2024;Lewis et al., 2024;Moftakhari et al., 2017).
Statistical approaches typically require explicit assumptions about data distribution to predict 昀氀ood data associated with a given return period (e.g., peak storm surge, river 昀氀ow, and rainfall) (Boumis et al., 2023;Maduwantha et al., 2024;Moftakhari et al., 2021;Zhong et al., 2024).
Machine learning techniques, particularly deep learning (DL) models, can capture complex nonlinear associations and hidden patterns from input data features; thus, enhancing the prediction accuracy (Foroumandi et al., 2024;Fu et al., 2022;Sattari et al., 2025). Recently, transfer learning techniques have been incorporated into DL model architectures to further enhance the predictive capability in areas with limited or scarce data (Obara and Nakamura, 2022;Seleem et al., 2023). Previous studies have noted accurate prediction of storm surges, extreme water levels, river 昀氀ow, among other 昀氀ood drivers (Daramola et al., 2025;Green et al., 2025;Hussain and Khan, 2020;Mahakur et al., 2025;McKeon and Piecuch, 2025;MuÞ noz et al., 2021;Samantaray et al., 2025;Tang et al., 2025;Tiggeloven et al., 2021). Nevertheless, DL models are mostly used for station-based 昀氀ood predictions with hazard and risk map generation focused on the peak 昀氀ood extent (Ayyad et al., 2022;Tedesco et al., 2024). This is because many conventional DL models excel at either spatial (e.g., Graph Convolution Networks, Convolutional Neural Networks) or temporal analysis of 昀氀ood events (e.g., Long Short-Term Memory networks, Gated Recurrent Units). As a result, they often lack the ability to process 昀氀ood dynamics in both space and time. To overcome this limitation, studies tend to use DL frameworks that combine spatial and temporal architectures, preferably employing static datasets for spatial analyses and dynamic features for temporal contribution (Farahmand et al., 2023;Fathi et al., 2025). Additionally, DL models designed for spatiotemporal analysis like Convolutional Long Short-Term Memory (ConvLSTM) may effectively capture general spatiotemporal patterns but struggle to accurately represent the varying conditions across the spatial domain. The latter results from the lack of inherent physical constraints in DL models to guide 昀氀ood dynamics in contrast to those present in hydrodynamic models. Furthermore, regional variations in spatial 昀氀ood dependencies are ignored when conducting current and future 昀氀ood risk assessments (Brunner et al., 2020).
Given the challenges described, a key question arises: which DL techniques can effectively capture hysteresis in system behavior across spatiotemporal domains? In 昀氀ood dynamics, hysteresis occurs when water levels depend not only on current conditions but also on prior states, producing lagged effects that vary in space and time (Wu et al., 2023). Spatially, water level variations at one location can induce changes in water levels of adjacent areas due to factors like bottom friction (roughness) effects on the propagation of 昀氀ood waves as well as morphological features that either attenuate or amplify water levels in coastal and estuarine systems (Hoitink and Jay, 2016;Prandle, 1985;Talke and Jay, 2020). Temporally, peak water level timing might differ across inland, transition, and coastal regions, resulting in peak water levels occurring earlier or later depending on distance from the 昀氀ood source, local topography and bathymetry, or 昀氀ood connectivity. This complex interplay challenges traditional DL models, which often fail to capture these lagged dependencies effectively (Brunner et al., 2020). As a result, such models may erroneously assume a single peak water level across a large domain (Yu et al., 2024). This oversight can miss critical timing differences in peak water levels, such as those along tidally-in昀氀uenced rivers or estuaries driven by storm surges, river 昀氀ow, and backwater effects (Hoitink and Jay, 2016;Sandbach et al., 2018). To model hysteresis accurately, advanced DL frameworks must integrate spatiotemporal dependencies across the model domain including the interconnected nature of 昀氀ooding processes and different peak water level timings.
To address these limitations, we propose architectural modi昀椀cations that enable DL models to better capture the spatiotemporal dynamics of compound 昀氀ood events, as simulated in hydrodynamic models. Specifically, our approach involves three key modi昀椀cations: (i) implement mechanisms that enable the model to selectively focus on and weigh the importance of features, enhancing their ability to understand and replicate intricate spatiotemporal dependencies; (ii) cluster the spatial domain regionally and modulate with temporally varying contributions based on water level dynamics at observational stations to re昀氀ect water level variability over time and space; and (iii) leverage dynamic input data from all 昀氀ood drivers (just like in hydrodynamic models) to facilitate the learning process of nonlinear interactions and improve the model's predictive accuracy. Our proposed DL model is applied to Galveston, TX, on the Gulf Coast of the United States, a region that has been signi昀椀cantly affected by torrential rainfall and cyclone-induced 昀氀ood events.
The Galveston Bay (G-Bay) is the seventh largest estuary in the US and serves as a vital hydrological link between Houston, Texas, and the Gulf of Mexico through a complex network of bayous, interior bays, and rivers (MuÞ noz et al., 2024). Extending 56 km in length and averaging 31 km in width, G-Bay, with a shallow depth of approximately 2 m, encompasses an expansive surface area of about 1600 km 2 (Fig. 1a). The combination of substantial freshwater discharge and dynamic tidal forces creates a sensitive hydrodynamic environment susceptible to rapid changes in water levels. Additionally, G-Bay's complex morphology resembles a bottle-like neck that connects the Buffalo Bayou River with the Bay. Such a morphologic feature exacerbates the impacts from coastal and inland 昀氀ood drivers, including storm surges and rainfall-runoff (MuÞ noz et al., 2022;Valle-Levinson et al., 2020), and makes it particularly susceptible to compound 昀氀ood events affecting both Galveston and Harris Counties.
For this study, we selected six of the most recent 昀氀ooding events that have affected G-Bay, including Hurricanes Ike (2008), Harvey (2017), Nicholas (2021), and Beryl (2024), as well as the torrential rainfall events of Memorial Day (2015) and Tax Day (2016). Making landfall near Galveston on September 13, 2008, as a Category 2 hurricane, Hurricane Ike brought maximum sustained winds of approximately 175 km/h and a storm surge exceeding 4.5 m (National Weather Service, 2008). The hurricane caused extensive 昀氀ooding in Galveston and Harris Counties. Torrential rainfall on May 23-24, 2015, around Memorial Day resulted in widespread 昀氀ooding across South-Central Texas. Southern Blanco County received a record of 0.25-0.33 m leading to Blanco River at Wimberley rising from nearly 1.5 m to over 12.5 m within hours, causing severe 昀氀ooding (National Weather Service, 2015). Similarly, about 0.61 m of rain fell in Houston and in Waller County from April 17-18, 2016, 昀椀lling local reservoirs to full capacity and hitting record high water levels (ABC News, 2019).
Harvey made landfall near Rockport, Texas, on August 25, 2017, as a Category 4 hurricane with winds of 215 km/h. The storm stalled over Southeast Texas, dropping over 1.27 m of rain in some areas, resulting in catastrophic 昀氀ooding, particularly in the Houston metropolitan area (National Environmental SatelliteDataand Information Service, 2024). Nicholas made landfall on September 14, 2021, near Sargent Beach, Texas, as a Category 1 hurricane with maximum sustained winds of 120 km/h, causing signi昀椀cant rainfall and 昀氀ooding in Southeast Texas and impacting communities around Galveston Bay (National Hurricane Center, 2022). Beryl made landfall near Matagorda, Texas, on July 8, 2024, as a Category 1 hurricane with winds of 129 km/h, leading to extensive 昀氀ooding in Galveston and Harris Counties (NASA earth observatory, 2024). The six events resulted in over 210 fatalities and approximately $200 billion in economic damages. This high susceptibility to compound 昀氀ooding renders Galveston Bay an ideal site for applying the proposed DL model.
We developed a DL framework that incorporates information from a physically-based model, namely Delft3D-FM, to estimate compound 昀氀ood dynamics during extreme events. Both the DL and the Delft3D-FM models utilize the same input features: observed water levels and river discharge at speci昀椀c stations, along with spatially distributed variables, such as precipitation, digital elevation model (DEM), atmospheric pressure, and wind speed (Fig. 2). These input features capture key spatiotemporal patterns from variables driving compound 昀氀ood dynamics, enabling both models to predict 昀氀ood depth and inundation extent. Since the DL model aims to replicate the model predictions of Delft3D-FM, the simulated water depth data of the latter model is the target feature of the proposed framework during training. The framework employs the LSTM model to analyze water level data from multiple tide-gauge stations (Fig. 2a and b). Through an attention mechanism, which is a means to solve the problem of information overload (Niu et al., 2021), the model identi昀椀es critical time steps that re昀氀ect complex 昀氀ood dynamics across the area. The attention vectors generated by this step are then integrated into the ConvLSTM model, which captures spatially distributed 昀氀ood patterns, allowing for nuanced time-dependent interactions across the domain. Enhanced by its own attention mechanism, the ConvLSTM model learns both the unique contributions of individual spatial features and their interactions over time and space (Fig. 2c and d). The attention vectors from the stations are then combined with the ConvLSTM output, in an attempt for the framework to improve 昀氀ood predictions (Fig. 2b and d). The framework's 昀椀nal output is the spatiotemporal evolution of 昀氀ood depth and inundation extent (Fig. 2e).
We leverage a previously developed Delft3D-FM model of G-Bay that simulates compound 昀氀ood dynamics associated with hurricane events (MuÞ noz et al., 2024).
The model employs an unstructured 昀椀nite-volume mesh grid composed of triangular cells that vary in size. The grid resolution ranges from 3 km at the ocean boundary to 5 m within Harris County, enabling high-resolution simulation of the complex interactions between G-Bay's physical morphologic and hydrodynamic processes. The model is driven by multiple external forcing conditions, including tidal harmonic constituents at the open ocean boundary, derived from the TPXO 8.0 global inverse tide model; hourly riverine in昀氀ows from USGS river gauges; and hourly wind speed and atmospheric pressure data from the ERA5 reanalysis dataset with a spatial resolution of ~31 km. Since rainfall data is a key 昀氀ood driver in inland areas, we interpolated these data using a dense network of rain gauges in Harris County. Following Sebastian et al. (2021) and MuÞ noz et al. ( 2024), we set the "inverse distance weight" as the interpolation method in ArcGIS with an output cell size of 1 km (e.g., shortest Euclidean distance between existing rain gauges), a search radius of 5 points, and a power function of 2. Moreover, we supplemented rain gauge data with "total precipitation" from ERA5 to estimate rainfall patterns in coastal areas beyond Harris County and over the Gulf of Mexico.
The G-Bay model has been rigorously calibrated and validated using all the above-mentioned compound 昀氀ood events as shown in MuÞ noz et al. (2024) as well as Figs. S3-S6 in the supplementary material. Legacy DEMs, including the 2006 Galveston, Texas Coastal Digital Elevation Model and the Continuously Updated Digital Elevation Model (CUDEM), are used to simulate the conditions during past and recent hurricane and storm events. Additionally, the National Land Cover Database (NLCD) provides annual land cover maps at a 30-m resolution, which are crucial for inferring spatially distributed roughness values during the calibration process. Model calibration involves adjusting the Manning's roughness coef昀椀cient for different land cover types to ensure that simulated water levels closely align with those of observed data. To further optimize the model's performance, 50 ensemble simulations were conducted using high-performance computing resources to 昀椀nd the optimal combination of roughness coef昀椀cients that minimize errors of simulated water level and inundation depth in terms of RMSE, KGE, and NSE altogether. Observed water levels and ground-truth inundation depths were obtained from NOAA's Tide & Currents and high-water marks from the USGS's Flood Event Viewer. The ensemble members consisted of a unique combination of plausible roughness values according to the NLCD land cover class. Such combinations were obtained from the Latin Hypercube Sampling technique and considered for the compound 昀氀ood events (Helton and Davis, 2003;MuÞ noz et al., 2022). This comprehensive calibration ensures that the Delft3D-FM model accurately reproduces water levels, inundation depth, and 昀氀ood extent in G-Bay under a variety of scenarios.
The proposed DL framework integrates LSTM and ConvLSTM to process point-based and spatially distributed data, respectively. These different data types are initially processed separately to allow the combined model to exert varying localized in昀氀uences across broader spatial features. The proposed model is trained and validated on sequential and representative 昀氀ooding events including Hurricane Ike (2008), Memorial Day (2015), and Tax Day ( 2016). Both point and spatial data of these events are processed in a 6-h sequence to effectively learn patterns over meaningful and practical time intervals without adding excessive complexity in the model. In this regard, the U.S. National Hurricane Center provides forecasts every 6 h that are eventually used to run physically-based models. We use a Bayesian optimization to 昀椀rst select the best combination of hyperparameters after 50 tuning trials, while the model training for the best model is con昀椀gured for up to 300 epochs. We utilize the Adam optimizer with mean squared error as the loss function. The model validation is monitored with two critical callbacks (model checkpoint and early stopping). The model checkpoint saves the best model whenever the validation loss improves, ensuring that the 昀椀nal model retains the best weights achieved during training. Early stopping halts the training process if the validation loss does not improve after 10 consecutive epochs, preventing over昀椀tting and ensuring that the best-performing model is saved. L2 regularization is also incorporated to prevent over昀椀tting by adding a penalty to the loss function based on the magnitude of the model's weights (Xie et al., 2022).
Hourly water level data from 21 observation stations provide localized temporal dynamics of the extreme 昀氀ood events. The time-indexed water level data are loaded into structured arrays, with each array corresponding to a different observation station. For each station i, the time series data is processed through two sequential LSTM layers, with an attention mechanism that captures unique temporal patterns and computes the relative importance of each timestep at the observation stations (Chaudhari et al., 2021).
The 昀椀rst LSTM layer processes the input sequence, generating hidden states for each timestep (Equation ( 1)). These hidden states serve as input to the second LSTM layer, which produces a 昀椀nal hidden state that encapsulates the model's comprehensive understanding of the entire sequence (Equation ( 2)) (Chaudhari et al., 2021). The hidden states at each timestep for both layers are computed as follows:
where h t-1,i is the hidden state from the previous timestep, x t,i is the input (water level data) at timestep t, h t,i is the hidden state at the current timestep, and h
(2)
T,i is the 昀椀nal hidden state of the second LSTM layer for station i.
ii. Temporal attention mechanism: A unique custom attention mechanism is applied to each observation station. For each station i, the attention score e t,i at timestep t is computed using Equation (3):
where W and b are learnable parameters within the attention mechanism.
Rather than performing a simple dot product between the hidden states of LSTM layers, the attention mechanism increases the attention score of the highest attention weights to enhance the model's focus on the most critical timesteps (e.g., peak water level). The attention weights aw t,i , are calculated by applying a softmax function to the attention scores. This is then followed by amplifying the in昀氀uence of the top 10 % of the attention scores using an emphasis factor (Daramola et al., 2025) (Equation ( 4)):
Both the percentage and emphasis factor are iteratively explored until the best model performance is achieved. The process improves the model's ability to accurately capture magnitude and timing of extreme water levels (Daramola et al., 2025). The attention weights aw t,i are utilized to compute the attention vector av i in Equation (5). This vector represents a weighted sum of the LSTM hidden states across all timesteps (Chaudhari et al., 2021):
Since each hidden state is a vector, the attention vector will have the same dimensionality as the individual hidden states at all stations. Finally, the attention vectors are extracted and used to modulate the output of the spatial features processed by the ConvLSTM. It is important to note that each attention vector is not a single, 昀椀xed value applied across the entire modeling period. Instead, it dynamically represents the importance of each timestep within a sequence of data, capturing the evolving conditions of the event.
Spatial features, including atmospheric pressure, wind speed, DEM, precipitation, and water depth, are processed as 50-m resolution Geo-TIFF images at hourly timesteps. The DEM, which remains static for the duration of each event (e.g., often 1-week), is replicated to match the temporal resolution of the other features. These spatial data are loaded into structured arrays while maintaining the spatial coordinate reference system (CRS) and transformation metadata, which are critical for conducting a geospatial analysis. To address the varying scale among features, minmaxscaler normalization is applied. Subsequently, the input data for the ConvLSTM is stacked along the channel dimension, where each channel corresponds to a spatial feature (i.e., atmospheric pressure, wind speed, DEM, precipitation), with water depth designated as the target or predicted variable. A mask is used to track invalid cells containing NaN values within the spatial feature map. This mask is utilized in a custom loss function and during the spatial attention mechanism to ensure that missing data do not adversely affect the model's training (Fig. S2 in supplementary material). After the minmaxscaler normalization, the values of valid cells across all features are scaled between 0.1 and 1.0, with all NaN cells set to 0. In the denormalization process, the predicted output is similarly masked, restoring the NaN values to ensure consistency with the original data structure.
i. Convolutional and recurrent operation: At each timestep t, the ConvLSTM applies convolutional 昀椀lters to the input spatial data to extract local spatial features (Equation ( 6)), using a kernel size of 3 × 3 pixels (or 150 × 150 m). The ConvLSTM updates its hidden states over time, like those of the LSTM (Equation ( 7)), but it also processes spatial information through convolutions (Gavahi et al., 2021(Gavahi et al., , 2023)). The convolution operation involves sliding a kernel (or 昀椀lter) over the input feature maps to produce a new output (MuÞ noz et al., 2021). Speci昀椀cally, we use ConvLSTM consisting of three layers for balance of complexity and ef昀椀ciency (Table 1). Note that adding too many layers can lead to over昀椀tting or excessive computational cost.
Layer normalization is applied after the 昀椀rst two ConvLSTM layers to stabilize the training process (Equation ( 8)). Unlike batch normalization, it operates independently of batch size by normalizing across the features of a single timestep, making it more robust where memory constraints result into a smaller batch size. By treating each timestep independently, layer normalization ensures that the model effectively learns from the unique characteristics of each timestep.
where X
(l)
t is the output feature map of the l-th layer at timestep t, h
t are the unnormalized and normalized hidden state at timestep t in layer l, W (l) and b (l) are the 昀椀lter weights and biases of the l-th convolutional layer, μ and σ 2 are the mean and variance of the activations in layer l, γ and β are learnable scaling and shifting parameters, and ϵ is a small constant for numerical stability.
ii. Spatial attention mechanism: To enhance the model's focus on 昀氀ooded regions, a Convolutional Block Attention Module (CBAM) (Woo et al., 2018) is integrated after each of the 昀椀rst two ConvLSTM layers. The CBAM applies both channel (feature) and spatial attention modules. First, the channel attention module enables the model to emphasize relevant features effectively (see supplementary material). Then, the spatial attention module emphasizes relevant spatial locations by computing attention weights α t (x, y) for each spatial location at timestep t (Equation ( 9)), which is used to modulate the output feature map (Equation ( 10)). The attention-modulated feature maps from both CBAM modules are averaged to produce the 昀椀nal modulated feature map (Equation ( 11)).
where X channel att t is the channel-attended feature map at location (x,y), * denotes convolution, X
(1),att t (x, y), and X
(2),att t (x, y) are the outputs from the 昀椀rst and second CBAM modules, respectively.
The ConvLSTM model outputs these spatiotemporal feature maps that encode both the spatial patterns via convolutions as well as temporal dependencies via the recurrent updates. The output from the 昀椀nal ConvLSTM layer (X T ) is 昀氀attened into a vector (s) that contains all the spatial-temporal information learned by the ConvLSTM up to the 昀椀nal timestep T.
The study domain is partitioned into clusters using Voronoi tessellation (Fukami et al., 2021), where the number of clusters (N) corresponds to number of observation stations (i 1 , i 2 , …, i N ) (Table S1 in the supplementary material). This segmentation enables the model to effectively capture the varying 昀氀ood dynamics associated with each station and its immediate surrounding region (Fig. 3). This assumption is based on the expectation that cells surrounding the stations within each cluster are likely to experience water level peaks at nearly the same time. For each observation station i, the Voronoi cluster V i is de昀椀ned as:
where p is a point in the spatial domain, R 2 refers to the 2-dimensional Euclidean space, 6 ⋅6 is the Euclidean distance, s i is the location of station i, and s j are the locations of other stations. Binary cluster masks are created to identify the cells belonging to each Voronoi cluster (Equation ( 13)). Using a mutual exclusivity technique, we ensure that spatial cells are associated with at most one cluster, preventing overlapping in昀氀uences. These masks are crucial for applying the temporal attention vectors on the spatial feature map in a localized manner during the model training.
iv. Applying localized temporal attention through selective assignment: After extracting the temporal attention vectors for each station, they are transformed and applied to modulate the spatial feature map (Table S2 in the supplementary material). This approach selectively assigns the attention vectors to modulate their speci昀椀c corresponding clusters, which could lead to spatially varying localized in昀氀uences. Each station provides an attention vector, denoted as av i , which captures key temporal patterns from water level dynamics (Equation ( 5)). Recall that the individual attention vector is only within a cell in spatial feature map. Hence, each attention vector is passed through a fully connected (dense) layer that projects it into a 1D tensor, flattened vector i , having a length equal to the total number of cells in the spatial feature map (H ×W) (Equation ( 14)). The flattened vector i is then reshaped into a 2D array, reshaped vector i , which aligns attention information with the grid-like shape of the spatial map (Equation ( 15)).
)
The Voronoi cluster mask, mask i (x, y), corresponding to a speci昀椀c spatial region of each station is now used to con昀椀ne the in昀氀uence of the attention vector on each reshaped vector i into a localized vector i (Equation ( 16)). This ensures that the localized vector i only has attention vector within the cluster, effectively masking out areas outside the cluster. As a result, only the values at points (x, y) within each cluster V i is retained. Lastly, a combined vector(x, y) is initialized with zero-values to match the shape of the spatial feature map. For each station i, we assign its localized vector i (x, y) to the corresponding cluster locations in combined vector(x, y) using a function (Equation ( 17)). Since the Voronoi clusters are mutually exclusive, each cell in the spatial feature map is in昀氀uenced by exactly one attention vector.
{ localized vector i (x, y)
Modulating the ConvLSTM output: The combined vector(x, y) is used to modulate the ConvLSTM output s through element-wise multiplication (Equation ( 18)). This modulation emphasizes or deemphasizes the ConvLSTM output based on high or low learned attention weights (Chenmin et al., 2024) across the spatial regions, respectively.
modulated output(x, y) = s(x, y) × combined vector(x, y) (
The modulated output is 昀氀attened and passed through fully connected layers z to generate the 昀椀nal output y (Equation ( 19)), which is then reshaped back to the spatial feature map's dimensions (Equation ( 20)):
The proposed framework offers a signi昀椀cant advancement in modeling the spatiotemporal complexities of cyclone-induced 昀氀ood dynamics, surpassing the capabilities of conventional ConvLSTM models and existing hybrid approaches (Fig. 4). While conventional ConvLSTM models effectively capture general spatiotemporal patterns, they are not suitable to account for the region-speci昀椀c temporal variability and spatial heterogeneity inherent in 昀氀ood events. This limitation becomes particularly evident in large and/or diverse regions, such as Galveston Bay, where 昀氀ood propagation and peak water level timing exhibit substantial variation between coastal, transition, and inland zones. To overcome these shortcomings, our approach incorporates a novel cluster-based temporal attention mechanism. This framework partitions the spatial domain into distinct clusters, assigning localized water level dynamics to each cluster. By aligning temporal observations with region-speci昀椀c spatial predictions, the model accurately captures the varying timing of peak water levels across extensive areas. Moreover, the proposed framework leverages dynamic inputs that evolve over time to provide comprehensive spatial and temporal information. In contrast, many hybrid models rely heavily on static features, such as elevation, 昀氀oodplain, land use, etc., to ensure spatial consistency (Farahmand et al., 2023;MuÞ noz et al., 2024;Valle-Levinson et al., 2020), often at the expense of capturing temporally varying 昀氀ood behaviors. This innovation enables the model to precisely replicate the varied timing of peak water levels across extensive spatial domains, potentially improving its robustness and predictive accuracy over conventional ConvLSTM and hybrid methods.
We evaluate the trained model's performance on unseen data to assess its ability to generalize and make accurate predictions of water depth and extent during 昀氀ooding events. Testing is conducted using process-based model simulations of Hurricane Harvey (2017), Hurricane Nicholas ( 2021) and Hurricane Beryl (2024). The input and output feature processing follow the same steps as in the model training, ensuring consistency in data handling. The cluster-based training approach and the application of localized attention through selective assignment are already built into the model. Nevertheless, the accuracy of the prediction relies on using the attention model to extract attention vectors from the test data (Equations ( 21) and ( 22)). This provides insight into how the model weighs different timesteps from the test events. This cluster-based attention mechanism allows the model to focus on different spatial regions during prediction, providing more re昀椀ned, localized predictions of inundation depth and 昀氀ood extent.
vectors = attention model ( water level test sequences ) (21) y pred = best model ( X test sequences , water level test sequences )
Recommended performance metrics for predicting inundation depth are computed for each 昀氀ood event, including the Root-Mean Square Error (RMSE), coef昀椀cient of determination (R 2 ), Kling-Gupta Ef昀椀ciency (KGE) and Nash-Sutcliffe Ef昀椀ciency (NSE) (Gupta et al., 2009;Mahakur et al., 2025;Nash and Sutcliffe, 1970;Samantaray et al., 2025). In addition, hit rate (H), false alarm ratio (F), critical success index (C), and error bias (E) are used to compute the accuracy of 昀氀ood map replication (Wing et al., 2022). H measures the ratio of correctly predicted 昀氀ood instances (true positives) to the total number of actual 昀氀ood instances, indicating the model's ability to identify 昀氀ooding when they occur. F indicates the proportion of false positives among all instances where the model predicted a 昀氀ood, measuring how often the model incorrectly predicts 昀氀ooding when there is none. C assesses the proportion of correctly predicted 昀氀ood events by considering all hits, false alarms, and misses, offering a comprehensive measure of the model's accuracy in predicting 昀氀oods. E re昀氀ects the ratio of predicted 昀氀ood instances (both true positives and false positives) to actual 昀氀ood instances (both true positives and false negatives), revealing whether the model has a tendency to overpredict or underpredict 昀氀ooding. Best results are closer to 1 except for F, which should be closer to 0. We calculated the residuals by taking the direct difference between the actual and predicted 昀氀ood maps of the water depth to evaluate prediction errors.
The Bayesian optimization process systematically searches the hyperparameter space across 50 trials, successfully identifying con昀椀gurations that minimize the validation loss. This approach iteratively tunes the hyperparameter values based on prior trial outcomes. For each trial, the optimal model is derived using the Adam optimizer with mean squared error as the loss function, achieving convergence within the 300-epoch limit. The model checkpointing callback saves the best model weights based on the lowest validation loss achieved during training, and the early stopping callback terminates training after 10 consecutive epochs without improvement in validation loss. This approach resulted in models that not only achieved the lowest validation loss per trial but also demonstrated strong generalization to unseen data, effectively capturing the underlying patterns in the training set. The consistent improvement in validation loss across trials demonstrates the effectiveness of this method in identifying a robust model con昀椀guration. Fig. 5 illustrates the convergence of the optimization process, showcasing the performance of the optimal hyperparameters (i.e., models) obtained across all 50 trials.
Most trials exhibit satisfactory training and validation loss curves, with rapid decreases during the initial epochs that stabilize near-zero around 50 epochs (Fig. 4a and b). Generally, trials with steeper initial declines and lower asymptotic losses tend to converge better and may have less risk of over昀椀tting (Singh et al., 2024). Most trials also stop early, around 40 to 50 epochs (Fig. 5c), suggesting they reach optimal validation loss relatively quickly, though a few continue up to 200 and 300 epochs, indicating more complex convergence patterns. The validation losses of most trials fall within a narrow range, with a signi昀椀cant concentration between 5E-5 and 1E-4 m 2 (Fig. 5d). The trial with the lowest validation loss (Trial 11) is achieved after 50 epochs (Fig. 5e). In this trial, the training-to-validation loss ratio initially oscillates above 1 in the early epochs (Fig. 5f), signaling some expected early divergence between the losses. However, this ratio stabilizes around epoch 10, trending closer to 1, which suggests a balanced performance between training and validation. Therefore, all trial ratios trending towards 1 underscore the model's consistency in generalizing across datasets (Fig. 5f), while ratios below 1 indicate over昀椀tting. Finally, we apply all models and evaluate their ability to accurately capture 昀氀ood events.
The attention vectors extracted by the LSTM inform how well the model integrates the in昀氀uence of water level variability in predicting 昀氀ood dynamics across the study area. Timesteps of high-water level variability (or abrupt changes) receive higher weights, indicating potential 昀氀ooding around the observation stations, while lower magnitudes suggest reduced 昀氀ood risk (Fig. 6). Based on the modulation approach, attention weights from individual stations are projected over the clusters containing those stations, allowing distinct and varying in-昀氀uences that correlate with the station's water levels. For instance, at timestep 130 h, a higher attention weight is applied to the cluster associated with station 4 (Galveston Bay Entrance, TX -NOAA 8771341) since the peak water level for Hurricane Ike occurs approximately 5 h earlier than at station 18 (Manchester, TX -NOAA 8770777) (Fig. 6a). However, attention weights for both station clusters at timestep 188 h during the Memorial Day 昀氀ooding are quite similar, re昀氀ecting the relatively low water levels at that time (Fig. 6b). This approach effectively accounts for lag-times in 昀氀ood dynamics across the larger domain, enhancing the model's ability to predict regional 昀氀ood variations accurately. It is important to note that the spatiotemporal analysis is conducted at the cell-size scale using the ConvLSTM model. Furthermore, the attention vector does not imply uniform 昀氀ooding across the entire cluster; rather, its in昀氀uence is applied to the spatially varied S. Daramola et al. Environmental Modelling and Software 191 (2025) 106499 values per timestep already estimated at individual cells by the ConvLSTM model. This means that each individual cell in the spatial map has a calculated value that varies across cells and at each timestep. The attention vector, although uniform over its cluster, emphasizes these varying values per cell with information about the events' evolution.
In addition, we compare the proposed framework with a model applying an attention vector to the entire domain to highlight its ability to accurately account for the timing differences in 昀氀ood dynamics over a large domain. The single attention vector re昀氀ects the combined in昀氀uence of all observation stations. The cell with the highest actual peak water depth across all timesteps within each cluster is selected and compared with its predicted data to observe any lag-or lead-times. Results of this analysis show a signi昀椀cant lag-time of about 5 h for the majority of clusters and a lead-time of about 2 h for the clusters when a single attention vector is used. However, correct peak timing is observed in two clusters without a color shade (Fig. 7a), suggesting where the maximum attention weight aligns with the peak 昀氀ood timing. This indicates that averaging attention across stations may only capture 昀氀ood dynamics correctly in clusters where timing coincides by chance. In contrast, incorporating cluster-based application of attention vectors demonstrates a signi昀椀cantly improved match with the actual 昀氀ood timing, with only 5 clusters within 1 h lead-and lag-times (Fig. 7b).
In a previous study, we demonstrated that deep learning models can be signi昀椀cantly enhanced by incorporating an attention mechanism (Daramola et al., 2025). This approach helps emphasize indicators that signal the presence and severity of 昀氀ooding events, rather than relying solely on repetitive patterns in training datasets or overlapping characteristics among various 昀氀ood drivers. Consequently, attention mechanisms were integrated into the architecture of the proposed model in this study to ensure its effectiveness during both the training and testing phases. Given the model's satisfactory performance, this technique is shown to bolster its ability to effectively identify 昀氀ood patterns across space and time. This improvement enables the proposed model to generalize effectively, indicating that the framework is robust and can be applied to other regions with varying distributions of observation stations or diverse hydrodynamic conditions.
Most of the DL models indicate that the proposed framework achieves satisfactory accuracy in predicting 昀氀ood dynamics in terms of 昀氀ood extent and inundation depth (Fig. S4 in the supplementary material). To showcase the model's ability, we present its performance for predicting the peak 昀氀ood depth maps of the three events. For Hurricane Beryl (Fig. 8), the peak 昀氀ood depth occurred on July 8, 2024, corresponding to the timestep with the maximum RMSE value. The best model slightly overpredicts 昀氀ood depth in coastal areas and underpredicts in inland areas, resulting in an overall slight underprediction with an error bias (E) lesser than 1. The false alarm ratio (F) identi昀椀es only 14 % of non-昀氀ooded areas are falsely classi昀椀ed as 昀氀ooded, with a critical success index (C) over 70 % of 昀氀ood maps (percentage of accurately predicted depth and extent). The model also achieved a high accuracy in identifying 昀氀ooded areas, with a hit rate (H) of 80 %. For the evolution of 昀氀ood depth in selected areas, the model demonstrates satisfactory performance across all timesteps with KGE, NSE, and R 2 metrics over 0.70.
At the peak 昀氀ooding on August 29, 2017, for Hurricane Harvey, the best model tends towards slight overprediction of the 昀氀ood map, with an error bias (E) greater than 1 (Fig. S10a and b, supplementary material), especially along the coast. However, 昀氀ood depth is underestimated in the inland areas (Fig. S10c). Similar to Hurricane Beryl, the 昀氀ood map has the maximum RMSE value. Nevertheless, C is over 60 %, indicating that most of the 昀氀ood map's depth and extent are accurately predicted. Additionally, over 80 % (H) of the 昀氀ooded areas are correctly identi昀椀ed. F is below 30 %, which means that the model is less likely to classify non-昀氀ooded areas as 昀氀ooded. Similarly, the model demonstrates satisfactory performance across all timesteps with KGE, NSE metrics over 0.70 for the evolution of 昀氀ood depth at selected areas. Just like the 昀椀rst two results, the time of peak 昀氀ooding (on September 14, 2021) corresponds to the timestep with the maximum RSME for Hurricane Nicholas. The best model tends to slightly underestimate 昀氀ooding in the inland region while overestimating 昀氀ood depth along the coast especially at the eastern area of the model domain (Fig. S11a and b, supplementary material). The model accurately predicts most of 昀氀ood map's depth and extent with a C of 64 %, while correctly identifying 昀氀ooded areas with 80 % accuracy (H). F is below 25 %, i.e., low tendency to incorrectly classify non-昀氀ooded areas as 昀氀ooded.
Regarding the peak 昀氀ood maps in Fig. 7, S10 and S11, the underestimation in inland areas is due to the omission of many interconnected rivers that could contribute to 昀氀ooding during these events. In contrast, the overestimation in coastal areas is caused by high attention vectors projected onto the clusters, as the associated stations are exposed to the ocean. Nevertheless, the model performances metrics are satisfactory for all 昀氀ood events and across all clusters (Fig. 9, S12 and S13).
The model utilizes a physically-based model simulation as the "ground truth" 昀氀ood extent and inundation depth during the compound 昀氀ood events. However, these model simulations slightly differ from actual observations introducing errors despite the thorough model calibration (MuÞ noz et al., 2024). The Delft3D-FM model is subject to several sources of uncertainty, including the initial condition, the forcing (or boundary) conditions, model parameters, and model structure (Abbaszadeh et al., 2022;MuÞ noz et al., 2022). Additionally, there may exist interpolation errors when resampling the spatially varying cell sizes of Delft3D-FM to a uniform 50-m DL model input's cell size. The resampling procedure was needed to enable correct model training and validation of the CNN architecture. Next, river discharges recorded at upstream river-gauge stations are applied along the river branches and used as an input feature in the training process. Ignoring the lag-time in peak 昀氀ow that may occur along the river branches might contribute to prediction errors, especially in large river deltas, estuaries, and bays. These uncertainties in the Delft3D-FM outputs can propagate to the DL model and affect its performance evaluation due to inherited biases, evident by the slight over-and underestimation of the comparison plots (Figs. S3-S6 in the supplementary material).
While the cluster-based attention vector for water level variability can account for the lag-time, some clusters may encompass areas with slight time variations, potentially affecting the model's accuracy. Nevertheless, all cells contained in each cluster are still spatially varying in magnitude, even though they follow the same timing pattern for the attention vector modulation during an event's evolution. When the average timing of all cells within the cluster is estimated, all the clusters have lead-and lag-times between 0 to approximately 2 h with respect to the actual peak timing. In other words, the more observation points for water levels, the more spatial timing variation will be captured; leading to a substantial improvement of the model's accuracy. Although incorporating attention mechanisms emphasizes certain weights and enhances generalization, some models produce straight-line predictions across all time steps due to L2 regularization effects in the architecture that penalize large weights (Figs. S7-S9). The latter occurs because a higher regularization coef昀椀cient can overly smooth the attention emphasis during predictions, preventing the model from capturing rapid 昀氀uctuations or noise in the data. Therefore, it is crucial to balance the strength of regularization technique within the architecture to enhance generalization without inducing under昀椀tting. To further reduce the error between model predictions and actual observations, future work will focus on introducing residual learning techniques (Tedesco et al., 2024;Zou et al., 2023). Residual learning enables the model to concentrate on learning the differences or residuals between its predictions and the true observations rather than attempting to model the entire mapping directly.
The performance metrics used in this study also have limitations. For example, RMSE is highly sensitive to large errors, so extreme 昀氀ood events can disproportionately in昀氀ate its value, potentially skewing perceptions of the model's overall performance. R 2 indicates explained variance but does not account for prediction bias, which is crucial for accurate 昀氀ood depth estimation. KGE and NSE, which are standard metrics in hydrology, help capture these biases during extreme events. Metrics like H, F, and CSI effectively assess 昀氀ood extent but may overlook the magnitude of depth predictions. E, however, is particularly signi昀椀cant for extreme events because it captures even small systematic over-and underpredictions. Despite their limitations, these metrics are used in a complementary manner and remain the most valid metrics for 昀氀ood analysis. Future studies could consider comparison with other state-of-the-art DL architectures and explore additional metrics or weighting schemes that prioritize performance during extreme conditions. The latter will ensure a balanced evaluation across all 昀氀ood scenarios. and E are 0 m, 1, 1, 0, 1, and 1, respectively. S. Daramola et al. Environmental Modelling and Software 191 (2025) 106499
In this study, we developed a DL framework to address the complexities of cyclone-induced compound 昀氀ood dynamics in Galveston Bay, TX. The framework combined LSTM and ConvLSTM model architectures, integrating spatial and temporal attention mechanisms to enhance the model's ability to capture nonlinear associations among 昀氀ood drivers such as precipitation, river discharge, and storm surge. By leveraging output data from a coastal hydrodynamic model (Delft3D-FM) and incorporating cluster-based attention vectors, the framework achieved high accuracy in replicating 昀氀ood dynamics in terms of inundation depth and 昀氀ood extent. Importantly, it accounts for the hysteresis in system behavior, leading accurate timing of 昀氀ood dynamics, particularly the propagation of 昀氀oods across a large domain, caused by extreme cyclonic events. We also show the capability of the DL framework to handle the spatial variability of dominant 昀氀ood drivers in coastal (storm surge) and inland areas (rainfall-runoff). This in turn demonstrates the potential of coupling DL architectures to enhance the prediction accuracy of cyclone-induced compound 昀氀ood dynamics, speci昀椀cally in complex coastal systems.
The model demonstrates a satisfactory performance in identifying 昀氀ooded areas, achieving an average hit rate of over 80 % while maintaining average false alarm rate below 25 %, even with the large domain size and resolution. Average KGE and NSE for the prediction of all events evolution are above 0.7. The study demonstrates that applying clusterbased temporal attention allows the model to focus selectively on different spatial regions, becoming more effective at capturing localized timing variations in 昀氀ood dynamics across a large domain. Despite these performances, incorporating more water level observation points and discharge observations along river channels could enhance spatial variability and improve the model's accuracy. The proposed coupled DL model holds promise for disaster preparedness and 昀氀ood mitigation in vulnerable coastal regions, offering a scalable approach adaptable to other cyclone-prone areas.