2.1. Related works

Previous studies have investigated meteotsunamis in different regions of the world, revealing both their common generation mechanisms and their distinct local features. Early works identified the principal resonant mechanisms responsible for amplification—Proudman resonance (Proudman, Reference Proudman1929), Greenspan resonance (Greenspan, Reference Greenspan1956), and harbor resonance (Denamiel et al., Reference Denamiel, Šepić and Vilibić2018). Numerical and observational analyses, such as those by Hibiya and Kajiura (Reference Hibiya and Kajiura1982) in Nagasaki Bay and Vilibić et al. (Reference Vilibić, Monserrat, Rabinovich and Mihanović2008) in the Adriatic, it showed how bathymetry and atmospheric forcing jointly shape meteotsunami behavior. More recent studies have explored meteotsunami dynamics in various basins, including the Adriatic Sea (Denamiel et al., Reference Denamiel, Šepić, Huan, Bolzer and Vilibić2019a, Reference Denamiel, Šepić, Ivanković and Vilibić2019b), the Yellow Sea (Kim et al., Reference Kim, Choi and Omira2022), the coast of Portugal (Kim and Omira, Reference Kim and Omira2021), and the United Kingdom (Thompson et al., Reference Thompson, Renzi, Sibley and Tappin2020; Lewis et al., Reference Lewis, Smyth, Williams, Neumann and Cloke2023). Together, these works demonstrate that meteotsunamis are more frequent and globally distributed than once believed (Gusiakov, Reference Gusiakov2021).

Numerical models based on the nonlinear shallow-water equations (NLSWE) remain the standard tool for simulating tsunami and meteotsunami propagation. Recent developments have also highlighted the potential of combining physical models with data-driven techniques to improve computational efficiency and interpretability (Fauzi and Mizutani, Reference Fauzi and Mizutani2020a; Andraud, Reference Andraud2025).

From the observational perspective, dense offshore networks such as DONET and S-Net in Japan (Baba et al., Reference Baba, Takahashi and Kaneda2014; Kaneda et al., Reference Kaneda, Kawaguchi, Araki, Matsumoto, Nakamura, Kamiya, Ariyoshi, Hori, Baba and Takahashi2015; Kubota et al., Reference Kubota, Saito, Chikasada and Sandanbata2021) and the DART system (Titov et al., Reference Titov, Gonzalez, Bernard, Eble, Mofjeld, Newman and Venturato2005) have significantly improved tsunami detection and early warning. However, their high cost and maintenance demands (Bernard and Titov, Reference Bernard and Titov2015) limit their feasibility for meteotsunami monitoring. This limitation motivates research into optimal sensor placement methods that can maximize information from a minimal number of stations.

In recent years, data-driven techniques have increasingly been applied to tsunami science. DA has been used to integrate model predictions with observations in real time (Kalnay, Reference Kalnay2002; Maeda et al., Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015; Byrne et al., Reference Byrne, Horsburgh and Williams2023). Similar hybrid approaches combining DA with machine-learning frameworks have also been explored for rapid tsunami wavefield reconstruction (Fauzi and Mizutani, Reference Fauzi and Mizutani2020b). For optimizing sensor networks, most tsunami studies have relied on Empirical Orthogonal Function (EOF) analysis to identify spatial patterns that capture the largest variance (Mulia et al., Reference Mulia, Gusman and Satake2017, Reference Mulia, Gusman, Williamson and Satake2019; Ren et al., Reference Ren, Wang, Wang, Zhao, Hu and Li2022). However, EOF assumes stationary linear modes, which may not fully describe the non-stationary and time-dependent characteristics of meteotsunamis forced by evolving atmospheric disturbances. To address this, DMD has emerged as a promising approach. DMD captures time-dependent structures and dominant frequencies within complex flow fields (H. Tu et al., Reference Tu, Rowley, Luchtenburg, Brunton and Nathan Kutz2014) and has recently been applied to atmospheric and oceanic datasets (Gugole and Franzke, Reference Gugole and Franzke2020; Li et al., Reference Li, Mémin, Tissot, Chapron, Crisan, Holm, Mémin and Radomska2023). Building upon these developments, the present study integrates DMD and DA to identify dynamically significant regions for offshore monitoring and to reconstruct meteotsunami wavefields from sparse observations.

2.2. Background concept

