Modelling the Impact of Climate Change on Coastal Flooding: Implications for Coastal Structures Design

In the present work, the impact of climate change on coastal flooding is investigated through a set of interoperable models developed by the authors, following a modular modelling approach and adapting the modelling sequence to two separate objectives with respect to inundation over large-scale areas and coastal protection structures’ design. The modelling toolbox used includes a large-scale wave propagation model, a storm-induced circulation model, and an advanced nearshore wave propagation model based on the higher order Boussinesq-type equations, all of which are presented in detail. Model capabilities are validated and applications are made for projected scenarios of climate change-induced wave and storm surge events, simulating coastal flooding over the lowlying areas of a semi-enclosed bay and testing the effects of different structures on a typical sandy beach (both in northern Greece). This work is among the few in relevant literature that incorporate a fully non-linear wave model to a modelling system aimed at representing coastal flooding. Results highlight the capabilities of the presented modelling approach and set the basis for a comprehensive evaluation of the use of advanced modelling tools for the design of coastal protection and adaptation measures against future climatic pressures.


Introduction
Climate change is expected to have significant effects on the intensity and frequency of occurrence of extreme weather events, consequently affecting sea levels, circulation patterns, currents and waves in oceans and seas around the world [1][2][3][4]. Moving from the open sea to the densely populated coastal zones, more frequent storm surges and higher waves will be experienced through a number of impacts such as beach/dune erosion and inundation of low-lying areas [5,6]. Increased flooding risks are projected to have dire effects on socioeconomic aspects at regional and global levels [7][8][9][10][11], thus dictating the need for effectively designed coastal protection and adaptation measures [12,13].
Coastal flooding is attributed to the combined effect of tides, surges and waves acting over a broad range of scales in space and time. As scales change, so do the interactions and relative importance of the above physical processes, with their connection growing stronger in the nearshore. There, water levels caused by tides and storm surges, combined with wave setup and onshore wave propagation, can lead to the overtopping of coastal structures and the inundation of low-lying coastal areas. This interplay between waves and water levels dictates the characteristics of the models needed for the accurate representation of related processes, while scale issues allow the use of modular modelling approaches based on varying-complexity nesting schemes and the coupling of interoperable models.
The modelling perspective of the above has been investigated over the last couple of decades by various researchers. The storm surge inundation model presented by Hubert plications for the two objectives described above is presented. Model results are presented and discussed in Section 5, while Section 6 presents the conclusions drawn from this work.

Model Description
Following the rationale described in Section 1, the present work retains the viewpoint of a modular modelling approach in order to simulate storm-induced coastal flooding. The modelling toolbox used includes a large-scale wave propagation model (WAVE_LS), a storm-induced circulation model (SICIR) and an advanced nearshore wave propagation model (WAVE_BQ). Based on the analysis of the interplay between waves and water levels at different scales, the modelling sequence is adapted to fit two separate objectives. The first one, aiming at simulating flooding over large/regional coastal areas, uses WAVE_LS and SICIR. In this, WAVE_LS provides the components of the radiation stress tensor to SICIR which is afterwards used to simulate coastal flooding. The second one, aiming at investigating the effect of the presence of coastal structures in the above context, uses WAVE_BQ as well. In this, SICIR results are used as offshore boundary conditions for WAVE_BQ, which is used to simulate wave overtopping and transmission behind the structures, and wave runup on beaches.

The Large-Scale Wave Propagation Model (WAVE_LS)
WAVE_LS is based on the directional wave energy balance equation [39,40]:

