Time-varying storm surges on Lorentz’s Wadden Sea networks

We present an idealized network model for storm surges in the Wadden Sea, specifically including a time-dependent wind forcing (wind speed and direction). This extends the classical work by H.A. Lorentz who only considered the equilibrium response to a steady wind forcing. The solutions obtained in the frequency domain for the linearized shallow-water equations in a channel are combined in an algebraic system for the network. The velocity scale that is used for the linearized friction coefficient is determined iteratively. The hindcast of the storm surge of 5 December 2013 produces credible time-varying results. The effects of storm and basin parameters on the peak surge elevation are the subject of a sensitivity analysis. The formulation in the frequency domain reveals which modes in the external forcing lead to the largest surge response at coastal stations. There appears to be a minimum storm duration, of about 3–4 h, that is required for a surge to attain its maximum elevation. The influence of the water levels at the North Sea inlets on the Wadden Sea surges decreases towards the shore. In contrast, the wind shearing generates its largest response near the shore, where the fetch length is at its maximum.


Introduction
Storm surges, the raised water levels induced by strong winds in coastal areas, pose a serious hazard of flooding and of loss of life and property. This is amplified by trends such as a growing population pressure, sea level rise, and increasing storminess projections due to climate change (Pugh 1987). Weather systems act upon the water by means of the atmospheric-pressure differences and the wind stresses acting on the free surface. The surge dynamics is further influenced by tides, non-linear tide-surge interactions, wave dynamics, bed interactions, and the physiographic features of the coastal area-see, e.g., Pugh (1987) for general information on storm surges.
Storm surges have their greatest impact on shallow seas, in embayments, and on shores of low-lying lands. A region combining these vulnerabilities in the Netherlands was the Zuiderzee (lit. Southern Sea), a fringe basin in the southern North Sea, which in the nineteenth and early twentieth centuries had been affected by 18 major floods (Rijkswaterstaat 1916). In the aftermath of the deadly flood of 13 January 1916, the political consensus was reached to separate permanently the Zuiderzee from the outer basin with a dam, which was completed in 1932. The placemarks DO and KZ of Fig. 1a indicate the 32-km-long Afsluitdijk (lit. closure dike) dividing the former Zuiderzee bay (the present-day lake IJssel) from the basin between the mainland and the tidal islands (the Wadden Sea). In preparation of the construction works, in 1918, a task panel, the State Committee on the Zuiderzee, was appointed to determine the change of peak elevations caused in the Wadden Sea by the diversion of the tidal and wind-driven currents. Its chairman, the Nobel-laureate H.A. Lorentz, asserted the necessity of a novel investigation based on first principles (e.g., Mazure 1963 andKox 2007).
The report of State Committee for the Zuiderzee (1926) provides a pioneering idealization of the physiographic complexity of a tidal basin. The extensive flats of the Wadden Sea are separated by deep channels originating at the inlets, also visible in Fig. 1a. Since most of the water The "tidal network" used by the State Committee for the Zuiderzee for simulating tidal flows (their Section 45). c The "storm network" used for simulating storm surges (their Section 89), with channels representing the Zuiderzee in blue (see Section 3.2.1) flows along these tidal channels, the State Committee for the Zuiderzee represented the western Wadden Sea as a network of equivalent channels having depth and width uniform over their length-specifically commented upon in Section 2.1. Whereas originally this simplified approach aimed to overcome limitations of computing, idealized basins have been used until recently for investigating phenomena in coastal dynamics, such as tides (e.g., Hill and Souza 2006;Alebregtse et al. 2013;Alebregtse and de Swart 2014;Alebregtse and de Swart 2016), storm surges (Stroband and Wijngaarde 1977), and tide-surge interactions (Prandle and Wolf 1978). Their enduring advantage over the models retaining the full complexity of physics and topographyfor example Zijl et al. (2013), Duran-Matute et al. (2014), and Duran-Matute et al. (2016)-lies in their computational efficiency. Idealized models do reproduce key physical processes at a limited computational cost, and provide accurate results (both in a quantitative and qualitative sense), for example usable for extensive sensitivity analyses against geometrical and physical modelling parameters driving the system's response-as presented in this article.
Further, in order to compute the one-dimensional flow inside the individual channels, the State Committee for the Zuiderzee (1926) linearized the shallow-water equations averaged over the channels' cross sections. To this end, the seabed friction was parametrized through a novel procedure based on energy-equivalence arguments, known as Lorentz linearization, described in Section 2.2.2.
Unfortunately, for the lack of adequate processing power at the beginning of the twentieth century, the tide and storm-surge simulations could only be performed separately. In particular, the storm-surge simulations consisted of calculating the steady-state water levels in equilibrium with an extreme wind having fixed speed and direction. Therefore, the State Committee for the Zuiderzee could not identify that both motion and storage of the surge water inside the Wadden Sea are modulated by the temporal variability of the wind field, as well as by the fluxes across the tidal inlets-see, for example, the reanalyses of Lipari et al. (2008), van Vledder (2009), andDuran-Matute et al. (2016). This implies that, in a semi-enclosed basin, the most severe surges are not necessarily generated by the storms with the highest wind speed. Time-varying storms can cause higher surges than steady-state storms do, even when the peak wind speed is the same (Lipari et al. 2008). Jallah and Bakker (1994) already coded a computerized transcription of the model and algorithms of the State Committee for the Zuiderzee. Here, inspired by Lorentz's seminal studies of the Wadden Sea surges, and drawing from the work of Chen et al. (2015Chen et al. ( , 2016, we have developed a new idealized network model allowing for time-varying external forcing, namely the storm surge elevation in the outer sea and the wind-stress field. The model is based on the linearized shallow-water equations applied to the flow in the network channels. Unlike the State Committee for the Zuiderzee, the equations are cast in the frequency domain after Fourier transformation of both input and output variables and the equations. Because of the linearity, the superposition of the solutions of the individual modes gives the unsteady solution in the time domain.
With this tool, we aim to provide insights on the transient behavior of storm surges within tidal basins. We will specifically investigate how basin characteristics, storm characteristics, and different forcing mechanisms affect the transient behavior of storm surges in tidal basins.
The model, the forcings, and the outline of the solution method are presented in Section 2. The simulations presented in Section 3 deal with the hindcast of the storm surge of 5 December 2013, the sensitivity analysis of the peak surge on the coast to selected basin and storm parameters, and a modal analysis of the surge response to the external forcing. Finally, Sections 4 and 5 contain the discussion and the conclusions.