The numerical modelling of meteotsunamis is commonly based on the NLSWE with atmospheric pressure as an external forcing term. These equations describe the depth-averaged conservation of mass and momentum, allowing simulation of long-wave propagation in coastal and shelf regions. The sea-surface displacement is governed by hydrostatic balance, responding to spatial and temporal variations in air pressure through the inverse-barometer effect. When the propagation speed of the atmospheric disturbance approaches the wave celerity $ c=\sqrt{gd} $ , where $ d $ is the local ocean depth, Proudman resonance occurs, producing resonant energy transfer from air to sea. Additional amplification mechanisms, such as Greenspan and harbor resonances, depend on coastal geometry and bathymetry (Denamiel et al., Reference Denamiel, Šepić and Vilibić2018).

Atmospheric pressure forcing is typically represented using idealized Gaussian-shaped perturbations (Anarde et al., Reference Anarde, Cheng, Tissier, Figlus and Horrillo2021) or by interpolating observational data (Kim and Omira, Reference Kim and Omira2021; Kim et al., Reference Kim, Choi and Omira2022). These formulations enable the reconstruction of realistic atmospheric inputs for numerical simulations and are widely used for assessing meteotsunami generation and propagation.

To extract dominant spatiotemporal patterns from the simulated fields, DMD is employed as the primary decomposition method within the framework. Coupling DMD with sequential DA (Kalnay, Reference Kalnay2002; Maeda et al., Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015) allows model states to be updated using sparse observations, enhancing both reconstruction accuracy and monitoring efficiency.

These physical and methodological concepts form the foundation of the present framework, which combines NLSWE-based simulation, DMD analysis, and DA to improve meteotsunami detection and forecasting in data-sparse coastal regions. In the following, we describe the 2022 Ireland meteotsunami and associated modelling choices.

3.1. Numerical model description: governing equations

We use a numerical model based on a non-linear shallow water equation (NLSWE) with atmospheric pressure as the source model. NLSWE-based models have been widely used in meteotsunami studies. For example, Vilibić et al. (Reference Vilibić, Monserrat, Rabinovich and Mihanović2008), Šepić et al. (Reference Šepić, Rabinovich and Sytov2018a), and Kim and Omira (Reference Kim and Omira2021) employed depth-averaged NLSWE models with idealized atmospheric pressure forcing to simulate meteotsunami dynamics, validating their results against field observations. We consider a spherical coordinate system ( $ R $ , $ \theta $ , $ \lambda $ ) with the origin at the Earth’s center. The Earth is assumed to be a sphere, hence the radius $ R $ is constant and equal to the Earth’s radius, whereas $ \theta $ is latitude (measured southwards) and $ \lambda $ is longitude (measured eastwards). The governing equations of NLSWE are given by the following expressions:

(1) $$ \frac{\partial D}{\partial t}+\frac{1}{R\sin \theta}\frac{\partial (uD)}{\partial \lambda }+\frac{1}{R}\frac{\partial (vD)}{\partial \theta }=0, $$

(2) $$ \frac{\partial u}{\partial t}+\frac{u}{R\sin \theta}\frac{\partial u}{\partial \lambda }+\frac{v}{R}\frac{\partial u}{\partial \theta }+\frac{g}{R\sin \theta}\frac{\partial h}{\partial \lambda }+{C}_f\frac{u\sqrt{u^2+{v}^2}}{D}+\frac{1}{\rho R\sin \theta}\frac{\partial P}{\partial \lambda }=0, $$

(3) $$ \frac{\partial v}{\partial t}+\frac{u}{R\sin \theta}\frac{\partial v}{\partial \lambda }+\frac{v}{R}\frac{\partial v}{\partial \theta }+\frac{g}{R}\frac{\partial h}{\partial \theta }+{C}_f\frac{v\sqrt{u^2+{v}^2}}{D}+\frac{1}{\rho R}\frac{\partial P}{\partial \theta }=0, $$

(4) $$ {C}_f=\frac{g{n}^2}{D^{\frac{1}{3}}}. $$