∂E ∂t
where E(f,θ|x,y,t) is the spectral density of frequency f and direction θ at a point (x,y) and time t; and c x , c y , c θ are the x, y, θ components of the group velocity c g , respectively, according to: c x = c g sin θ c x = c g cos θ c θ = − c g c cos θ ∂c ∂x − sin θ ∂c ∂y (2) In Equation (2), c is the wave celerity. In Equation (1), D is the dissipation of wave energy expressed as: where H m is the maximum wave height according to [41]; ρ is the water density; and Q b is the probability of a wave breaking at a certain depth, expressed as (1 − Q b )/(lnQ b ) = (H rms /H m ) 2 according to [42], where H rms is the root-mean-square wave height. Wave diffraction is incorporated to the model by replacing c x , c y and c θ with C x , C y and C θ according to [40]: where κ is the wave number and δ is expressed as: The model is capable of simulating wave propagation in large coastal areas with complicated bathymetries, describing the phenomena of refraction, bottom diffraction and breaking. The numerical solution of Equation (1) is based on an implicit finite difference scheme and the model's output includes the four components of the radiation stress tensor (Sxx, Syy, Sxy = Syx), calculated by the well-known expressions valid for progressive waves, as in [43].

The Storm-Induced Circulation Model (SICIR)
The storm-induced circulation model is based on the depth-averaged wind-induced circulation equations, following [44]: where ζ is the water surface elevation above the mean water level; d is the still water depth; h is the total water depth (h = d + ζ); U, V are the depth-averaged velocity components along the xand ydirections respectively; g is the gravitational acceleration; ν h is the horizontal eddy viscosity coefficient and f c is the Coriolis coefficient. The terms τ sx , τ sy are the shear stress components at the water surface along the xand ydirections respectively, which represent the vertical boundary condition, expressed as: where k is the surface friction coefficient (in kg/m 3 , typically of the order of 10 −6 ; here we assume k = 10 −6 ÷ 3·10 −6 ), and W x , W y are the wind speed components along the xand ydirections (in m/s, at 10 m above sea level) respectively. The bed friction terms (τ bx , τ by ) are calculated based on the formulae proposed by [45]: where f b is the bottom friction factor, σ T is the standard deviation of the oscillatory horizontal velocity, and |U| = (U 2 + V 2 ) 0.5 . The horizontal eddy viscosity coefficient is expressed by the well-known Smagorinsky model, used for the representation of the damping by eddies smaller than the computational grid size, as:

∂U ∂y
where is the mixing length, approximated as equal to half the grid cell size dx [46]. Differential Equations (6)- (8) are approximated by finite difference equations according to the explicit scheme developed by [44]. Finally, regarding coastal inundation, the process is simulated using the "dry bed" boundary condition which, according to [47], can be written as the following set of pairs of conditions for any given grid point (i,j): where h cr is a terminal depth below which drying is assumed to occur (here this depth is set to h cr = 0.001 m).

The Advanced Nearshore Wave Propagation Model (WAVE_BQ)
Over the years, the classical Boussinesq equations have been extended by incorporating higher order non-linear terms, which can describe better the propagation of highly nonlinear waves in the shoaling zone. Nowadays, models based on Boussinesq-type equations (BTEs) are widely used to simulate waves transforming in the nearshore-up to the swash zone-and their interactions with various types of coastal protection. Exemplary reference is made to the recent works of [48][49][50][51][52][53], with [54] presenting a novel meshless numerical scheme for the solution of BTEs. A thorough overview of Boussinesq-type models can be found in [55].
WAVE_BQ is based on the higher order Boussinesq-type equations for breaking and nonbreaking waves as expressed in equations: where M u is defined as: and G as: where ζ is the wave surface elevation; the subscript "t" denotes differentiation with respect to time; U is the horizontal velocity vector U = (U,V); τ b = (τ bx , τ by ) is calculated from Equations (11) and (12), with the wave-current bottom friction factor calculated as in [56]; δ is the roller thickness, determined geometrically according to [57]; E is the eddy viscosity term, calculated according to [58]; and u o is the bottom velocity vector u o = (u o , v o ), with u o and v o being the instantaneous bottom velocities along the xand ydirections respectively. Following [56], wave breaking is initiated using breaking angle ϕ b = 30 • , which then gradually changes to its terminal value ϕ b = 10 • . The presented set of BTEs is accurate to the third order O(ε 2 ,εσ 2 ,σ 4 ) [59] and their numerical solution is based on the accurate higher-order numerical scheme of [60]. Coastal inundation is simulated as in SICIR (see Section 2.2, Equation (16) and [47]). The model is capable of simulating the phenomena of shoaling, refraction, breaking, diffraction, reflection and wave-structure interaction, as well as nonlinear wave-wave interaction. Regarding model capabilities in simulating the non-linear evolution of unidirectional or multidirectional wave fields in the nearshore, one can refer to [61] (see also [62] on the issue); regarding wave-structure interaction and energy transmission, one can refer to [63]. Further details on the model and its implementation to diverse coastal engineering applications can be found in [56,64,65].