Network geometry
The study area covers the western sub-basin of the Wadden Sea from the Texel Inlet in the west to the Frisian Inlet in the east, as shown in Fig. 1a. The State Committee for the Zuiderzee schematized the routing of water along the tidal channels with two networks of rectangular channels, with each individual rectangular channel having uniform depth (h) and width (b) over their length (l) (Fig. 1b,c). The finer "tidal network" (Fig. 1b) was used for the sole simulations of tidal dynamics. Figure 1c shows the coarser "storm network" conceived for the simulation of the storm surge. We have borrowed both networks for our analyses in Section 3. Both networks of channels consist of a number of nodes with channels between them. It is at the nodes that boundary conditions are imposed on the channels (see Section 2.2.3). This way the channels are linked together and the influence of the open sea or coast is accounted for. The characteristics of the two networks are summarized in Table 1; for the exact data, see Section 45 (tidal network) and Section 89 (storm network) in the report of the State Committee for the Zuiderzee (1926). A key characteristic of channels is that their width (b) is much larger than their depth h, i.e., h/b 1 (maximum value in our networks h/b = 0.03) and that the cross-channel variation in width is small.

Governing equations
The hydrodynamic model simulating unsteady winddriven flows in a network consists of a system of onedimensional, cross-sectionally averaged (since h b and the cross-channel variation is small, see Section 2.1), linearized shallow-water equations written for each j -indexed channel: where t is time, x denotes the position on the channel axis, h j is the constant channel depth with respect to the undisturbed water level, ζ j (x, t) is the corresponding free surface elevation, and u j (x, t) is the cross-sectionally averaged flow velocity. Further, τ lin b,j (x, t) is the linearized bed shear stress, further specified and discussed in Section 2.2.2, while τ w (t) andθ j (t) are the wind shear stress and the angle between the wind direction and the positive direction of the channel axis, both further specified and discussed in Section 2.2.4. Finally, g = 9.81 m s -2 is the gravitational acceleration, and ρ = 10 3 kg m -3 is the density of water. Henceforth, we will refer to the cross-sectionally averaged velocity as velocity, and to the free-surface elevation with respect to the undisturbed water level as elevation.