Manning’s roughness coefficient $ n $ is used to calculate the non-dimensional frictional coefficient $ {C}_f $ . A constant value of $ n=0.025 $ m^1/3s⁻¹ is applied homogeneously in the whole computational domain. This value is widely used in ocean wave propagation modeling (Imamura et al., Reference Imamura, Yalçiner and Ozyurt2006). In equations (1)–(4), $ u $ and $ v $ are depth-averaged horizontal velocity fluxes along the latitude and longitude axes, respectively, $ P $ is air pressure, $ t $ is time, $ g $ is the gravitational constant, and $ D $ is the total water depth, which is a summation of water depth $ d $ and sea surface amplitude $ h $ . We assume that seawater density $ \rho $ ~ 1025 kg/m³. The edge of the computational domain is used as an open boundary and the Coriolis effect is neglected in this study. The computational time step should satisfy the Courant–Friedrichs–Lewy (CFL) condition, whereby the time step $ \Delta t $ must be equal to or smaller than the time required by a wave to travel on the spatial grid size.

3.2. Air pressure forcing

Air pressure modeling is a critical variable in meteotsunami simulations and represents a key distinction from seismic tsunami modeling. In seismically generated tsunamis, once the hypocenter coordinates are known, the source model is typically estimated using waveform inversion of either oceanic or seismic waves (Satake, Reference Satake1989; Gusman and Tanioka, Reference Gusman and Tanioka2013). This is possible not only due to the availability of observing sensors but also because seismic waves travel through the Earth’s layers and can be detected globally (Obara, Reference Obara2007; Gee and Leith, Reference Gee and Leith2011; Di Giacomo et al., Reference Di Giacomo, Storchak, Safronova, Ozgo, Harris, Verney and Bondar2014). However, estimating a tsunamigenic atmospheric perturbation, such as an air pressure jump, is challenging due to the limited spatial coverage of observational stations. As a consequence, there is significant uncertainty regarding the shape, propagation direction, and speed of the pressure perturbation. Gaussian shapes are commonly used to represent the shape of the pressure jump (Anarde et al., Reference Anarde, Cheng, Tissier, Figlus and Horrillo2021). A more realistic source model can be constructed by interpolating between observed pressure values (Kim and Omira, Reference Kim and Omira2021; Kim et al., Reference Kim, Choi and Omira2022).

Due to the limited information on the spatial distribution of air pressure, this study employs a Gaussian-shaped profile to estimate the characteristics of air–pressure disturbances.

(5) $$ {P}_{i,j,t}=\sum \limits_{k=1}^N{P}_k\mathit{\exp}\left[{\left(-\frac{x_{i,j} sin\phi +{y}_{i,j} cos\phi -{V}_ot}{\frac{\sigma_k}{2}}\right)}^2\right], $$

where $ P $ is air pressure on a two-dimensional grid, $ {P}_k $ is the initial amplitude of the $ k $ th pressure component (with a total of $ N $ components), $ x $ and $ y $ are longitude and latitude coordinates, $ \phi $ is the propagation direction of the pressure disturbance measured counterclockwise from the $ x $ -axis, $ {V}_o $ is the propagation speed, $ t $ is the numerical time step, and $ \sigma $ is the spread of the kth component, which controls the width of the Gaussian distribution. Similar to Anarde et al. (Reference Anarde, Cheng, Tissier, Figlus and Horrillo2021), we create a synthesized representation of the physical pressure perturbation by combining individual components into a smooth, continuous two-dimensional distribution. A total of four crests and troughs, each with distinct $ {P}_0 $ and $ \sigma $ values, are identified from the pressure records. The spatial distribution of the AWSs allows for a quantitative estimation of the direction and propagation speed of the air-pressure disturbance. Based on the locations of the AWSs (Figure 1) and the clear sequential pattern of maximum pressure jumps over time (Figure 3), we conclude that the pressure disturbances propagated in a northeastward direction, in agreement with Renzi et al. (Reference Renzi, Bergin, Kokina, Pelaez-Zapata, Giles and Dias2023). We assume that the pressure disturbance propagates at a constant speed and angle. This assumption is widely used in the previous study (Šepić et al., Reference Šepić, Vilibić, Rabinovich and Tinti2018b; Anarde et al., Reference Anarde, Cheng, Tissier, Figlus and Horrillo2021; Kim and Omira, Reference Kim and Omira2021).To determine the propagation speed and direction, the process typically begins with hypothetical values, followed by multiple trial-and-error iterations to best fit the observations. In this study, we first identified the time of maximum air-pressure perturbation recorded at each AWS along the track of the synoptic disturbance. The average propagation speed was then calculated as the mean distance between Sherkin Island and the other AWSs divided by the corresponding time lag between their peak-pressure records, yielding an average speed of 55.4 m/s. The propagation direction (44° from the west) was estimated from the linear alignment of these stations and the spatiotemporal sequence of the pressure peaks. The parameter estimation workflow is summarized in Figure 4, illustrating the stepwise process from AWS observations to the synthesized pressure field.