Coastal Flooding
The capability of the presented advanced nearshore wave propagation model in the representation of coastal flooding is validated through the comparison with the twodimensional (cross-shore) experimental data by [38]. Roeber et al. [38] Figure 2 shows the comparison between measurements and model results for this test, as a series of snapshots of surface profile evolution. Again, as for the first test, measured and computed data are in very good agreement at all transformation stages. The model successfully captures wave shoaling over the relatively gentle slope, wave breaking on top of the reef crest, as well as the propagation of the wave bore (clearly identified by the bore front) over the reef flat.

Wave Transmission at Coastal Structures
Following the rationale presented in Section 1, coastal structures design is treated in this work as a means to countermeasure storm-induced flooding in present and future climates. With the focus set on the transmissivity of low-crested structures, the capability of the presented advanced nearshore wave propagation model in the representation of wave energy reduction behind LCS is validated through comparison with the formulation presented by Goda and Ahrens [34] for K t . Furthermore, and considering that due to the interplay between wave overtopping and wave infiltration energy reduction is expected to be lower behind impermeable structures than behind permeable ones, the former type of LCS is examined in this work so as to test the lower bound of their effectiveness.
The work of Goda and Ahrens [34] was based on the concept of the summation of wave energy transmitted over and through LCS. The processes were treated separately and two separate transmission coefficients were proposed, namely K t,over and K t,thru , respectively, with the summation concept applied for the derivation of K t using the approach of [66]. The formulation was validated using a versatile set of experimental datasets (851 tests in total), yielding a determination coefficient of r 2 = 0.865.  Figure 2 shows the comparison between measurements and model results for this test, as a series of snapshots of surface profile evolution. Again, as for the first test, measured and computed data are in very good agreement at all transformation stages. The model successfully captures wave shoaling over the relatively gentle slope, wave breaking on top of the reef crest, as well as the propagation of the wave bore (clearly identified by the bore front) over the reef flat. According to [34] K t for impermeable structures reduces to K t,over , which was estimated through empirical fitting to the design diagram of [35] and can be expressed as: where R c is the crest freeboard and a is expressed as: In Equation (22) L 0 is the deep water wavelength and B eff is the effective width of the structure, measured at still water level for emerged structures, at the level of 10% below the crest for zero freeboard structures and at the level of 20% below the crest for submerged structures [34].

