Evaluating essential processes and forecast requirements for meteotsunami-induced coastal flooding

Meteotsunamis pose a unique threat to coastal communities and often lead to damage of coastal infrastructure, deluge of nearby property, and loss of life and injury. The Great Lakes are a known hot-spot of meteotsunami activity and serve as an important region for investigation of essential hydrodynamic processes and model forecast requirements in meteotsunami-induced coastal flooding. For this work, we developed an advanced hydrodynamic model and evaluate key model attributes and dynamic processes, including: (1) coastal model grid resolution and wetting and drying process in low-lying zones, (2) coastal infrastructure, including breakwaters and associated submerging and overtopping processes, (3) annual/seasonal (ambient) water level change, and (4) wind wave-current coupling. Numerical experiments are designed to evaluate the importance of these attributes to meteotsunami modeling, including a “representative storm” scenario in the context of regional climate change in which a meteotsunami wave is generated under high ambient lake-level conditions with a preferable wind direction and speed for wind-wave growth. Results demonstrate that accurate representation of coastal topography and fully resolving associated hydrodynamic processes are critical to forecasting the realistic hazards associated with meteotsunami events. As most of existing coastal forecast systems generally do not resolve many of these features due to insufficient model grid resolution or lack of essential model attributes, this work shows that calibrating or assessing existing forecast models against coastal water level gauges alone may result in underestimating the meteotsunami hazard, particularly when gauging stations are sparse and located behind harbor breakwaters or inside estuaries, which represent dampened or otherwise unrepresentative pictures of meteotsunami intensity. This work is the first hydrodynamic modeling of meteotsunami-induced coastal flooding for the Great Lakes, and serves as a template to guide where resources may be most beneficial in forecast system development and implementation.


Introduction
The Laurentian Great Lakes are influenced by a variety of mid-latitude weather systems, from extratropical cyclones to localized convective storms. Combined with their sea-like characteristics, such as the immense sizes, great depths and steep bathymetric gradients, the Great Lakes have long been referred to as ''inland seas", prone to rolling waves, strong currents, storm surges and other hydrodynamic hazards. Over the last two decades, the Great Lakes water levels have swung from record lows to extreme highs (Gronewold and Rood 2019). Higher water levels, along with more intense storms (Feng et al. 2016;Jabbari et al. 2021) due to the hydrologic intensification that accompanies climatic warming trends, have further exacerbated coastal flooding hazards that endanger boaters, beachgoers, and caused more severe damage to coastal infrastructure, communities and ecosystems. Recent record-breaking high lake levels in 2020 across the Great Lakes calls for the urgent need for a capable modeling framework to predict coastal flooding events and to better prepare coastal communities for emergency management and development planning.
While various mechanisms, such as storm surges and seiches, can result in coastal flooding in the Great Lakes and coastal oceans, meteotsunamis are an important phenomenon that have caused disastrous damage to coastal property and loss of life due to their significant runup and associated strong currents (As-Salek and Schwab 2004;Šepić et al. 2015;Linares et al. 2019;Vilibić et al. 2021). These meteorologically induced water waves are similar to seismic tsunamis in spatial and temporal characteristics, limited to the frequency band of wave periods between 2 minutes to 2 hours, but are mainly caused by atmospheric pressure and wind perturbations associated with fast-moving weather events including severe thunderstorms, squalls and storm fronts (Vilibić et al. 2021). These atmosphericdisturbance-generated waves are often amplified by different resonance mechanisms such as the Proudman resonance, Greenspan resonance, shelf resonance and harbor resonance (Monserrat et al. 2006). Recent studies showed meteotsunami events with heights larger than one foot, a potentially dangerous magnitude, occur an average of 106 times per year, which is much higher than previous estimates, throughout the Great Lakes region, flooding coastal communities and causing dangerous rip currents (Bechle et al. 2016;Linares et al. 2019).
While early studies (e.g., Edwing et al. 1954;Donn 1959) suggested that meteotsunamis in the Great Lakes were primarily driven by atmospheric pressure perturbations; more recent studies identified that, depending on meteorological conditions, both atmospheric pressure and wind perturbations can be essential factors to influence meteotsunami magnitudes in the Great Lakes (Bechle and Wu 2014;Linares et al. 2019). Recently, a highamplitude pressure-driven meteotsunami occurred on April 13, 2018, when meteotsunami waves struck the Michigan coastline near Ludington, Michigan (Anderson and Mann 2020). Sitting on the east coast of Lake Michigan, the Ludington shoreline is characterized by sandy dunes and beaches with shallow nearshore water. The Ludington Harbor is protected by two harbor breakwaters open to the west, with maximum water depths of approximately 7 m in the harbor. On the shore, the Ludington region has a wide low-lying zone with sandy beaches and dunes (elevation lower than 3 m relative to low water datum) (Fig. 1). During this event, the harbor breakwaters were overtopped with the incident waves, and flood waters inundated the shoreline and nearby city streets. Damage to public docks and lakefront cottages were also reported.
Using surface meteorological conditions generated from a Weather Research and Forecasting (WRF; Powers et al. 2017) model simulation, Anderson and Mann (2020) examined 1 3 the performance of the National Oceanic and Atmospheric Administration (NOAA) Lake Michigan-Huron Operational Forecast System (LMHOFS; Kelley et al. 2020), an operational forecast model based on Finite Volume Community Ocean Model (FVCOM; Chen et al. 2006), in simulating the meteotsunami event. While the simulated results by LMHOFS show good agreement with observation in terms of timing of the meteotsunami arrival at several coastal locations, the present modeling framework falls short in a few important aspects relevant for coastal hazard prediction. First, with the model grid resolution of roughly 500 m in the shoreline, the LMHOFS model mesh only extends to outside of the breakwaters while the nearest water level gauge to Ludington is located inside the Ludington estuary (Fig. 2), which makes it infeasible to truly calibrate and validate the model's capability to predict water levels inside the harbor and Ludington estuary. Second, The green dot indicates the closest model grid to the gauge station 9,087,023 (red dot). The water level at the green dot was selected to compare with the gauged water level at station 9,087,023 in the original LMHOFS in Anderson and Mann (2020) the relatively coarse model grid resolution and the lack of inland coverage of the mesh and necessary numerical modules for nearshore processes do not permit the model to capture critical hydrodynamic details during this meteotsunami event, such as the wave signal into the harbor/estuary, the dynamics related to wetting and drying, wave runup, and inland inundation conditions. These deficiencies can result in a mischaracterization of meteotsunami intensity and impact, whereby an early warning system could provide incomplete or inaccurate description of the coastal hazard. Resolving these deficiencies and understanding the full extent of the meteotsunami impact are integral to developing a robust coastal flood forecast system.
In this study, we developed a high-resolution coastal model for the portion of east coastal water (86.42°W-86.60°W; 43.64°N-44.08°N) that covers the Ludington region ( Fig. 2) to resolve coastal structures and coastal topography for accurate simulation of nearshore hydrodynamics, including water overtopping the breakwaters, wetting and drying processes over the low-lying land, and wave-water interactions. Taking advantage of the unstructured grid mesh, we integrated the high-resolution coastal model mesh into LMHOFS (Hereafter referred to as Ludington-High Resolution or Ludington-HR). In such a way, the Ludington-HR model also covers the entire Lake Michigan-Huron ( Fig. 2a) as LMHOFS, thus providing representation of large-scale background circulation and remote forcing that drives Ludington nearshore water movement. This study revisits the 2018 meteotsunami event using the Ludington-HR and aims to identify and evaluate the required model capacities for forecasting dangerous nearshore conditions and flooding potential. The remaining sections are organized as follows. The Ludington-HR and the design of the numerical experiments are described in Sect. 2. In Sect. 3, the model performance in simulating the meteotsunami event and resulting flooding are presented. Impacts on the flooding of various factors, including model mesh grid resolution, coastal engineering structure, water level change, wave-current interaction, and atmospheric pressure perturbation, are examined. The conclusions are summarized in Sect. 4.