Figure 4. Atmospheric pressure parameter estimation workflow.

To illustrate the resonance condition, Figure 1 includes contours of the local long-wave phase speed $ c=\sqrt{gd} $ . The region where the atmospheric disturbance speed matches the oceanic wave celerity indicates where Proudman resonance is expected to occur—along the 50–60 m/s contour southwest of Ireland. Figure 5 presents the comparison between simulated and observed air pressure at the four AWS locations. To quantify the agreement between observed and simulated air pressure, we use the normalized cross-correlation coefficient (NCC), where values closer to 1 indicate a strong correlation. Additionally, we evaluate the root mean square error (RMSE). Overall, the simulated air pressure demonstrates good agreement with observations, yielding an NCC value of 0.83 and an RMSE of 0.38 hPa across all AWS locations.

Figure 5. Comparisons between simulated and observed air pressure at four AWSs.

3.3. Meteotsunami generation and propagation

At the beginning of the numerical simulation (at $ t=0 $ ), the air pressure disturbance is transformed into a two-dimensional grid to create an initial sea surface using hydrostatic approximation,

(6) $$ h=-\frac{\Delta P}{\rho g}. $$

The topography and bathymetry dataset are obtained from the GEBCO 2024 grid with a resolution of 15 arcsec (~450 m), see Figure 1. To satisfy the CFL condition, the numerical time step $ \Delta t $ is set to 1 s. We run the 12 h simulation time using simulated air pressure from 18 June 2022, 09:00 to 21:00 UTC to observe the dynamics of meteotsunami during the generation and propagation phases.

The snapshots of pressure and wave propagation over the computational domain are shown in Figure 6, whereas the comparison between observed and simulated waves at four coastal gauges is presented in Figure 7. Figure 6 shows that meteotsunami propagation is directly guided by pressure forcing along its propagation path. Our numerical model also captures wave amplification due to local bathymetry. Initially, the model generates only a few centimeters of sea surface displacement in the offshore region, which then propagates toward the coast, increasing in magnitude as a result of nearshore amplification processes. As shown in Figure 7, the simulated waveforms exhibit good agreement with observations at Castletownbere, although the model slightly overestimates both the maximum and minimum wave amplitudes. The simulated maximum amplitude is 0.17 m, overestimating observations only by 0.05 m, and the overall waveform pattern closely aligns with the observed data. Similarly, comparisons at Union Hall and Ballycotton show a strong correlation, despite a minor temporal shift in the timing of the maximum wave amplitude. The maximum wave amplitudes are 0.12 m at Union Hall and 0.15 m at Ballycotton, with overestimations of 0.015 m and underestimations of 0.01 m, respectively. At Wexford, however, the model overestimates the observed wave heights, suggesting potential discrepancies in model accuracy at this location, which will be discussed in Section 4. The maximum simulated amplitude at Wexford is 0.1 m, compared to 0.03 m in observations, indicating an overestimation of 0.07 m.

Figure 6. Snapshots of simulated air pressure and the influenced sea surface over computational domain.

Figure 7. Comparisons between simulated and observed waveforms at four tide gauges.

3.4. Computational setup

The numerical simulations were performed using Fortran for the nonlinear shallow-water solver and Python for the post-processing analyses, including DMD, K-means clustering, and DA. The Fortran code was compiled and executed on the Linux-based HPC cluster at Northumbria University. Each run used a single compute node (Intel Xeon 2.2 GHz CPU), and a 12-hour simulation of the 18 June 2022 event required about three hours of wall-clock time.

The computational complexity of the numerical solver scales linearly with the number of spatial grid points $ N $ expressed as $ O(N) $ , while the DMD analysis scales with both the number of spatial points and temporal snapshots $ M $ , approximately $ O\left(N\times M\right) $ . Data assimilation using Optimal Interpolation operates at $ O(N) $ per assimilation cycle. The numerical solver structure is summarized in Figure 8, and all simulation parameters are provided in Table 1 to ensure reproducibility. Detailed computational specifications—including hardware configuration, compiler version, and optimization flags, and software environment—are also reported in Table 1.

Figure 8. Nonlinear shallow-water equation solver structure.

Table 1. Numerical simulation parameters used in the meteotsunami modeling

4.1. Dynamic mode decomposition and sensitivity analysis