Wave Transmission at Coastal Structures
Following the rationale presented in Section 1, coastal structures design is treated in this work as a means to countermeasure storm-induced flooding in present and future climates. With the focus set on the transmissivity of low-crested structures, the capability of the presented advanced nearshore wave propagation model in the representation of wave energy reduction behind LCS is validated through comparison with the formulation presented by Goda and Ahrens [34] for Kt. Furthermore, and considering that due to the interplay between wave overtopping and wave infiltration energy reduction is expected to be lower behind impermeable structures than behind permeable ones, the former type of LCS is examined in this work so as to test the lower bound of their effectiveness.  Figure 3 shows a sketch of the governing parameters involved in wave transmission at emerged LCS. These are: the incident and transmitted significant wave height, H i and H t ; the peak period T p ; the wave steepness s op = 2πH i /g(T p ) 2 ; the crest freeboard and width, R c and B; the structure height and seaward slope, h c and tanα; and the breaker parameter ξ op = tanα/(s op ) 0.5 . Model runs were performed for the series of input parameters presented in Table 1 testing combinations of two incident waves and four different crest freeboards, which resulted in relative freeboards R c /H i ranging from 0.13 to 1.00 and relative effective structure widths B eff /L 0 ranging from 0.075 to 0.150 (tanα = 0.4 for all runs). The discretization steps used in model (WAVE_BQ) runs were dx = 0.125 m in space and dt = 0.005 s in time. Figure 4 shows the comparison between K t values as estimated using Equation (21) and as resulted from model runs, for the Tests of Table 1. Data are in good agreement overall, with Test 4 representing a liming case that results in K t = 0. The model generally predicts lower transmission coefficients than Equation (21). The divergences are mainly observed for higher waves and lower relative freeboard values, a result that is expected considering the concept behind the formulation of [34] and the physics behind the phenomenon. ξop = tanα/(sop) . Model runs were performed for the series of input parameters presented in Table 1 testing combinations of two incident waves and four different crest freeboards, which resulted in relative freeboards Rc/Hi ranging from 0.13 to 1.00 and relative effective structure widths Beff/L0 ranging from 0.075 to 0.150 (tanα = 0.4 for all runs). The discretization steps used in model (WAVE_BQ) runs were dx = 0.125 m in space and dt = 0.005 s in time. Figure 4 shows the comparison between Kt values as estimated using Equation (21) and as resulted from model runs, for the Tests of Table 1. Data are in good agreement overall, with Test 4 representing a liming case that results in Kt = 0. The model generally predicts lower transmission coefficients than Equation (21). The divergences are mainly observed for higher waves and lower relative freeboard values, a result that is expected considering the concept behind the formulation of [34] and the physics behind the phenomenon.

Large-Scale Applications
The modelling approach for the simulation of flooding over large/regional coastal areas, as presented in Section 2, was applied to the area of the Bay of Thessaloniki, i.e., the northern part of the Thermaikos Gulf, located in the northwestern Aegean Sea, Greece. Thessaloniki is the second largest city in Greece, with a population of approx. 1 million people residing in its metropolitan area. The Port of Thessaloniki is also the second largest in the country, handling a total throughput of approx. 13 million tonnes per year and more than 400,000 containers/TEUs (2018 data; [67]). Two terminals are located west of the Port's 6th pier, both of significant regional importance. The AGET terminal, whose cement and cement products production and handling facilities include a jetty that extends approximately 600 m into the sea; and the Liquid Fuels Terminal, whose facilities