Lorentz's linearization of the bottom friction
For the tidal simulations of the State Committee for the Zuiderzee, Lorentz proposed a linearization of the bottom shear stress τ b . Rather than using a quadratic formulation, such as where χ j is the Chézy smoothness coefficient, the parametrization with a linearized friction coefficient r j paved the way to a closed-form solution of the shallowwater equations. Lorentz required the friction coefficient r j to be such that, over the tidal cycle of the M 2 constituent, in each channel, the quadratic and linear stresses yield the same energy dissipation, hence setting with tidal period T = 2π/ω M 2 and ω M 2 the angular velocity of the monochromatic tide. For a harmonic signal for the velocity, u j = U j cos(ω M 2 t), the equality (5) leads to the expression for the friction coefficient with the amplitude of the tidal velocity, U j , providing a natural scale for it. For simulating aperiodic storm surges, no energy-based argument can be applied in a straightforward manner. To circumvent this, the State Committee for the Zuiderzee adopted a steady equilibrium approach, that allowed them to retain the quadratic friction formulation for a single moment in time. In contrast, we have implemented a linearized friction parametrization using the peak velocity attained in each channel in the simulated time, u j,peak , as velocity scale; this implies a mild form of non-linearity, discussed in Section 2.3. Our linearized friction coefficient for the storm surges reads: where a Manning formulation captures the explicit sensitivity of friction to the channel depth (μ is the Manning roughness coefficient). The corresponding linearized and quadratic friction parametrizations are then equal when the velocity in the channel is at its peak value. Further, the time-invariant part of our solution (Section 2.3) resembles the approach of the State Committee for the Zuiderzee closely.
In their surge simulations, the State Committee for the Zuiderzee used a constant Chézy coefficient χ = 50 m 1/2 s −1 , whereas we use a constant Manning coefficient μ = 0.0242 s m −1/3 . This corresponds to coastal waters with characteristic grain sizes of D 50 = 0.2 mm and D 90 = 0.5 mm with a typical depth of 5 m (Barua 2017). These conditions are typically found in tidal basins, such as in the Borndiep basin near Ameland (van Straaten 1954), and in the basin of Spiekeroog (Flemming and Ziegler 1995).

Boundary conditions
The boundary conditions for each channel are assigned at the n-indexed network nodes, and concern either the connection of the network with the open sea, or interior nodes where channels meet, or the coast. At the open sea, the time series of water elevations is imposed, say at x = 0 where the appropriate expression for the time series ζ sea (t) is specified later in Section 2.3. In the interior nodes n, where q channels meet, at all times the elevations must be equal for all channels and the sum of inflows must equal that of outflows: Here, a total of p channels cause inflows into node n and q − p channels cause outflows. Finally, the water velocity vanishes, at the coast, say at x = l:

Wind forcing
The model is forced by a spatially uniform, time-dependent wind stress of magnitude τ w (t) and direction θ(t). The wind stress τ w (N m −2 ) is represented as in Pugh (1987): with u w the wind speed at a standard 10 m height (m s −1 ), ρ air = 1.225 kg m −3 the density of air, and C w = 5.2 × 10 −4 u 0.44 w a dimensionless wind-drag coefficient, following Safaie (1984).