We use dynamic mode decomposition (DMD) for analyzing spatiotemporal patterns, and here we apply it to the 2022 Ireland meteotsunami case, where high-quality data offer a valuable opportunity to demonstrate its strengths. Currently, on the Irish side of the Celtic Sea, only two offshore buoys (M3 and M5) are located near the study area, with a 1-h sampling interval and positioned relatively far apart. This spatial and temporal sparsity is insufficient to capture meteotsunami dynamics. To improve detection and monitoring, observational stations should be strategically located in regions where the most energetic dynamics occur.

Previous studies have utilized empirical orthogonal function (EOF) analysis to determine the most energetic locations for seismically generated tsunamis (Mulia et al., Reference Mulia, Gusman and Satake2017, Reference Mulia, Gusman, Williamson and Satake2019; Ren et al., Reference Ren, Wang, Wang, Zhao, Hu and Li2022). EOF, essentially a principal component analysis, identifies static spatial patterns that explain the highest variance in the dataset. However, it assumes linear modes, which may not always adequately represent the evolution of wave propagation. In seismically generated tsunamis, wave propagation is primarily driven by water height gradients, which influence velocity changes. However, in meteotsunamis, waves propagate due to both water height gradients and atmospheric pressure forcing (Eqs. [2] and [3]). As a result, meteotsunamis exhibit distinct spatiotemporal variations compared to seismically generated tsunamis, which require a different approach.

To address this limitation, we employ DMD, which captures time-dependent patterns and frequencies within the dataset (Tu et al., Reference Tu, Rowley, Luchtenburg, Brunton and Nathan Kutz2014). Unlike EOF, DMD extracts modes that evolve dynamically over time, offering a more accurate representation of wave evolution. Here, we derive the DMD modes of meteotsunami amplitudes resulting from numerical simulations. For the DMD computation, tsunami amplitudes within the model domain are stored at 5-minute intervals. We use the first 4 h of the meteotsunami simulation, yielding a total $ m=48 $ snapshots. The DMD algorithm decomposes the spatiotemporal meteotsunami data into dynamic modes, each associated with specific frequencies and growth rates. Similar mode decomposition approaches have been used to isolate physically significant instabilities in nonlinear wave simulations, such as the dipole-shaped superharmonic unstable modes linked to wave breaking (Mansar et al., Reference Mansar, Turner, Bridges and Dias2025).

First, the simulated sea surface dataset $ \left\{{h}_1,\dots, {h}_m\right\} $ is arranged into two sequential matrices,

(7) $$ {X}_1=\left[{h}_1,\dots, {h}_{m-1}\right], $$

(8) $$ {X}_2=\left[{h}_1,\dots, {h}_m\right], $$

where each ith column represents a sea surface elevation at a given time interval. The goal of DMD is to find a best-fit linear operator $ A $ such that,

(9) $$ {X}_2\approx A{X}_1. $$

Since $ A $ is typically high-dimensional, a singular value decomposition is applied to $ X $ , yielding a reduced-order representation. The matrix $ X $ is factorized as

(10) $$ X=U\Sigma {V}^T, $$

where $ U $ and $ V $ are orthonormal matrices, and $ \Sigma $ is a diagonal matrix containing singular values. To obtain a low-rank approximation of $ A $ , we compute

(11) $$ \overset{\sim }{A}={U}^T{X}_2V{\varSigma}^{-1}, $$

where $ \overset{\sim }{A} $ is a reduced-order representation that captures the dominant system dynamics while mitigating computational complexity. Next, we solve the eigenvalue problem,

(12) $$ \overset{\sim }{A}W=W\Lambda, $$

where $ \Lambda $ contains the DMD eigenvalues, which define the oscillation frequencies and growth rates, and $ W $ contains the corresponding eigenvectors, representing temporal evolution. Finally, the true DMD modes, which describe the spatial structures of wave dynamics, are reconstructed using

(13) $$ \varphi ={X}_2V{\Sigma}^{-1}W. $$

These DMD modes reveal the dominant spatiotemporal patterns of meteotsunami evolution, providing a more dynamic characterization than traditional EOF analysis. By employing DMD, we capture the complex interactions between pressure forcing and ocean response over time, offering a more comprehensive approach for identifying key regions where offshore observational stations should be deployed.