Large-Scale Applications
The modelling approach for the simulation of flooding over large/regional coastal areas, as presented in Section 2, was applied to the area of the Bay of Thessaloniki, i.e., the northern part of the Thermaikos Gulf, located in the northwestern Aegean Sea, Greece. Thessaloniki is the second largest city in Greece, with a population of approx. 1 million people residing in its metropolitan area. The Port of Thessaloniki is also the second largest in the country, handling a total throughput of approx. 13 million tonnes per year and more than 400,000 containers/TEUs (2018 data; [67]). Two terminals are located west of the Port's 6th pier, both of significant regional importance. The AGET terminal, whose cement and cement products production and handling facilities include a jetty that extends approximately 600 m into the sea; and the Liquid Fuels Terminal, whose facilities include and offshore jetty located at approximately 800 m from the coast. Figure 5 shows the geographic location and a satellite image of the study area; the six piers of the Port of Thessaloniki are identified, along with the location of the aforementioned facilities. The model domain was bounded to the south by the virtual East-West line connecting the Mikro Emvolo Cape to the western coast of the Bay (see Figure 6). The intermittently dry/wet area at the western part of the Bay (green dotted area in Figure 6) was modelled as a dry flat for the scenarios run in this work; this is justified by the consideration that the specific area-even when wet at highest tide-is covered by no more than a few centimetres of water (it is noted that small topographic variations over the flat do exist). The domain also included the projected final geometry of the 6th pier of the Port of Thessaloniki (grey crossed area in Figure 6), while it should also be noted that the artificial coast of the Bay of Thessaloniki (Port of Thessaloniki and waterfront eastwards of the Port up to the Mikro Emvolo Cape) was modelled as a solid boundary (i.e., no flooding allowed). This choice served the applications' computational efficiency, as the waves and storm surges expected in the study area (see Table 2 in the following) are covered by the design of the Port's piers and the waterfront's seawalls. The bathymetric data were extracted from digitized nautical charts acquired from the Hellenic Navy Hydrographic Service, while the topographic information was extracted by Digital Elevation Models of the NASA Shuttle Radar Topography Mission, acquired through AppEEARS [69,70].
The models were run for two scenarios of climate change-induced wave and storm surge events, representative for the study area, based on the results and analysis presented by [71]. The first scenario (henceforth denoted by LS1) envisaged a southern wave of significant wave height Hs = 1.58 m and peak period Tp = 4.60 s. The second scenario (LS2) envisaged the same wave combined with a storm surge of height SSH = 0.30 m (see Table 2). The discretization steps used in model runs were dx = 10.0 m in space and dt = 0.05 s in time for WAVE_LS, and dx = 10.0 m in space and dt = 0.125 s in time for SICIR. It should be noted that dx in large-scale wave propagation models is generally of the order of a wavelength L or less [39,40]. However, the choice of dx also depends on bathymetry and model domain characteristics; rapid bathymetric changes or diffraction effects would require smaller values, such as dx = L/5 [40]. In a modular modelling approach, the spatial step of the circulation model would have to be the same as previously selected due to model interoperability. Finally, the temporal step dt is controlled by the Courant number criterion Cr < 1 (Cr = cdt/dx, where c = (gh) 1/2 ). In this work, the choices for the spatial and temporal discretization of WAVE_LS and SICIR followed the above rules. The model domain was bounded to the south by the virtual East-West line connecting the Mikro Emvolo Cape to the western coast of the Bay (see Figure 6). The intermittently dry/wet area at the western part of the Bay (green dotted area in Figure 6) was modelled as a dry flat for the scenarios run in this work; this is justified by the consideration that the specific area-even when wet at highest tide-is covered by no more than a few centimetres of water (it is noted that small topographic variations over the flat do exist). The domain also included the projected final geometry of the 6th pier of the Port of Thessaloniki (grey crossed area in Figure 6), while it should also be noted that the artificial coast of the Bay of Thessaloniki (Port of Thessaloniki and waterfront eastwards of the Port up to the Mikro Emvolo Cape) was modelled as a solid boundary (i.e., no flooding allowed). This choice served the applications' computational efficiency, as the waves and storm surges expected in the study area (see Table 2 in the following) are covered by the design of the Port's piers and the waterfront's seawalls. The bathymetric data were extracted from digitized nautical charts acquired from the Hellenic Navy Hydrographic Service, while the topographic information was extracted by Digital Elevation Models of the NASA Shuttle Radar Topography Mission, acquired through AppEEARS [69,70]. The models were run for two scenarios of climate change-induced wave and storm surge events, representative for the study area, based on the results and analysis presented by [71]. The first scenario (henceforth denoted by LS1) envisaged a southern wave of significant wave height H s = 1.58 m and peak period T p = 4.60 s. The second scenario (LS2) envisaged the same wave combined with a storm surge of height SSH = 0.30 m (see Table 2). The discretization steps used in model runs were dx = 10.0 m in space and dt = 0.05 s in time for WAVE_LS, and dx = 10.0 m in space and dt = 0.125 s in time for SICIR. It should be noted that dx in large-scale wave propagation models is generally of the order of a wavelength L or less [39,40]. However, the choice of dx also depends on bathymetry and model domain characteristics; rapid bathymetric changes or diffraction effects would require smaller values, such as dx = L/5 [40]. In a modular modelling approach, the spatial step of the circulation model would have to be the same as previously selected due to model interoperability. Finally, the temporal step dt is controlled by the Courant number criterion Cr < 1 (Cr = cdt/dx, where c = (gh) 1/2 ). In this work, the choices for the spatial and temporal discretization of WAVE_LS and SICIR followed the above rules.