Outline of solution method
Drawing from the work of Chen et al. (2015Chen et al. ( , 2016, a temporal Fourier transform is applied to the governing (1) and (2), which are solved for the individual modes in the frequency domain. The superposition of the individual solutions then gives the solution to the full problem in the time domain.
The time-dependent wind stress is written as a superposition of m-indexed modes having 2M + 1 equally spaced angular frequencies ω m : with the constant complex amplitudes W m (N m −2 ). For τ w to be a real-valued quantity, we require W −m to be equal to its complex conjugate (denoted by an overbar), i.e., W −m = W m . T simul is the simulated time window, after which the transformed wind-stress signal (13) repeats itself by construction because all modes are periodic. The lowest resolved angular frequency is thus given by ω min = 2π/T simul . At the other end of the spectrum, the highest resolved angular frequency is ω max = Mω min , where the truncation number, M, determines the temporal resolution of the signal. We note that Eq. 13 includes a mode with number m = 0 that represents the time-invariant mode. For realistic simulations of a given storm event with a physically determined duration, T event , we then require T simul to be sufficiently large compared to T event to prevent that spurious sequences of returning events cause significant interferences.
Likewise, the Fourier expansions for the velocity and elevation in each channel are having space-varying complex amplitudes U j,m (x) and Z j,m (x). This same formulation holds for the boundary conditions with assigned elevations, as in Eq. 8. The full solutions of each mode are given in Appendix A.2.

Network
The solution for an entire network (see Fig. 1b,c) is obtained by solving the linear system of the channel equations and the corresponding boundary conditions at the nodes. First, for each mode, the flow solution for the network is derived. The summation of the 2M + 1 network-wide solutions then gives the evolution over the simulation period of all flow quantities at the nodes. The distribution of elevations and velocities inside each channel are obtained from summing up the expressions (20) and (21) in Appendix A.2.
Finally, the determination of the velocity scale u j,peak for the bottom friction coefficient in formula (7) implies that the above procedure is nested in a loop: starting from an estimated guess and using an under-relaxation procedure, the velocity scale (u j,peak ) is adjusted iteratively until it corresponds to the actual solution (u j ). We require that the residual R between u 2 j,peak and u 2 j should not exceed 10 −3 ; with R 2 = (u 2 j,peak − u 2 j ) and threshold R < 10 −3 .

Simulations
Our analyses are based on the simulation of two weather events: -The storm of 5 December 2013, known as Sinterklaas storm or Xavier storm (Watermanagementcentrum Nederland 2013). We have considered the measured wind-velocity and direction at the station of Vlieland, located as in Fig. 1a. The time series are shown in Fig. 2a. The wind direction during the storm was almost NW-ly. In our simulations, both the magnitude and direction are unsteady; -An artificial episode with a schematized wind stress pattern (Fig.2b), for evaluating the influence that the storm ramp-up time T ramp , storm duration T event , and wind direction θ have on the free-surface set-up and set-down at the coast. The ramp-up time T ramp is the duration of the ramp-up stage (same as ramp-down). The duration of the storm event T event is defined as the period between halfway the ramp-up and halfway the ramp-down of the storm, and should be larger or equal to the ramp-up time (i.e., T event ≥ T ramp ). The ramp-up time can be varied (in the range 0 ≤ T ramp ≤ T event ) without changing the total "amount" of wind stress experienced by the system, as illustrated by the blue shaded areas. The peak wind stress is 1.25 N m −2 , and the direction is fixed.

Hindcast
To gain confidence in our model, we first present a hindcast of the 2013 Sinterklaas storm here. We will check both the qualitative and the quantitative performance of our model, by requiring that the simulated water levels ζ modelled do not show large phase lags compared to the measured water levels ζ measured and that the maximum water level during the surge lies within the 20% error range. The storm surge of 5 December 2013 has been simulated on both networks of Fig. 1b,c by applying the measured time signals of Fig. 2a as a time-varying and spatially uniform wind field over the Frisian inlet, this means using the measurement data from the nearby measurement location. For the Eierland inlet, this means using the measurements from Vlieland Haven, and for the Borndiep inlet, this means using the data from Wierumergronden. The modelling error introduced by not having direct data for these two inlets is minor, due to the limited importance of these two inlets on the elevation deeper in the Wadden Sea (see Fig. 6).
The scatter plots in Fig. 3 show that the simulated storm surge elevations have a maximum error of around ±20% of the measurements, for all stations (columns) and both networks (rows). In most cases, the simulated surges tend to underestimate the measurements, that is the dots during the surge lie below the perfect agreement line, where ζ modelled is lower than ζ measured . The maximum errors during the surge lie, in fact, within 0.6 m with the tidal network and 0.8 m with the storm network, whereas the error in the maximum surge level is even smaller, being below 0.1 m. The line connecting the 1-h-spaced data points renders the temporal development: there, magnitude errors appear as distances from the plot bisector, and phase errors appear as loops.
Using the denser tidal network of Fig. 1b (Fig. 3, top) leads to slightly better quantitative agreement during the surge, than the coarser storm network of Fig. 1c (Fig. 3,  bottom). This is clearly observable for the surge at stations Den Oever (DO) and Kornwerderzand (KZ), since the dots in the top panel are closer to the perfect agreement line. The agreement at the station of Harlingen is possibly affected by the fact that the networks do not allow flow towards the eastern Wadden Sea. Also note that the network nodes that can be associated to each stations are not exactly the same in either network. Before and after the surge (when the water levels are lower), the storm network (Fig. 3, bottom) performs better than the tidal network (Fig. 3, top). There is a reasonably good qualitative agreement between the modelled water levels ζ modelled and the measured ones ζ measured and the error in peak elevation. Nonetheless, the performance of the model is qualitatively correct and quantitatively acceptable at all three stations, remarkably in the lack of any ad hoc calibration.
Having considered that the tidal network gives only slightly better results, at the cost of accounting for double the number of channels, we proceed to use the storm network in the remainder of our analysis. bottom row: surge network (Fig. 1c). The black bisector indicates a perfect match between model results and measurements. The red lines indicate a 20% error. Numerical parameters: T simul =10 days, M =127

Basin parameters
We study the influence of the basin on the storm surge through three parameters: sea-level rise, basin size, and friction. The results of varying these parameters on the peak of the storm surge of 5 December 2013 are shown in Fig. 4. The output stations also include a location at the dam, for the basin size parameter. Our results here are based on the fully forced model, i.e., forced by both wind stress over the domain and elevation signals at the inlets. So, we assume that the variations imposed here do not affect the elevation signals at the open boundaries. The sea level rise has been imposed by uniformly increasing the depth of the network channels. Its influence, shown in the left panel, results in an almost linear increase of elevations by the coast. Interestingly, the increase of the peak elevations by the coast is slightly smaller than the assigned sea level rise, which indicates that the combined effect cannot be reduced to the addition of surge and sea level rise.
Next, the middle panel shows the peak elevations against different southward dam displacements from its real position. Moving the location of the Afsluitdijk southwards increases both the wind's fetch length and the basin's water storage. The network has been extended by adding channels representing the Zuiderzee (as shown in Fig. 1c, in blue). A larger basin size increases the peak elevations at the fictional dam location because of the increased fetch, while reducing those at Den Oever and Kornwerderzand (at the real dam position) because of the basin's wider extent. The peak elevations at Harlingen, further away from the bay entrance, are less sensitive to the repositioning of the dam.
Finally, the right panel shows the influence of the Manning's roughness coefficient, μ, with a single value applied to the whole network. Higher roughness coefficient leads to lower elevations, since increased friction hampers the inflow of surge water into the network. The range of roughness coefficients ranges from unrealistically smooth (to identify resonance peaks); (i.e., μ = 0.006 s m −1/3  (Chow 1959). The range of the roughness coefficients used here is the same as that in the frequency response analysis of Section 3.3.1. Within the range of realistic parameters (i.e., μ = 0.018 s m −1/3 to μ = 0.042 s m −1/3 ) the impact of the uncertainty of the roughness coefficient, on the surge level at the coast is below 0.5 m. It should be noted that different estimates can be obtained when the roughness has a spatial distribution of its own.

Storm parameters
The parameters defining the artificial storm profile of Fig. 2b have been varied to highlight some influences on the peak elevations. These parameters are the wind direction θ , the peak duration T event , and the ramp-up time T ramp . Unlike in Section 3.2.1, here we have isolated the effect of wind stress from the other forcing in our model (i.e., we take the elevation signal at the inlets equal to 0). This is also a pragmatic step, since we do not know the elevation at the inlets during the artificial wind event. This approach is justified by the linearity of our model; the solution to the fully forced model is the superposition of the individual solutions of the separate forcings, except for the frictional iteration. Figure 5 shows the peak water elevations at Kornwerderzand, where we distinguish between set-up and setdown. There, a NW-ly wind (θ = 315 • ) causes the highest set-up of 75 cm, followed by the N-ly and W-ly ones (θ = 0 • , 55 cm; θ = 270 • , 55 cm). This is expected because of the downwind position of this station and because of the long fetch length in the basin. In contrast, S-ly and E-ly winds tend to push the water out of the basin towards the open sea, lowering the elevation at the coast, hence causing its set-down. This is observed for winds from the SE (θ = 135 • , -75 cm), S (θ = 180 • , -55 cm), and E (θ = 90 • , -55 cm). Panels a, g, and h also indicate that a minimum storm duration is required before the maximum peak surge is reached. This minimum duration is also observed in panels c through e for a set-down. For Kornwerderzand, this minimum duration for maximum set-up is of about 3-4 h. For locations Den Oever and Harlingen, the peak elevation (ζ peak ) is somewhat lower (55 and 45 cm, respectively) and is reached around the same time (after 3-4 h). The influence of ramp-up time is limited to minor variations in the resulting surge, at least for a fixed wind direction.

Frequency response analysis
The formulation of a time-varying process in the frequency domain makes it possible to dissect the causal relationship between external forcings and resulting flow fields into the constituting individual modes. To this end, we consider as many forcing scenarios as there are modes: in the m-th scenario, the m-th mode in the forcing has unit amplitude, while all other modes are 0. Again, as in Section 3.2.2, this is motivated by the model's linearity. The corresponding solution highlights the degree of sensitivity of the overall response to that unimodal unit forcing. Therefore, we can identify the modes in the forcings (external hydrography, wind shear) that are conducive to higher responses (elevations) at any selected location (the coastal stations). The forcing factors considered separately are the elevations assigned at the network open boundaries (the tidal inlets), and the direction of the wind stress vector, these simulations have been carried out on the storm network. Figure 6 shows the elevation response to monochromatic unit boundary conditions at each of the five tidal inlets. The wind stress is nil and the velocity scale (u j,peak ) used in Eq. 7 is fixed at 1 m/s. The response is quantified by the amplitude of the complex amplitude of the elevation, Z m . With the different Manning friction coefficients (reference value is μ = 0.0242sm −1/3 ), the elevations at the Texel and Vlie inlets generate the largest response at Den Oever and Kornwerderzand, the stations nearby. All these locations are in the westernmost part of the basin, unlike the Borndiep and Frisian Inlets that are somewhat further east and exert no influence at the output stations. The narrow Eierland Inlet, also in the west, has a smaller influence. For the Texel Inlet, the modes at angular velocities of around ω = 3 × 10 −4 rad s −1 (T = 8.7 h) determine the highest response. For the Vlie Inlet, this occurs for angular velocities around or just under the same value, depending on the station.

Forcing of elevations at the tidal inlets
Higher friction coefficients (μ = 0.033 s m −1/3 and μ = 0.042 s m −1/3 ) result in a lowering of the response for all cases. This is due to the increased friction holding back the flow of water and reducing the response at the measuring stations further in the basin. Lower friction coefficients (μ = 0.015 s m −1/3 and μ = 0.006 s m −1/3 ) lead to strong variations in the frequency response at the observation locations. There are clear response peaks around ω = 3 × 10 −4 rad s −1 ; these are observed at all three locations for forcings originating at one of three inlets (Marsdiep inlet, Eierland inlet, and Vlie inlet). In practice, this maximum response suggests that the basin experiences resonance due to waves of this precise frequency. For the two other inlets (Borndiep inlet and Frisian inlet), resonance peaks can be observed at different frequencies. We also observe that the effect of the friction coefficient is small for modes with a low angular frequency. It is at these low angular frequencies that the power of a storm is concentrated (gray bars in Fig. 6), since they are low-frequency events (i.e., non-recurring) in our simulation window. So, the effect of changes in roughness coefficient will be limited for the storm surge. Figure 7 shows the corresponding phases expressed in the frequency domain corresponding to the amplitudes in Fig. 6. The variations in phase are the largest if observation location and inlet are the furthest away (e.g., DO and Frisian Inlet) and vice versa (e.g., HL and Vlie Inlet). This indicates that distance between inlet and observation location is again important to the frequency response. A second result is that changes in friction coefficient μ result in a phase shift of the  Fig. 6, but now plotting the phase (φ) instead of the modulus of the onshore elevation response for different values of the roughness coefficient μ the phase changes faster, as can be seen in Fig. 7. The opposite effect occurs for lower friction coefficients. These result in a slower phase change overall (i.e., over the entire spectrum), although the response becomes more susceptible to fluctuations due to resonance.
The range of the roughness coefficients used here is the same as that in the sensitivity analysis for the roughness coefficient in Section 3.2.1. Figure 8 shows the elevation response at the output stations to a unimodal unit wind stress. The water elevations at the boundaries are 0. The strongest responses occur for SSE-ly (157 • ) and NNW-ly winds (337 • ) with the greatest sensitivity to forcing modes of ω = 2 × 10 −4 rad s −1 . In the Wadden Sea, NW-ly wind directions tend to cause elevation set-up, the SE-ly set-down, as seen in Fig. 5. The station most sensitive to wind influence appears to be Kornwerderzand, while the station that is the least sensitive is location Harlingen. At higher frequencies, the response is much smaller and independent from the wind direction. The frequency at which this happens is around 4 × 10 −4 rad s −1 . This corresponds to a wave period of just over 4 h. So, the effect of a fast oscillating wind stress is much smaller than that of slower oscillating wind stress. Due to the longer event time of storms (on the timescale of days), they are mainly composed of slowly oscillating wind stresses. So, for the basin to generate the largest responses, a minimum storm event duration is required. This is in agreement with the finding in Section 3.2.2. Figure 9 shows the phase response at the output stations, in a similar way to Fig. 8. For lower frequencies, the phase response mainly shows two opposite responses. Thus, the wind amplitude peaks in Fig. 8 act in an opposite direction, since the phase is opposite. This corresponds to winds from opposite directions, causing an opposite effect. As shown in Fig. 5, winds from the SE have an opposite response to those from the NW.

Discussion
The network-based idealized model for the Wadden Sea described in this study appears to be well-suited for gaining insight in the behavior of unsteady storm surges in a semienclosed tidal basin. The hindcast of the storm surge of 5 December 2013 presented in Section 3.1 agrees with measurements within 25% of magnitude and with small phase errors at the peak time, all the simplifications of its construction and settings notwithstanding. Here, we will discuss the model assumptions and the factors that affect storm surges. dynamics implies the neglect of tide-surge interactions, the relevance of which has been noted by Spencer et al. (2015) in their study of the same 2013 storm surge on the English coast: "Storm surge impacts are not simply linearly related to maximum elevation but rely on more complex, nonlinear interactions between tide-surge condition." Nonetheless, there is reasonable agreement between measured and simulated elevations for the 5 December 2013 storm surge. Furthermore, Prandle and Wolf (1978) signal that on the Thames estuary the dominant interaction mechanism between tides and storm surges is non-linear (quadratic) friction. Instead, our model implements a channel-wise, time-invariant, linearized parametrization of bottom-friction using the peak velocity as scale. Also, this friction coefficient overestimates bottom friction before and after the storm, when the actual velocities are smaller. Overcoming this limitation is already the focus of ongoing research (Roos et al. 2017).

Critique of the model assumptions
Along the same lines, Horsburgh and Wilson (2007) explained the surge clustering at the time of rising tide mathematically as the consequence of a tidal phase shift combined with the modulation of surge generation due to water depth. Another assumption of linear dynamics that water elevations are small in comparison to the water depth becomes less realistic under storm conditions.
Here, we should notice that Lorentz's networks do not allow the water to flow further into the eastern Wadden Sea, which naturally occurs and can be an important factor in the basin-wide surge dynamics, hence of which storms actually generate severe surges (Lipari et al. 2008 andDuran-Matute et al. 2016). The importance of this is not reflected in our model results: we found that the two easternmost inlets only have a limited effect on the rest of the basin (see Fig. 6).

Determinants of the surge elevations
Storm surges in tidal basins are caused by a local wind effect, a set-up of water at the tidal inlets (e.g., as an indirect wind effect), and atmospheric pressure effects (neglected in this study). Our linear model is capable of simulating storm surges in tidal basins and, by linearity, allows us to study the effect of the different forcing mechanisms separately. Both the wind and elevation in the tidal inlets show a stronger response at lower frequencies (and thus longer periods), at locations Den Oever, Kornwerderzand, and Harlingen in the Wadden Sea. The elevation at the Marsdiep and Vlie inlets (the two dominant inlets) shows peak responses around ω = 3×10 −4 rad s −1 , with lowering responses at lower frequencies (for μ = 0.0242 s m −1/3 ). The wind response shows distinct peak values at frequencies in the range of ω = 0 rad s −1 to ω = 4 × 10 −4 rad s −1 . At higher frequencies, the response is significantly lower. Waves with a long period are only present in storms with a sufficiently long duration. Therefore, storm surges in tidal basins require a minimum duration before they attain their maximum elevation.
We studied this minimum duration by focusing only on the locally generated storm surge (i.e., the local wind effect).
We found that a minimum duration in the order or several hours (3-4 h) is required for the wind stress in the domain to generate its maximum surge. The largest set-up is obtained for winds from the NW; and the largest set-down, for winds from the SE.
The minimum duration required for reaching the peak surge is caused by the non-instantaneous adaptation of the water to changes in the forcing (e.g., a sudden increase of wind stress, or an increase of water in the inlets). It is important to realize that storms in real life never have a constant wind stress, but rather a time-varying wind stress (as can be seen in Fig. 2 for the 2013 Sinterklaas storm). Therefore, during a storm, the system constantly adjusts to the evolving time-dependent forcing instead of aiming at steady-state equilibrium.
The influence of the elevations imposed at the tidal inlets on the surge at the coast can be ascribed to two factors: the inlet size and the distance from the tidal inlet, given the network. The narrower Eierland Inlet has a trailing influence on the water elevations at the output stations Den Oever and Harlingen, regardless of the proximity to it. Then, the response generated by the conditions set at the tidal inlets fades with the distance from them. For example, the Texel Inlet influences Den Oever more than elsewhere, and likewise for the Vlie Inlet and Harlingen. Furthermore, there is a small response at Den Oever, Kornwerderzand, and Harlingen to water level variations in the Borndiep Inlet and Frisian Inlet. The wind forcing shows an opposite relationship with the distance from the boundary shoreward because of the fetch length. Locations with a shorter fetch length (when considering a NW storm), such as Den Oever and Harlingen, show smaller surge response to a unit wind stress than locations with a larger fetch length like Kornwerderzand.

Conclusion
We have developed a new idealized network model for storm surges in the Wadden Sea, inspired by the model of the State Committee for the Zuiderzee (1926), and extended with a time-dependent wind stress (both in magnitude and direction). We probed the validity of our approach by simulating the 5 December 2013 Sinterklaas storm, leading to sufficient confidence in our model. The effect of basin characteristics on the transient behavior of storm surges was studied using our full model (i.e., including wind forcing over the domain and elevation signals forced at the inlets). The effect of storm characteristics on storm surges was studied using only wind forcing.
To study the effect of basin characteristics on storm surges, we investigated the impact of changes in sea level, basin extension, and friction coefficient. We found that sea level rise, expressed in deepened channels, results in an increase of the surge height, slightly less than the amount of sea-level rise imposed. The effect of a basin extension, up to the historic situation before the closure of the Zuiderzee, was found to be that the surge height increases at the back of the basin. In the middle of the basin, near the present-day location of the Afsluitdijk, we found lower elevation. The sensitivity of the model to the roughness coefficient showed that an increase in bottom roughness results in a lower surge and vice versa.
The effect of storm characteristics on the transient behavior of storm surges was studied by varying the rampup time, duration, and wind direction for an artificial wind event. We found that the effect of the ramp-up time is limited, whereas a minimum duration is required for a surge to attain its maximum level. The wind direction has a clear effect on the set-up, winds from a NW-ly direction result in the highest surge, followed by winds from a N-ly and W-ly direction. A set-down is obtained when the wind direction is from a E-ly direction to a S-ly direction, with a maximum set-down for winds from a SE-ly direction.
The effect of different forcing mechanisms on the transient behavior of storm surges was studied by dissecting the solution of our model. This has been done for both the wind forcing and the elevation signal that is forced at the tidal inlets. For the elevation signal, we found that the response in the basin is affected by the distance from the inlet, and by the size of the inlet. Larger inlets have a larger impact on the basin, and inlets have a larger impact on the part of the basin that is close to them. The wind forcing was shown to have the largest effect near the coastal boundary of our basin. There the fetch is the largest, resulting in the largest set-up of water against the coast.