Furthermore, since the DMD results consist of multiple spatial points, many of which may be redundant, we apply K-means clustering (Lloyd, Reference Lloyd1982) as an optimization step to group similar modes and remove redundancies. K-means clustering optimizes the selection of dominant modes by minimizing the intra-cluster variance, which is given by the objective function:

(14) $$ J=\sum \limits_{i=1}^k\sum \limits_{x\in {C}_i}{\left\Vert x-{\mu}_i\right\Vert}^2, $$

where $ k $ is the number of clusters, $ {C}_i $ represents the set of data points in cluster $ i $ , $ {\mu}_i $ is the centroid of cluster $ i $ , and $ {\left\Vert x-{\mu}_i\right\Vert}^2 $ denotes the squared Euclidean distance between a data point $ x $ and its cluster centroid. By iteratively updating centroids and reassigning points, K-means partitions the DMD modes into a set of representative clusters, reducing noise and redundancy in the final analysis. This process ensures that only the most significant and distinct spatiotemporal patterns are retained, improving the interpretability and efficiency of meteotsunami dynamic analysis.

The results of the first two spatial modes and their corresponding extrema are presented in Figure 9a,b. We focus on these two leading modes, as they capture the majority of the variance and represent the most energetic meteotsunami dynamics in the study area, primarily associated with meteotsunami propagation. While the first mode likely represents the largest-scale, slowest-changing features in the dataset, the second mode captures secondary variations, such as reflected waves, smaller oscillations, and bathymetric influences. The energy variance for these modes is similar, with the first mode representing 25.04% and the second mode for 22.68% of the total variance. Although subsequent modes also contribute to the overall variance (e.g., 15.96% for the third mode and 13.58% for the fourth mode), their contributions are significantly smaller, and higher modes primarily represent less energetic processes. For comparison, Figure 9d,e display the spatial modes and extrema derived from the traditional EOF analysis.

Figure 9. Spatial distributions and station selections based on DMD and EOF analysis. (a) First DMD mode and (b) second DMD mode, with extrema indicated by black circles. (c) Selected extrema from DMD modes identified as potential observational stations (black circles). (d) First EOF mode and (e) second EOF mode, with extrema indicated by black circles. (f) Selected extrema from EOF modes used as potential observational stations (black circles). (g) Predefined station locations resembling the S-Net array configuration, with approximately 50 km spacing. (h) Evolution of the greedy optimization process for DMD- and EOF-based station selection. Numbers above and below the markers indicate the stations removed at each iteration.

To identify optimal observation locations, we evaluated the extrema of modes 1 and 2 through an iterative, greedy optimization based on DA sensitivity testing. In each iteration, one candidate extremum was removed, and the DA performance was recalculated. The extrema whose exclusion caused the largest decrease in the overall assimilation score (i.e., the most dynamically significant locations) were retained, while those whose removal had minimal impact on the score were discarded. The performance metric was quantified using a scalar skill score that combines improvements in both amplitude accuracy and waveform coherence:

(15) $$ Score=w\left(1-\frac{RMSE}{RMSE_f}\right)+\left(1-w\right)\left(\frac{NCC-{NCC}_f}{1-{NCC}_{\mathrm{f}}}\right), $$

where $ {RMSE}_f $ and $ {NCC}_f $ denote the forward modeling baselines, and $ w=0.5 $ assigns equal weight to amplitude and phase coherence. Higher values of the score indicate better assimilation performance.

The evolution of this optimization process for both DMD and EOF extrema is shown in Figure 9h, where the score changes with the successive removal of candidate stations. This greedy elimination continued until all non-significant extrema were removed. In addition, extrema located less than 30 km from the coastline (D1, E1, and E13) were excluded. The remaining candidate stations should have a space of approximately 85–100 km apart, ensuring a balanced spatial distribution that enhances meteotsunami monitoring and forecasting (Heidarzadeh et al., Reference Heidarzadeh, Wang, Satake and Mulia2019; Ren et al., Reference Ren, Wang, Wang, Zhao, Hu and Li2022).

The final DMD-based configuration initially identified the last seven candidate stations (D2, D3, D11, D12, D15, D16, and D20), where the assimilation score generally decreased as stations were removed. However, stations D11 and D12 were manually excluded to maintain adequate spacing between observation points. The resulting five optimized stations—D2, D3, D15, D16, and D20—correspond to the final five iterations in Figure 9h (see Figure 9c for their spatial locations).