Applications for Coastal Structures Design Evaluation
The modelling approach for the investigation of the effect of coastal structures on flooding, as presented in Section 2, was applied to a typical sandy beach located in the East Macedonia and Thrace Region (northern Greece; see Figure 7). The low-lying coastal areas in the region are not densely populated (the coastal cities Kavala and Alexandroupoli do not compare to the size of Thessaloniki), but do host numerous touristic activities and sites of significant ecological importance (Nestos Delta, Lakes Vistonida and Ismarida wetlands). Furthermore, they are exposed to larger waves in both present and future climates, making them more susceptible to coastal flooding. Table 3 presents the layouts of the selected beach tested using the advanced nearshore wave propagation model (layouts and profile are presented in Section 5.2. The layouts include one of the unprotected beach, one of the beach protected by a breakwater, and one of the beach protected by a sea dike (henceforth denoted by L1, L2 and L3, respectively). The models were run for two scenarios of climate change-induced wave and storm surge events, representative for the study area, based on the results and analysis

Applications for Coastal Structures Design Evaluation
The modelling approach for the investigation of the effect of coastal structures on flooding, as presented in Section 2, was applied to a typical sandy beach located in the East Macedonia and Thrace Region (northern Greece; see Figure 7). The low-lying coastal areas in the region are not densely populated (the coastal cities Kavala and Alexandroupoli do not compare to the size of Thessaloniki), but do host numerous touristic activities and sites of significant ecological importance (Nestos Delta, Lakes Vistonida and Ismarida wetlands). Furthermore, they are exposed to larger waves in both present and future climates, making them more susceptible to coastal flooding. Table 3 presents the layouts of the selected beach tested using the advanced nearshore wave propagation model (layouts and profile are presented in Section 5.2). The layouts include one of the unprotected beach, one of the beach protected by a breakwater, and one of the beach protected by a sea dike (henceforth denoted by L1, L2 and L3, respectively). The models were run for two scenarios of climate change-induced wave and storm surge events, representative for the study area, based on the results and analysis presented by [71]. The first scenario (henceforth denoted by CS1) envisaged a southern wave of significant wave height H s = 5.0 m and peak period T p = 8.0 s. The second scenario (CS2) envisaged the same wave combined with a storm surge of height SSH = 0.30 m (see Table 4). The discretization steps used in model runs were dx = 0.125 m in space and dt = 0.005 s in time for WAVE_BQ (see Section 4.1 for WAVE_LS and SICIR).     Figure 8b shows the respective results for LS2.The models performed satisfactorily in both cases, resulting in smooth flooding contours that follow the modelled topography (it can be noted again that small topographic variations over the flat existed). For scenario LS1 the flooded area mainly covers parts of the aforementioned dry flat, with the flooded area approximately equal to 1.2 km 2 . For scenario LS2 the flooding extends to the low-lying coastal areas west of the port's 6th pier as well, with the flooded area approximately equal to 2.7 km 2 (increase of approximately 125%). The AGET Terminal appears to be mostly impacted in the case of LS2, while impacts on the operations of both the AGET and liquid fuels jetties are to be expected in either scenario.

Large-Scale Applications
Results are indicative of the models capabilities and the potential effects coastal flooding would have on the low-lying areas and operational facilities of coastal cities. Considering that Thessaloniki is located in a relatively protected bay where waves and storm surges are expected to be moderate (as reflected in the selected scenarios in Table  1), it is reasonable to expect more severe flooding in cities exposed to higher storm events, where the overtopping of coastal defence will result in flooding of the urban fabric as well.     Figure 8b shows the respective results for LS2. The models performed satisfactorily in both cases, resulting in smooth flooding contours that follow the modelled topography (it can be noted again that small topographic variations over the flat existed). For scenario LS1 the flooded area mainly covers parts of the aforementioned dry flat, with the flooded area approximately equal to 1.2 km 2 . For scenario LS2 the flooding extends to the low-lying coastal areas west of the port's 6th pier as well, with the flooded area approximately equal to 2.7 km 2 (increase of approximately 125%). The AGET Terminal appears to be mostly impacted in the case of LS2, while impacts on the operations of both the AGET and liquid fuels jetties are to be expected in either scenario.   Tables 3 and 4), as snapshots of surface profile evolution at times T and T + T/2. In Figure 9, coloured "x" symbols denote the landward limit of the flooding extent; this information, expressed as horizontal distance, is comparatively presented for all tests in Figure 10.

Applications for Coastal Structures Design Evaluation
As expected, the unprotected beach layout (L1) results in the highest flooding, with storm surge causing an increase in flooding extent by approximately 30% (+18 m). The LCS layout (L2) protects the beach quite satisfactorily for SC1 despite the quite low relative freeboard of the structure (see Tables 3 and 4), reducing the flooding extent by approximately 51% with respect to L1 (−31 m). However, the addition of storm surge in SC2 Results are indicative of the models capabilities and the potential effects coastal flooding would have on the low-lying areas and operational facilities of coastal cities. Considering that Thessaloniki is located in a relatively protected bay where waves and storm surges are expected to be moderate (as reflected in the selected scenarios in Table 1), it is reasonable to expect more severe flooding in cities exposed to higher storm events, where the overtopping of coastal defence will result in flooding of the urban fabric as well.  Tables 3 and 4), as snapshots of surface profile evolution at times T and T + T/2. In Figure 9, coloured "x" symbols denote the landward limit of the flooding extent; this information, expressed as horizontal distance, is comparatively presented for all tests in Figure 10. exceeding 7 m for CS1 (i.e., a reduction of approximately 89% with respect to L1 and of approximately 77% with respect to L2). For SC2, however, the flooding extent more than triples as the wave runs up the slope, reaches the crest of the dike, and propagates inland over the flat. It can be noted that the flat is covered by no more than a few centimetres of water in this test and the flooding extent is still smaller even from L2-SC1; nevertheless, this result is indicative of the function of such structures when facing combinations of increased waves and water levels. Coloured "x" symbols denote the landward limit of the flooding extent.   Flooding extent is measured from the shoreline at SWL (horizontal distance).