Hydrodynamic model
Both LMHOFS and the Ludington-HR are developed based on the hydrodynamic model FVCOM. FVCOM is a three-dimensional (3D) free surface, primitive-equation model that solves the momentum, continuity, temperature, salinity, and density equations and is closed physically and mathematically using the 2.5 level turbulence submodel for vertical mixing (Mellor and Yamada 1982) and the Smagorinsky formulation for horizontal diffusion (Smagorinsky 1963). FVCOM is solved numerically using the finite volume method in the integral form of the primitive equations over a horizontal unstructured triangular grid mesh, and the vertical dimension is represented by generalized terrain-following coordinates. The advantage of an unstructured grid for shoreline fitting and the flexibility of local mesh refinements makes it popular in applications to coastal waters and Great Lakes, in both stand-alone hydrodynamic modeling Anderson et al. 2018;Huang et al. 2019;Ye et al. 2019Ye et al. 2020 and coupled with other models such as water quality models (Xue et al. 2014a;Rowe et al. 2017) or regional climate models (Xue et al. 2014b(Xue et al. 2017Xue and Eltahir 2015).

3
The model mesh and bathymetry of Ludington-HR were developed based on those of LMHOFS and include local refinements. Specifically, the LMHOFS has an unstructured horizontal grid of 90,806 nodes and 171,377 elements, with horizontal resolution ranging from ~ 200 to 500 m near the coast to 2.5 km offshore (Anderson et al. 2018). In comparison, the Ludington-HR refines the original model grid of LMHOFS in the east coast around Ludington with a much higher horizontal resolution of 10 m (Fig. 2b), increasing the grid elements from 1,880 to 96,916 in the refinement area. The unstructured design allows the finer mesh to gradually relax to the LMHOFS grid for the seamless regional refinement (Fig. 2a). Vertically, LMHOFS and Ludington-HR share the same 21 uniform, terrain-following, sigma layers. In the Ludington-HR model, the bathymetry and topography in the Ludington region were updated with data from NOAA Electronic Navigational Charts (NOAA ENC), (https:// encdi rect. noaa. gov/) that contain high-resolution water depth information including inside channels and harbors, and 2012 USACE NCMP Topobathy Lidar data: Lake Michigan (https:// coast. noaa. gov/ digit alcoa st/ data/) nearshore topography, which has 1 m resolution that is sufficient to resolve the complexity at the land-water interface.
In addition to the higher model grid resolution, bathymetry, and shoreline, the Ludington-HR model is configured with two extra modules-the dike-groin module and the wetting/drying module-both are critical for accurate coastal inundation simulations. The dike-groin module was first developed and introduced by Ge et al. (2012), enabling the Ludington-HR model to represent the breakwaters by allowing water exchange over the breakwaters and blocking flow below the submerged structures. Triangular elements were generated along a breakwater on both sides. Figure 3a shows a segment of the breakwater (red line) and its surrounding model triangular elements. In FVCOM, the scalar variables such as water surface elevation (ζ) are designated at the triangle vertices. They are calculated by net flux through the Tracer Control Element (TCE), a section enclosed by the surrounding triangle centroids and the middle points of triangle sides (light-blue shaded regions). For a breakwater, the TCE is divided into two elements (Fig. 3b), and calculates For the water below the breakwater, the treatment is similar to the solid boundary condition, which ensures no flux normal to the wall. When the breakwater is temporarily submerged, the fluxes of the elements of the water layer above the breakwater are combined to determine the water level and currents. When water overtopping occurs with the rising water level on one side of the breakwater (e.g., surges, tides, and tsunamis), the total volume of spilled water and associated water levels on both sides of the breakwater are calculated under the local and global mass conservation (see Ge et al. 2012 for detailed model development and validation).
To simulate the water transport flooding onto and draining out of the low-lying coastal zone, the wetting/drying module (Chen et al. 2008) was also enabled. In this method, a viscous sublayer with a thickness D min is used to avoid the occurrence of a singularity when the local water depth approaches zero during the wetting and drying process. When the model simulated water column at a given triangle node is less than D min = 5 cm, the node is treated as dry. Subsequently, for a triangular cell with three nodes i, j, and k; the wet/dry condition is determined by the wet/dry conditions of the three nodes using the following criteria: where H(i,j,k) and ζ(i,j,k) are the bathymetry (negative value overland) and the surface water elevation at nodes i, j, and k. When a triangular cell is treated as dry, the velocity (u, v), which is configured at the centroid of this triangle. (Fig. 3) is specified to be zero and no flux is allowed through the three side boundaries of this triangle. This triangular cell is then removed from the flux calculation in the TCEs. The wetting/drying treatment has been validated for both idealized and realistic estuarine cases with detailed discussion by Chen et al. (2008).
Lastly, the complexity of coastal and nearshore hydrodynamics also lies in the fact that multiple dominant processes interact with each other and form strong nonlinearity. For example, the wave radiation stresses influence nearshore currents and water levels, and water level fluctuations also affect wave propagation and dissipation. To account for wave contributions to coastal flooding, we dynamically coupled the hydrodynamic model FVCOM and the wave spectral model SWAN (Simulating Waves Nearshore).
SWAN is a third-generation spectral wave model developed at Delft University of Technology that computes random, short-crested wind-generated waves in coastal regions and inland waters (http:// swanm odel. sourc eforge. net/). It solves the evolution equation of wave action density and accounts for various wave energy sources and sinks, including wave generation by wind, wave decay due to white capping, bottom friction, and depth-induced wave breaking, as well as energy redistribution through nonlinear wave-wave interactions. The model has been recognized as a reliable coastal community wave model which has been widely used for wave hindcasting and forecasting in coastal and inland waters (Rogers et al. 2003(Rogers et al. 2007Niroomandi et al. 2018). In the SWAN wave simulations, the computational mesh was curvilinear and consisted of 679 × 1073 grid cells, which gave a horizontal resolution of 10-20 m around the Ludington region (Fig. 4). The spectral domain was discretized into 12 directions with 30° intervals and 31 frequency bands from 0.0521 to 1.0. The breakwaters were modeled by a subgrid approach in SWAN as line structures, as they usually have a transversal area that is too small to be resolved by the model grid. Breakwaters will reduce the wave height of waves propagating through or over the structures and cause waves to be reflected. These effects were accounted for in the simulations. The wave transmission coefficient was calculated by the Goda formula (Goda et al. 1967). The reflection coefficient was set to be a constant of 0.5. Wave diffraction process was activated in the simulation.
In the coupled FVCOM-SWAN framework, the OASIS-MCT coupler (Craig et al. 2017) is implemented to exchange information between the two models. By doing so, modifications of the model code structure in each model were minimized and it is also much more efficient with respect to keeping model constituents updated to the relatively new versions (FVCOM4.1 and SWAN v41.01). In the coupled system, FVCOM and SWAN are integrated forward simultaneously, and the coupler passes the SWAN-simulated significant wave height, mean wave direction, mean wavelength, and peak wave period to FVCOM to calculate the radiation stress for resolving the wave-induced momentum. The coupler passes the FVCOM simulated free surface elevation and currents to the SWAN for the instantaneous water depth and relative wind speed for wave calculation. Also, in the hydrodynamic stand-alone simulations, the bottom stresses ( bx , by ) are calculated from a quadratic expression bx , by = C d are the x and y components of bottom current velocities. The drag coefficient C d is formulated by matching a logarithmic bottom boundary layer to the model at a height Z ab above the bottom as , 0.0025) where K = 0.4 is the von Kármán constant and Z 0 is the bottom roughness parameter. In the coupled FVCOM-SWAN simulations, turbulent wave-current bottom boundary layer (BBL) flows and combined bottom shear stresses due to the presence of waves and wave-current interactions are calculated with the BBL model proposed by Madsen (1994).

Atmospheric forcing
The Ludington-HR model was spun up from March 1, 2018, initialized from the NOAA LMHOFS, and driven by hourly meteorological output from High-Resolution Rapid Refresh (HRRR). The Ludington-HR simulation during the meteotsunami event on April 13, 2018, was driven by surface meteorological forcing from a WRF simulation with 7-km grid spacing and a two-minute output interval to capture the high-frequency variation of the barometric pressure gradient during the event, as described in Anderson and Mann (2020). However, we found an adjustment of underestimated barometric pressure from the WRF simulation is necessary to capture the pressure-driven meteotsunami event as described by the water level analysis below.
As shown in Fig. 2, due to the coarse model resolution, Anderson and Mann (2020) had to use the water level at the nearest LMHOFS model node (the green dot in Fig. 2) to compare with the gauged water level at National Ocean Service (NOS) observation station 9,087,023 (the red dot in Fig. 2) located far behind harbor breakwaters inside the drowned river mouth estuary. This geographic discrepancy makes it impossible to truly verify the LMHOFS performance in simulating coastal water levels as well as those inside the estuary and harbor. Notice that the gauged water level at station 9,087,023 had a maximum value of 177.1 m (i.e., a 0.3 m water level rise relative to the lake level of 176.8 m prior to the event) and LMHOFS predicted a maximum water level of 177.3 m during the event. It would appear that the LMHOFS overestimated the water level rise if not for the geographic difference between the LMHOFS model node and gauged location (Fig. 2). In fact, the LMHOFS simulation in Anderson and Mann (2020) has underestimated the water level, based on several sources of evidence. First, while recorded water level rise were just up to 0.3 m at station 9,087,023, the National Weather Service (NWS) office in Grand Rapids received reports of water level fluctuations of 2 m recorded just outside of the Ludington harbor, where divers were performing maintenance on a water intake at the time of the event. Second, the Lake Level Viewer operated by NOAA (https:// coast. noaa. gov/ llv/) provides the relationship of the static water level and the coastal inundation in the great lakes (Fig. 5). It shows that at a water level of 177.3 m as simulated by LMHOFS, no coastal inundation would occur. Even when the water level reaches 177.9 m, the inundation would occur but not be severe enough to flood onto the street, as reported for this event and documented with photographic and video evidence. Using the Lake Level Viewer, the water level would need to reach 178.2 m in order to create the flooding extent reported from local authorities. This is also consistent with the reported water level fluctuations of ~ 2 m outside the harbor.
Using the above water level and inundation analysis, along with the fact that the barometric pressure jump generated by the WRF simulation (dotted red line, Fig. 6a) is noticeably smaller than observation (black line, Fig. 6a), we assume the underestimated water level from the LMHOFS simulation in Anderson and Mann (2020) is due to an underestimated atmospheric pressure in the meteorological forcing (Fig. 6). From the atmospheric perspective, the underestimation of the atmospheric pressure jump in this high-amplitude pressure-driven meteotsunami event was the main reason that was responsible for underestimated magnitudes of meteotsunami waves traveling in the lake, which would not be able to induce coastal flooding as opposed to the actual condition. The purpose of this study is to, from the perspective of lake hydrodynamic processes and modeling, identify those critical hydrodynamic model attributes that are essential to resolve meteotsunami-induced flooding. Therefore, an adjustment of the atmospheric pressure was made as follows, where P t i is the atmospheric pressure at model grid, time t , and P t i is the temporal mean of P t i during the time 12:00-16:30 (GMT). Hence, ΔP t i is the temporal variation of P t i relative to its meanP t i . The variation of the adjusted pressure AP t i is amplifed by a factor of AF(= 2.9) so that the atmospheric pressure jump is elevated to the observed magnitude (blue line, Fig. 6a). The original and amplified (hereafter referred to as AP1 and AP2) spatial patterns of atmospheric pressure are shown in Fig. 6b, c, both show the traveling atmospheric inertia-gravity waves. We acknowledge the imperfection of this empirical adjustment, yet it serves well for our purpose to stay focused on evaluating essential hydrodynamic processes and identifying key hydrodynamic model attributes without being diverted to re-develop or recalibrate the atmospheric forecasting model, which is an undoubtedly important component in a real forecasting system but far beyond the scope of this study.

Numerical experimental design
To analyze the impacts of model resolution, coastal structure, water level change, windinduced waves, and atmospheric pressure perturbation on the event, a control run and six process-oriented numerical experiments were designed. The control run is a hydrodynamic stand-alone simulation that incorporates the best hydrodynamic model configuration available, as described in Sects. 2.1 and 2.2, including the high-resolution model grid, updated model bathymetry and topography, enabled dike-groin module and wetting-drying module, and adjusted barometric pressure forcing. The breakwaters are configured with a crest elevation of 2 m above the low water datum (176 m). The design of four sets of experiments is briefly summarized below for an overview. Further elaboration on the design of each case are presented in Sects. 3.2, 3.3 and 3.4: 1. Impact of breakwaters: case (BW0) was configured the same as the control run but without the breakwaters; case (BW1) was configured the same as the control run but with increased crest elevation to 3 m above the low water datum (176 m). 2. Impact of lake level: case (LL) was configured as the same as the control run, but the lake mean lake level was increased from 176.8 m in the control run to 177.6 m, representing the high water level observed in 2020 and serving as a sensitivity analysis of the inundation to natural lake level variation. 3. Impact of wind-induced waves: the wave-current coupled run (WC) is the case that dynamically couples the hydrodynamic simulation of control run (including breakwaters) with SWAN (including breakwaters). 4. A "representative storm" scenario: it integrates the meteotsunami, high water level, as well as favorable wind for wind-wave development along the east coast. Two cases were configured without (RS1) and with (RS2) wave-current coupling. Both were tested with the idealized southwesterly wind to favor the wind-induced wave growth around the east coast combining with the high water level observed in 2020 as in the case LL.
A summary of the configuration of these experiments is presented in Table 1.

Evaluation of simulated coastal flooding
The coastal flooding area simulated by the Ludington-LR model shows a close agreement with the reported flooding area. While no news reports or observed flood extent maps exist, videos of the meteotsunami event that were posted on social media were able to capture major flooded areas. Figure 7a shows the model simulated maximum flood extent around the Ludington north breakwater entrance, which is highly consistent with video snapshots of the flood water by local photographers (Fig. 7b-d). Both model results and recorded video snapshots show that the southern beach of Ludington Stearns Park was severely flooded with rising lake level (Fig. 7c) and the flood waters inundated the streets near West Ludington Avenue (Fig. 7d). This corroborated with an account from a local photographer that reported, "Water was also flooding the beach and the end of Ludington Avenue" (https:// www. mlive. com/ news/ 2018/ 04/ lake_ michi gan_ pier_ compl etely. html). The model simulation also successfully captured the flood waters that intruded around the local highland to the lighthouse pier entrance west of Ludington Skate Park (Fig. 7b). Furthermore, the model predicted flooding areas correspond quite well with the coastal inundation map generated from the NOAA Lake Level Viewer under the condition when the water level rises to 178.2 m (Fig. 5d). In fact, the model predicted maximum water level rise was around 178.1-178.3 m across this region (e.g., Figs. 9 and 10). The meteotsunami waves during this event were observed with wave periods between 18 and 24 min based on the observed water level fluctuations at nearby coastal gauge stations (Anderson and Mann 2020), which led to rapid water level changes that caused the Ludington breakwaters to be submerged and re-emerged in a short time period of 10 min. The most widely reported information about the Ludington meteotsunami on April 13, 2018 is a set of two photos of the Ludington North Breakwater (Fig. 8a, b), which show how the rise in lake level from the meteotsunami wave completely covered the north breakwater and retreated below the structure again approximately 10 min later. The same phenomenon was captured by the Ludington-HR model (Fig. 8c, d), demonstrating the model's ability to simulate overtopping on coastal structures in both magnitude and phase.
The Ludington-HR resolves local morphological features such as the Ludington Harbor, the river channel, the connecting Pere Marquette Lake, and the drowned river mouth estuary where the gauge station is located. Hence, the model was able to reproduce the realistic  Water level change (contour) and water fluxes per unit length (arrows) in the cases of BW0 (left panels), control run (middle panels), and BW1(right panels). Note that case BW0 has no breakwaters, and the original breakwater locations were marked in red in left panels. Water flux per unit length is calculated as vertically averaged velocity times local water depth (m 2 /s) differences in water level fluctuation between the gauged station, which is far inside the harbor, and the water level outside of the harbor in Lake Michigan. Outside of the Ludington Harbor, the model simulated a peak water level of 178.1 m (Fig. 9, blue line), which is considerably higher than the original LMHOFS predicted high water level of 177.3 m (Fig. 5a). Meanwhile, the model reproduces the much smaller water level fluctuation at the gauge location, where the peak water level only reached 177.2 m as a result of energy dissipation in the shallow estuary. This again highlights the limitations and potential mischaracterization of the meteotsunami hazard if an early warning system fails to adequately resolve the coastal topography and nearshore dynamics.

Impact of breakwaters
The Ludington Harbor is protected by two converging breakwaters (Fig. 2b). The north and south breakwaters have a length of roughly 549 m and 518 m respectively, and create a harbor entrance that is approximately 168 m wide (https:// www. lre. usace. army. mil/ Missi ons/ Operations/Ludington-Harbor-MI/). These two breakwaters were constructed to protect the wind-induced waves (see Sect. 3.4 and Fig. 14) and sediment deposition into the harbor and to maintain the depth and width of the river channel for boating. The meteotsunami waves caused the water level to rise in 10 min on a regional scale of 50-80 km along the east coast of Lake Michigan. As meteotsunami wave periods are much longer than those of the wind-induced waves, the breakwaters provide limited protection from meteotsunamiinduced coastal inundation (Fig. 10), and in some cases, enclosed harbors can even serve to amplify the meteotsunami wave height. A barotropic pressure gradient force due to a water level difference drove the water to deluge the Ludington Harbor through the harbor entrance and through overtopping of the breakwaters. Case BW0, which is configured the same as the control run only with the breakwaters removed, was designed to examine the impact of the existence of the breakwaters on the meteotsunami-induced coastal inundation. We focused on the comparison of water transport, water level rise and flooding when the largest meteotsunami wave (the second wave) hit Ludington that resulted in flooding and submerged breakwaters between 16:10 and 16:40 (GMT) with and without the breakwaters (Fig. 10). While water levels started to rise at 16:10 (GMT) in all cases (Fig. 11a) during the 16:10-16:20 (GMT), the impact of breakwaters on slowing down the water rise inside the harbor is noticeable (Fig. 10a 1,2,3 ). With the breakwaters present, the water level inside the harbor was between 177.1 and 177.2 m, while outside the harbor the water level reached 177.4 at 16:16 (GMT) and flooding occurred along the south beach. The water level further increased to 177.9-178.0 m outside the harbor at 16:18 (GMT), but overtopping had not yet occurred in the cases with breakwaters included. The water level inside the harbor also continued increasing to 177.6-177.7 m and began to flood the nearshore streets (Fig. 10b 1,2,3 ). At this point in time, without the breakwaters (BW0), the water level inside the harbor would be 10-20 cm higher and would be similar to levels outside the harbor (Fig. 10c 1,2,3 ). At 16:20 (GMT), the water level outside the harbor reached its maximum of 178.1-178.3 m, resulting in a strong water level gradient that increased transport into the harbor and caused overtopping in the control run, reaching a peak level of 178.1-178.2 m for the inner harbor (Fig. 10c 1,2 ). In both the BW0 and control cases, flooding extent reached W. Ludington Ave, though the breakwaters had the effect of inducing a strong transport in the northeast corner of the Fig. 11 Time series of averaged water level inside the Ludington Harbor in the control run, BW0 and BW1 run (a); time series of water transport through the opening between the breakwaters, overtopping, and through the channel (the connection waterway through which water flows into and out of Pere Market Lake) in the control run (b) and BW1 run (c) and BW0 run (d). Positive values represent water transport into Ludington Harbor 1 3 harbor, which set up a higher water level near the breakwater and slightly increased flooding ( Fig. 10d 1,2 ).
The vertically integrated water transport is the driving mechanism responsible for water level change over a region. Results from case BW1 show that increasing the breakwater crest elevation to 3 m above the low water datum of 176 m would significantly reduce overtopping (Fig. 10c 3 and green lines in Fig. 11b, c) yet enhance the water transport through the opening due to the higher water setup outside the harbor. As a result, net water transport into the harbor is similar in the control run and BW1 case (black lines in Fig. 11b,  c), and the water level rise and flooding extent is only slightly reduced in the BW1 case (Figs. 11a, 10d 2,3 ). On closer look, the water level started to rise at 16:08 (GMT) with the largest meteotsunami wave, which came after two prior, smaller meteotsunami waves. Consequently, the net transport (i.e., the water transport into the Harbor through the breakwater opening and overtopping minus the transport out of the harbor into Pere Marquette Lake through the channel) to Ludington Harbor turned to positive and the water level in the harbor continued to rise, reaching a peak water level at 16:22 (GMT). Thereafter, the water level decreased to its low level at 16:36 (GMT). In the control run, the net water transport was between 16:08 and 16:22 is 8.5027e + 05 m 3 (including 2.0913e + 05 m 3 overtopping, accounting for 24.5% of net transport). In the case BW1, while the overtopping is significantly reduced, the net water transport through the opening is 7.5475e + 05 m 3 (including 0.2740e + 05 m 3 overtopping, accounting for 3.6% of net transport), which is 11% smaller than that in control run. Similarly, the general patterns of net transport were similar (black lines in Fig. 11b-d) in the cases with and without the breakwaters.

Impact of lake level
Water levels in the Great Lakes have been characterized by significant fluctuations on decadal, interannual and seasonal scales. The primary drivers of water levels in the Great Lakes are runoff, over-lake precipitation, and evaporation; collectively called the net basin supply (NBS). The Lake Michigan level declined from a relatively higher water level of 177.2 m in 1997 to a level below the long-term mean in 1999 and remained low until 2014 (as low as 175.57 m in 2013). However, the water level has increased rapidly since then. Over just six years, the lake water level has risen by ~ 2 m (Fig. 12a). In 2020, water level broke the monthly record high from January through August (Fig. 12b). Recent studies suggest that the water level rise was caused by the combination of increased precipitation and decreased lake evaporation since 2013-2014 (Gronewold et al. 2021). Regional climate projections suggest the trend of rising water may continue into the future (Notaro et al. 2015;Kayastha et al. 2021).
Case LL was designed as a "likely scenario" to examine the vulnerability of the region to natural lake level variation. The LL case assumes the same meteotsunami to occur at a high water level of 177.6 m, which is 0.8 m higher than the water level in April 2018, when the meteotsunami occurred. Such a water level is roughly 15 cm higher than the highest monthly mean water level in 2020, and equivalent to a high water level observed in southeast Lake Michigan near Chikaming Township on September 17, 2020. Results in case LL show that the flooding would exacerbate significantly if the event were occurring at a higher mean lake water level. Not only the W. Ludington Ave, but several streets to the south would also suffer from severe flooding (Fig. 13). In addition, the beach protected by the south breakwaters would also be flooded. Compared to the control run, the increased lake level in case LL changes the water transport pathways. The majority of the water transport would be through both overtopping and the breakwaters' opening. In fact, transport through the opening would decrease by ~ 41%, while that through overtopping would increase by ~ 48%.

Impact of storminess
In order to assess coastal flooding potential under possible extreme conditions in the context of climate change, we designed a "representative storm" scenario that considers the combined threats of meteotsunami and wind waves during a high lake level. In coastal regions, wave-current interactions are likely to have significant impacts on coastal flooding during storm events. Waves not only enhance coastal flooding through wave runup and overtopping, but also increase mean water level due to wave setup (Olabarrieta et al. 2011). In Lake Michigan, climate change in the past decades has caused more severe storms and rapidly rising water levels (Wuebbles et al. 2021), which allow large waves to attack the shore directly and pose greater threats to coastal systems. During the April 2018 meteotsunami event, a northwest wind was dominant over the entire lake, which produced small Fig. 12 Observed monthly mean lake wide average water level of Lake Michigan-Huron: long-term series (a) and monthly water levels (b) 1 3 Fig. 13 The water flood map simulated in the control run (top), LL case (middle), and the difference in the flood area (bottom) at peak water level timestamp 1 3 waves with significant wave height of 0.3-0.4 m (not shown) in Ludington due to the short fetch. As a result, the wave-current coupled simulation revealed that wind waves had negligible impacts during the event. However, the "representative storm" scenario assumes the same atmospheric pressure forcing as in the control run but in conjunction with a southwest wind with a peak speed of 20 m/s prevailing over the entire lake. This is the dominant wind direction during April and May, when meteotsunamis are most likely to occur in Lake Michigan (Bechle et al. 2016).
To study the wave impacts on coastal flooding, we conducted two simulations without (RS1) and with (RS2) wave-current coupling. Figure 14 shows the simulated wave height distributions in the Ludington Harbor as well as its neighboring coastal regions during the meteotsunami event. The 20 m/s southwest wind generates extremely large waves in eastern Lake Michigan with significant wave heights over 4 m. At the entrance of the Ludington Harbor, wave heights nearly reach 3 m. Inside the harbor, the wave height is significantly reduced to ~ 1.0 m near the shore because of the breakwaters, which protect the harbor from direct wind-wave attacks.
The wave impacts on coastal flooding are revealed from the difference in simulated water depths from a decoupled hydrodynamic-wave simulation (RS1) and a coupled wavecurrent simulation (RS2). The wave setup and setdown are indicated by the difference of the water depths from two simulations (Fig. 15). Overall, the waves caused an increase in water depth in the nearshore, resulting in more severe coastal flooding. The highest wave setup appeared around southern beach at Sterns Park, where the water depth increase could be higher than 0.1 m. Inside the harbor, due to smaller wave heights, the wave setup was generally lower. The water depth at the flooded streets slightly increased with wave impacts.
In addition to wave setup, wave runup on beaches could further increase coastal flooding. The upper limit of the runup, which is the maximum elevation of uprush above the still water level, determines the active beach profile and the inundated area. Prediction of wave runup requires a phase-resolving model (Ma et al. 2014) that can simulate nonlinear wave transformation and breaking and is capable of capturing a moving shoreline. The SWAN wave model employed in the current study was not aimed to predict wave runup. However, wave runup on plane, impermeable beaches could be estimated by the predictive equations proposed by Mase (1989), who found that the maximum runup ( R m ) was a function of surf similarity parameter, given by R m ∕H o = 2.32 o 0.77 , where H o is the deep water significant wave height, and o is the surf similarity parameter calculated using the deep water wave parameters. For instance, on the northern beach outside of the Ludington Harbor, the deep water significant wave height was predicted to be about 2.8 m, peak wave period was 7.4 s, and the offshore water depth was about 9.0 m during the "representative storm" event. Given the beach slope of 0.03, the maximum wave runup was estimated to be 1.63 m, which could further exacerbate the inundation.

Summary
In this study, we investigate the essential processes that contribute to the coastal hazards associated with meteotsunami-induced flooding. During a meteotsunami event, the rapid change in water level at the shoreline can lead to damage of coastal infrastructure, deluge of nearby property, and induction of dangerous currents in the nearshore. Most of the existing real-time coastal forecast systems in the world do not resolve these components, Fig. 14 Simulated significant wave heights in the east coast of Lake Michigan (a) and a zoom-in view of significant wave heights and wave directions in Ludington Harbor and surrounding region (b) in the RS2 run and thus the true scale of the hazards associated with meteotsunamis are not represented. Therefore, as meteotsunami forecasting is still in its early development, or non-existent in most regions of the world, it is critical to understand the modeling requirements for resolving these processes to develop robust forecast systems (Angove et al. 2021;Vilibić et al. 2021).
The Great Lakes are a known hot-spot of meteotsunami activity (Bechle et al. 2016;Vilibić et al. 2021), yet there is currently no available forecast system or detection system in place for meteotsunami conditions. As such, the Great Lakes serve as an important region for investigation of essential processes in meteotsunami flooding and forecast requirements. For this work, we evaluate four key hydrodynamic model attributes: (1) coastal model grid resolution and wetting and drying process in low-lying zones, (2) coastal infrastructure, including breakwaters and associated submerging and overtopping processes, (3) annual/ seasonal (ambient) water level change, and (4) wind wave-current coupling. A series of sensitivity analyses are carried out to evaluate the importance of these attributes to meteotsunami modeling, including a "representative storm" scenario in which a meteotsunami wave is generated under high ambient lake-level conditions with a preferable wind direction and speed for wind-wave growth. These results are placed in context with an existing forecast model used for hydrodynamic prediction in Lake Michigan.
In comparison to the existing Lake Michigan-Huron Operational Forecast System (LMHOFS), a high-resolution version of the model grid that incorporates breakwaters and shoreline topography (Ludington-HR) yields an improved water level simulation as measured against observations at the nearby Ludington water level gauge. Furthermore, the Ludington-HR model resolves coastal flooding into nearby beaches and city streets that agrees with eyewitness reports as well as photographs and video taken during the event. Although the harbor breakwaters add increased numerical complexity, our results show that they modulate the wave amplitude and exacerbate coastal flooding, particularly when overtopping occurs. In fact, overtopping itself is often a primary contributor to meteotsunami fatalities, and thus is critical to hazard assessment. The results of the Ludington-HR control case show breakwater overtopping that is corroborated by photographs of the event (Fig. 8a). While ocean coasts face sea level rise, the Great Lakes undergo large scale interannual and seasonal lake level fluctuations on the order of 1-2 m. In contrast to tidal fluctuations, the persistence of ambient lake level conditions can impact an entire meteotsunami season (April-July) or multiple years. Connections between ambient lake level and other hazards like coastal erosion have been documented; however, the link between lake level and meteotsunami impact has not been previously explored. The results in the LL case show that high ambient lake-level conditions exacerbate flooding extent into the coast and breakwater overtopping. The increased overtopping of the harbor breakwaters also yields a significant shift in water transport pathways into the harbor, where flow over the breakwaters increased by 48% and flow through the harbor entrance decreased by 41%, as compared to the control case.
The impact of wave-current interaction in the nearshore is critical to characterizing beach hazards and coastal flooding. Previous studies have demonstrated the link between wave conditions and dangerous currents during a meteotsunami event (Linares, et al. 2019). Here, we show how wave-current interaction impacts coastal flooding using results from a "representative storm" scenario that has favorable wind direction and speed for windwave growth during a meteotsunami event. The coupled wave-current simulation reveals that while flooding extent into the nearshore is only slightly greater in the coupled case, the intensity or depth of flooded waters increases by up to 0.2 m, which can be important to property damage and material transport in flooded areas. In addition, wave runup can further exacerbate the inundation.
Finally, we note that these four key model attributes we discussed above are from the perspective of lake hydrodynamic processes and modeling. Another critical factor in predicting meteotsunami-induced flooding is the accuracy of meteorological forcing, which comes to play from the atmospheric perspective. In this meteotsunami event, the adjustment of atmospheric pressure plays a fundamental role in improving the simulation of the magnitude of meteotsunami waves. In this study, our focus is to evaluate essential hydrodynamic processes and identify key hydrodynamic model attributes. However, it must be noted that the effort must also be dedicated to improving the atmosphere modeling accuracy, in addition to advancing hydrodynamic model attributes and features, to enhance the capability of real-time forecasting systems.
Overall, this study uses a numerical modeling approach to evaluate physical processes that are essential to characterizing shoreline meteotsunami impacts. These results demonstrate that accurate representation of coastal infrastructure and topography, ambient water level conditions, and wind-wave conditions can be critical to forecasting the realistic hazards associated with meteotsunami events. This work shows that calibrating or assessing existing forecast models against coastal water level gauges alone may result in underestimating the meteotsunami hazard, as gauged levels can represent dampened or otherwise unrepresentative pictures of meteotsunami intensity, particularly when gauging stations are sparse and located behind harbor breakwaters or inside estuaries. Creating such highresolution modeling systems for real-time applications with all these hydrodynamic model attributes, particularly along the entire coastline, requires a massive amount of computational resources that may not be feasible in the near-term model development. A potential alternative is to identify meteotsunami-prone locations for the implementation of the before-mentioned high-resolution model system. While existing coastal forecast systems generally do not resolve many of these features, this work serves as a template to guide where resources may be most beneficial in model development and implementation in concert with the relentless growth in computational power and fast evolution in earth system models.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.