Similarly, the EOF-based configuration retains three key stations—E2, E5, and E23—as shown in Figure 9f. Although the score shows a minor increase when E2 was excluded (Figure 9h), this likely reflects the sensitivity of the greedy procedure rather than a physically meaningful improvement. Station E2 is spatially well separated from the other EOF extrema and extends the observational coverage toward the eastern shelf, ensuring that the optimized network samples a broader portion of the meteotsunami propagation path. Therefore, E2 was retained to maintain spatial balance and avoid excessive clustering of observation points in the southern region.

4.2. Meteotsunami data assimilation

To model the meteotsunami monitoring system, we apply data assimilation using optimized station locations derived from DMD and K-means clustering to assess their performance. Sequential data assimilation (DA) techniques have been widely used in weather simulation (Kalnay, Reference Kalnay2002). Maeda et al. (Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015) introduced a DA approach for real-time tsunami wavefield simulation using the optimal interpolation method. Byrne et al. (Reference Byrne, Horsburgh and Williams2023) applied a similar technique to transform tide gauge records into a regional storm surge model. Optimal interpolation is computationally less demanding than more advanced methods, such as the ensemble Kalman filter (Yang et al., Reference Yang, Dunham, Barnier and Almquist2019), provided that the system remains linear. Despite its simplicity, this DA approach has demonstrated excellent agreement with both observed tsunami data and synthetic cases. To the best of our knowledge, this is the first study to implement data assimilation for meteotsunamis, with the goal of identifying potential offshore observational station locations.

In the numerical simulation, the meteotsunami wavefield at the nth time step is represented as $ {\boldsymbol{x}}_{\boldsymbol{n}}\left(\eta \left(n\Delta t,x,y\right),M\left(n\Delta t,x,y\right),N\left(n\Delta t,x,y\right)\right) $ , where $ \eta $ is the meteotsunami height, $ M $ and $ N $ are the depth-integrated meteotsunami fluxes in the $ x $ and $ y $ directions and $ \Delta t $ is the time step. The DA method can be expressed as

(16) $$ {\boldsymbol{x}}_n^f=\boldsymbol{F}{\boldsymbol{x}}_{n-1}^a, $$

(17) $$ {\boldsymbol{x}}_n^a={\boldsymbol{x}}_n^f+\boldsymbol{W}\left({\boldsymbol{y}}_n-\boldsymbol{H}{\boldsymbol{x}}_n^b\right). $$

At each time step, the forecasted meteotsunami wavefield, $ {\boldsymbol{x}}_n^f $ , is obtained by solving the non-linear long wave equation using the assimilated wavefield from the previous time step, $ {\boldsymbol{x}}_{n-1}^a $ . Here, $ \boldsymbol{F} $ is the propagation matrix corresponding to the two-dimensional NLSWE model. The residual—defined as the difference between the observed meteotsunami amplitude at the observational stations ( $ {\boldsymbol{y}}_n $ ) and the simulation output ( $ \boldsymbol{H}{\boldsymbol{x}}_n^b $ )—is calculated, where $ \boldsymbol{H} $ is a vector with entries of 1 at the observational stations and 0 elsewhere, used to extract the forecasted meteotsunami height. These residuals are then multiplied by a smoothing matrix, $ \boldsymbol{W} $ , which is crucial for bringing the assimilated wavefield closer to the observed meteotsunami field. The smoothing matrix is determined by solving the following equation:

(18) $$ \boldsymbol{W}=\boldsymbol{P}{\boldsymbol{H}}^{\boldsymbol{T}}+{\left(\boldsymbol{HP}{\boldsymbol{H}}^{\boldsymbol{T}}+\boldsymbol{R}\right)}^{-\mathbf{1}}, $$