General Discussion
Elaborating further on the presented modelling system's performance in simulating coastal flooding, particular insights can be drawn from the following.
Model setup was based on the successful model validation for the representation of coastal flooding and wave transmission at coastal structures, as presented in Sections 3.1 and 3.2. The satisfactory agreement between model predictions, experimental data and a well-known formulation supported the implementation of the proposed modelling approach to the set of applications presented in Sections 4.1 and 4.2. These were designed in order to fit two separate modelling objectives, testing all three models of the modelling toolbox and exploring all aspects that are deemed to distinguish this work in relevant literature: (a) the incorporation of an advanced nearshore wave propagation model in the modelling sequence; (b) the accurate representation of the inundation process through the "dry bed" boundary condition; and (c) the capability to test the effect of coastal protection on coastal flooding by simulating wave overtopping and transmission at coastal structures.
The results presented in Sections 5.1 and 5.2 are indicative of the models' capabilities, while highlighting the need and importance of relevant works for coastal cities with exposed coastal facilities, as well as for other low-lying coastal areas of ecological, cultural or touristic importance. Future work can explore the limitations of the modelling approach presented by investigating three-dimensional wave-structure interactions and their impact on the identification of flooded areas; work on the same path can also extend Flooding extent is measured from the shoreline at SWL (horizontal distance).
As expected, the unprotected beach layout (L1) results in the highest flooding, with storm surge causing an increase in flooding extent by approximately 30% (+18 m). The LCS layout (L2) protects the beach quite satisfactorily for SC1 despite the quite low relative freeboard of the structure (see Tables 3 and 4), reducing the flooding extent by approximately 51% with respect to L1 (−31 m). However, the addition of storm surge in SC2 practically reduces the emerged crest height and thus cancels a large part of the protection provided by the LCS, resulting in a flooding extent that is increased by approximately 50% (+30 m). It is worth noting that flooding for L2-CS2 is quite close to that of the unprotected beach with no storm surge (L1-CS1), a useful insight for the design of LCS in a changing climate. Finally, the Dike layout (L3) results in the lowest flooding extent, not exceeding 7 m for CS1 (i.e., a reduction of approximately 89% with respect to L1 and of approximately 77% with respect to L2). For SC2, however, the flooding extent more than triples as the wave runs up the slope, reaches the crest of the dike, and propagates inland over the flat. It can be noted that the flat is covered by no more than a few centimetres of water in this test and the flooding extent is still smaller even from L2-SC1; nevertheless, this result is indicative of the function of such structures when facing combinations of increased waves and water levels.