where $ \boldsymbol{P}=\left\langle {\varepsilon}^f{\varepsilon}^{fT}\right\rangle $ and $ \boldsymbol{R}=\left\langle {\varepsilon}^O{\varepsilon}^{OT}\right\rangle $ are the error covariance matrices of the forward simulation and the observations, respectively. $ {\varepsilon}^f $ and $ {\varepsilon}^O $ are the errors of the forward simulation and observations, while $ {\varepsilon}^{fT} $ and $ {\varepsilon}^{OT} $ are the corresponding transpose matrices. In practice, $ \boldsymbol{P} $ is modelled using an isotropic Gaussian correlation $ {P}_{ij}={\sigma}_p^2\mathit{\exp}\left[-\frac{d_{ij}^2}{L^2}\right] $ , where $ {d}_{ij}^2 $ is the distance between grid points $ i $ and $ j $ , $ {\sigma}_p^2 $ is the background-error variance, and $ L $ is the characteristic distance controlling the spatial influence of each observation. Observation errors are assumed to be uncorrelated, so $ R={\sigma}_r^2I $ . Following Maeda et al. (Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015), the relative error ratio is expressed by the nondimensional parameter $ \rho $ , defined as the root-mean-square observation error normalized by the forecast error ( $ \frac{R}{P}={\rho}_d^2 $ ). Similarly, we adopted $ {\rho}_d=1 $ implying equal confidence in model and observations. Both covariance matrices assume Gaussian-distributed errors with a characteristic distance of 20 km (Maeda et al., Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015; Wang et al., Reference Wang, Satake, Maeda and Gusman2018). The assimilation cadence was identical to the numerical time step ( $ \Delta t $ ) used in the forward model. By iteratively applying Eqs. (16) and (17), the tsunami wavefield is effectively assimilated, enabling the generation of forecasted meteotsunami waveforms at any location within the model domain during and after the assimilation process. The assimilation was performed sequentially at each model time step corresponding to the available observation interval, with an assimilation cadence of 1 min. The complete data assimilation procedure, including the optimal interpolation scheme based on Maeda et al. (Reference Maeda, Obara, Shinohara, Kanazawa and Uehira2015), is visualized in Figure 10. This sequential algorithm updates the wavefield at each time step using observations from the optimized station network.

Figure 10. Flowchart of the sequential data assimilation procedure using optimal interpolation.

For a clear performance comparison, we run DA simulations using three sets of observational station configurations: DMD-generated, EOF-generated, and predefined stations. The predefined stations, shown in Figure 9g, represent a dense network of observational stations with a spatial interval of approximately 50 km as well as for the approximate closest distance to the coast, resembling the dense S-Net network. We then compare the resulting waveforms at tide gauges and the maximum meteotsunami amplitudes to evaluate the performance of each approach. As shown in Figure 11, the assimilated waveforms from the DMD-generated, EOF-generated, and predefined stations exhibit good agreement with forward modeling results and observations. Even though EOF-generated stations exhibit similar statistical performance to DMD-generated, EOF tends to underestimate wave amplitudes. For instance, at Castletownbere, the maximum wave amplitudes are better captured by DMD-generated stations. The NCC values at this gauge are 0.58, 0.6, and 0.63 for the DMD, EOF-optimized, and predefined station scenarios, respectively, with a similar RMSE of 0.05 m. Similarly, the assimilated waveforms from DMD-generated stations at Union Hall, Ballycotton, and Wexford demonstrate good agreement with observations, with NCC values of 0.51, 0.56, and 0.60 and RMSE values of 0.06 m, 0.08 m, and 0.02 m, respectively. Detailed information including DA configurations are summarized in Table 2.

Figure 11. Comparisons between simulated, observed, and assimilated waveforms at four tide gauges.

Table 2. Data assimilation parameters and performance metrics for all station configurations

We also compared the assimilated maximum wave amplitudes to those from the meteotsunami forward simulation, since survey data for maximum wave height along the southern coast of Ireland is not available. The predefined station network shows good agreement with the simulated results, likely due to its wide spatial coverage, which enables better capture of meteotsunami dynamics. Comparison between assimilation and simulation near the coastline (Figure 12) further confirms that the predefined scenario closely aligns with the forward modeling. Similarly, the DMD-generated stations also show comparable maximum wave amplitudes, although some underestimation is observed in the longitude range between −7.5°W and − 6.5°W, which is likely due to a lower density of observation points in that area. Nevertheless, this underestimation is minimal and does not significantly impact the overall accuracy of the assimilation results. In contrast, the EOF-based stations exhibit more pronounced underestimation in both the −7.5°W to −6.5°W and the southwestern region between 51°N to 52°N and − 11°W to −9°W, which is attributed to the absence of EOF-generated stations in these areas.

Figure 12. Comparison of maximum meteotsunami amplitudes along the southern coast of Ireland from (a) the forward modeling simulation, (b) data assimilation (DA) using five DMD-optimized observational stations, (c) DA using three EOF-optimized observational stations, and (d) DA using 16 predefined observational stations. Insets above and to the right of (b), (c), and (d) show maximum wave amplitudes at coastal locations, highlighting differences in spatial accuracy.