General Discussion
Elaborating further on the presented modelling system's performance in simulating coastal flooding, particular insights can be drawn from the following.
Model setup was based on the successful model validation for the representation of coastal flooding and wave transmission at coastal structures, as presented in Sections 3.1 and 3.2. The satisfactory agreement between model predictions, experimental data and a well-known formulation supported the implementation of the proposed modelling approach to the set of applications presented in Sections 4.1 and 4.2. These were designed in order to fit two separate modelling objectives, testing all three models of the modelling toolbox and exploring all aspects that are deemed to distinguish this work in relevant literature: (a) the incorporation of an advanced nearshore wave propagation model in the modelling sequence; (b) the accurate representation of the inundation process through the "dry bed" boundary condition; and (c) the capability to test the effect of coastal protection on coastal flooding by simulating wave overtopping and transmission at coastal structures.
The results presented in Sections 5.1 and 5.2 are indicative of the models' capabilities, while highlighting the need and importance of relevant works for coastal cities with exposed coastal facilities, as well as for other low-lying coastal areas of ecological, cultural or touristic importance. Future work can explore the limitations of the modelling approach presented by investigating three-dimensional wave-structure interactions and their impact on the identification of flooded areas; work on the same path can also extend to investigate coastal flooding of the urban fabric and combine flooding extent results with damage cost models (e.g., [72,73]).

Conclusions
Coastal flooding is on the rise in many areas around the world. The risk of coastal flooding is expected to further increase in the future, as tides, surges and waves will be significantly affected by climate change and the consequent increase of extreme weather events. Numerical modelling can be used to capture the interplay between waves and water levels acting over a broad range of scales in space and time, with relevant modelling systems producing results of direct use for effective design of coastal protection and adaptation measures.
This work presents an in-house, versatile tool for the representation of coastal flooding, which follows a modular modelling approach and consists of three models: a large-scale wave propagation model, a storm-induced circulation model and an advanced nearshore wave propagation model based on the higher order Boussinesq-type equations for breaking and non-breaking waves. The numerical model at the higher-resolution end of the modelling sequence is validated through comparison with the experimental data by Roeber et al. [38] (coastal flooding) and with the formulation of Goda and Ahrens [34] (wave transmission at coastal structures), performing very well for all tests. The proposed modelling approach was applied to simulate coastal flooding: (a) over the low-lying areas of the Bay of Thessaloniki, a semi-enclosed bay located in northern Greece that hosts a major coastal city and significant economic activities; and (b) on a typical sandy beach of the East Macedonia and Thrace Region (also in northern Greece), testing two layouts of coastal protection structures in order to evaluate their effectiveness. Applications for both cases were made for projected scenarios of climate change-induced wave and storm surge events, representative of the respective study areas. Results highlighted the potential effects coastal flooding would have on the low-lying areas and operational facilities of coastal cities, as well as the function of typical coastal protection structures in the above context. This paper is among the few in relevant literature that incorporate an advanced wave model to a modelling system aimed at representing coastal flooding, along with the accurate representation of the inundation process through a straightforward boundary condition and the capability to test different layouts in order to assess the effect of coastal protection. Drawing insights from the full body of work presented here and relevant literature, this study is considered as a solid example of how the use of advanced modelling tools can facilitate strategic planning and adaptation in coastal areas, shaping the pathways towards a successful response to climate change.
Author Contributions: Conceptualization, methodology, model development, model validation, model applications, formal analysis, visualization, writing, A.G.S. and T.V.K. Both authors have read and agreed to the published version of the manuscript.