Unsteady Linearisation of Bed Shear Stress for Idealised Storm Surge Modelling

Abstract

The modelling of time-varying shallow flows, such as tides and storm surges, is complicated by the nonlinear dependency of bed shear stress on flow speed. For tidal flows, Lorentz’s linearisation circumvents nonlinearity by specifying a (steady) friction coefficient r based on a tide-averaged criterion of energy equivalence. However, this approach is not suitable for phenomena with episodic and irregular forcings such as storm surges. Here, we studied the implications of applying Lorentz’s energy criterion in an instantaneous sense, so that an unsteady friction coefficient $r\left(t\right)$ adjusts to the temporal development of natural wind-driven flows. This new bed-stress parametrisation was implemented in an idealised model of a single channel, forced by time-varying signals of wind stress (acting over the entire domain) and surface elevation (at the channel mouth). The solution method combines analytical solutions of the cross-sectionally averaged linearised shallow-water equations, obtained in the frequency domain, with an iterative procedure to determine $r\left(t\right)$. Model results, compared with a reference finite-difference solution retaining the quadratic bed shear stress, show that this new approach accurately captures the qualitative and quantitative aspects of the surge dynamics (height and timing of surge peaks, sloshing, friction-induced tide-surge interaction) for both synthetic and realistic wind forcings.

1. Introduction

Bottom friction is important for shallow water flows, e.g., those associated with tides and storm surges in the marine environment, as discussed in [1]. The nonlinear dependency of bed shear stress on flow speed contributes to various phenomena such as tide-surge interactions [2,3] and the deformation of tidal curves [1]. In many natural systems, bottom friction even acts as the dominant source of nonlinearity, as is the case for the Thames Estuary, UK; see [2].

At the same time, this nonlinearity may complicate solution techniques in process-based models. This is particularly the case for the subclass of idealised or exploratory process-based models, e.g., [4], which aim to study physical phenomena in isolation. These models usually schematise processes and geometry, thus enabling one to determine the analytical solution to (parts of) the problem, which in turn allows for rapid extensive sensitivity analyses.

The nonlinearity of the bed shear stress can be circumvented by a linearisation that simplifies the solution while being physically acceptable for the problem under consideration. Assuming one-dimensional channel flow, with cross-sectionally averaged velocity u and bed shear stress ${\tau }_{\mathrm{b}}$, the linearisation can be written as

$$\frac{{\tau }_{\mathrm{b}}}{\rho }=c_{\mathrm{d}}\left|u\right|u\approx ru,$$

(1)

where $\rho$ is the water density, $c_{\mathrm{d}}$ is the dimensionless drag coefficient of the conventional quadratic stress parametrisation, here assumed to be constant, and r is the linear friction coefficient (in $\mathrm{m}{\mathrm{s}}^{-1}$).

For tidal flows, Lorentz [5,6] first proposed such a linearisation. He specified the (steady) friction coefficient r by requiring that the energy dissipated by the bottom friction per unit of bed area over a tidal cycle T, $\int {\tau }_{\mathrm{b}}\mathrm{d}t$, is equal for both parametrisations in Equation (1). Assuming a monochromatic tidal flow signal $u=Ucos\omega t$ with angular frequency $\omega =2\pi /T$ and (yet unknown) amplitude U, this gives

$$r=\frac{8c_{\mathrm{d}}U}{3\pi }(\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{t}\mathrm{i}\mathrm{d}\mathrm{a}\mathrm{l} \mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{s}).$$

(2)

The combination of Equations (1) and (2) is known as Lorentz’s linearisation, also referred to as equivalent linearisation [7]. Two caveats are in order, e.g., [8].

• Equation (2) can be applied locally, leading to an r-value valid at a single point in the model domain. However, to facilitate analytical solution, the friction coefficient r is often assumed to be spatially uniform over the model domain (or over subdomains); the energy criterion introduced above should then be applied in a spatially integrated sense. This implies that Equation (2) must be replaced with a more complex expression, involving spatial and tidal averaging; also see Section 2.
• The solution to the hydrodynamic problem, consisting of flow velocity and surface elevation as functions of space and time, obviously depends on the friction coefficient r. However, the solution also depends on the velocity amplitude U through Equation (2) and thus again on the solution itself, since U is also a flow variable. This cyclic dependency can be tackled by an iterative procedure that converges to an r-value for which U agrees with the velocity amplitude of the solution. We argue that, when allowing for $r=r\left(u\right)$, the ‘linear’ stress parametrisation in Equation (2) is in fact still nonlinear in u.

Lorentz [6] applied the linearisation in his pioneering channel network model of the Dutch Wadden Sea, and found good agreement with measurements, by the standards of his time. Linearisation has also been applied in other tide-related, idealised models, e.g., of tidal inlet systems [8], large-scale tidal basins [9] and sandbank dynamics [10]. Experimental support for Lorentz’s linearisation has been provided for tidal inlet systems [11]. The linearisation can also be applied to flows with multiple tidal constituents, e.g., with a dominant constituent for which Equation (2) continues to hold and with a weaker constituent, which turns out to feel proportionately more friction [12,13]. We conclude that the linearisation of the bed shear stress in Equation (1) has been extensively applied to tidal flows.

For the flow in storm surges, however, such a linearisation is less straightforward and to the best of our knowledge, has to date not been considered. Unlike tides, the forcing of storms is episodic and irregular (Figure 1), making it difficult to identify a single velocity scale U representative for the entire event. The steady r-value that would result from applying Lorentz’s energy criterion over the time lapse of a storm event is likely to underestimate (or overestimate) friction during the stages of strong (or weak) flow. In fact, Lorentz [6] avoided these inaccuracies with a separate approach for storm surges, retaining the quadratic bed shear stress of Equation (1) and computing the steady-state equilibrium surge caused by a chosen peak value of the wind stress. Clearly, this approach could not capture the unsteady nature of both forcing and response in a storm event.

The goal of the present study was to devise and test a linearisation of the bed shear stress with a linear friction coefficient adjusting to the temporal development of natural wind-driven flows. To do so, we adopted Equation (1) and applied Lorentz’s energy criterion in an instantaneous manner, effectively turning it into a power criterion. Our approach naturally implies an unsteady friction coefficient

$$r=r\left(t\right)(\mathrm{for} \mathrm{storm} \mathrm{surges}),$$

(3)

allowing $r\left(t\right)$ to be large at times when the flow is strong and small when it is weak. The two caveats on equivalent linearisations, which continue to apply in the unsteady setting, will be dealt with by spatial integration and an iterative solution method.

We then tested and analysed our model for the idealised case of a single channel with an open and a closed end, forced by time-varying signals of the wind stress on the surface and the water level at the open boundary. Specifically, we compared the results obtained with: (i) our new approach using an unsteady $r\left(t\right)$; (ii) the same approach but now with steady r-values; (iii) a reference finite-difference solver with the conventional quadratic parametrisation of the bed shear stress. We focused on both the qualitative and quantitative properties of these solutions, including the reproduction of tide-surge interactions.

In this paper, Section 2 presents the hydrodynamic model together with the power criterion. Then, Section 3 explains the solution procedure, carried out in the frequency domain based on [14,15]. Model results are presented in Section 4, the discussion in Section 5, and the conclusions in Section 6.

2. Model Formulation

Consider a one-dimensional rectangular channel of length l, uniform width b and uniform undisturbed depth h (Figure 2). With x representing the along-channel coordinate, the channel has an open boundary at $x=0$, the ‘mouth’, and a closed boundary at $x=l$, the ‘head’. Let $\zeta (x,t)$ and $u(x,t)$ denote the width-averaged free-surface elevation and the cross-sectionally averaged flow velocity, respectively, both defined as a function of the along-channel coordinate x and time t. Positive velocities are directed towards the channel head.

Assuming that $\left|\zeta \right|\ll h\ll b$ and adopting the linearisation posed by Equations (1) and (3), the cross-sectionally averaged mass and momentum balances in linearised form become

$$\begin{array}{cc}\hfill \frac{\partial \zeta }{\partial t}+h\frac{\partial u}{\partial x}& =0,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial t}+\frac{r\left(t\right)u}{h}& =-g\frac{\partial \zeta }{\partial x}+\frac{{\tau }_{\mathrm{w}}\left(t\right)}{\rho h}.\hfill \end{array}$$

(5)

Here, g is the gravitational acceleration, $\rho$ is the water density, and ${\tau }_{\mathrm{w}}\left(t\right)$ is the along-channel component of the wind stress vector, which is time-dependent yet spatially uniform. The wind stress signal either comes from observed wind speeds (Figure 1) or represents a synthetic storm (Figure 3). It is treated as an external forcing term, implying that we do not account for the potential feedback of the water surface on the wind drag. The spatial uniformity of ${\tau }_{\mathrm{w}}$ is justified if the size of the channel is much smaller than the spatial scales of the atmospheric forcing such as, e.g., in the Wadden Sea; see [6]. In the notation of Equation (5), we emphasised the time-dependencies of the friction coefficient and the wind stress. Our time-varying wind stress suffices to represent an episodic and irregular external forcing, so we neglect atmospheric pressure gradients (but they could be readily added). Since this study concentrates on the bed shear stress, all other sources of nonlinearity (advective terms, contribution of the free-surface elevation to the water depth in the continuity and friction terms) were ignored from the outset.

The boundary conditions

$$\zeta (0,t)=f\left(t\right),u(l,t)=0,$$

(6)

prescribe a surface elevation signal $f\left(t\right)$ at the channel mouth and zero velocity at the channel head. As the initial conditions ($t=0$), we require either still water conditions or conditions in dynamic equilibrium with a tidal elevation signal at the channel mouth.

Finally, a novel friction coefficient $r\left(t\right)$ is specified by equating the instantaneous power dissipated by both parametrisations in Equation (1), integrated over the entire channel:

$$r\left(t\right)=c_{\mathrm{d}}\frac{{\int }_0^l\left|u(x,t)\right|u{(x,t)}^2\mathrm{d}x}{{\int }_0^{l}u{(x,t)}^2\mathrm{d}x}(\mathrm{for} \mathrm{storm} \mathrm{surges}).$$

(7)

Here, having performed spatial integration, we address the first caveat on equivalent linearisations in Section 1.

Importantly, the closure relationship in Equation (7) is the instantaneous (hence unsteady) counterpart of Lorentz’s classical linearisation with steady r if the latter is applied non-locally, so including spatial integration. The corresponding linearisation à la Lorentz is thus given by

$$r=c_{\mathrm{d}}\frac{{\int }_0^T{\int }_0^l\left|u(x,t)\right|u{(x,t)}^2\mathrm{d}x\mathrm{d}t}{{\int }_0^T{\int }_0^{l}u{(x,t)}^2\mathrm{d}x\mathrm{d}t}(\mathrm{for} \mathrm{tidal} \mathrm{flows}).$$

(8)

Equation (8) reduces to Lorentz’s expression $r=\frac{8}{3\pi }{c}_{\mathrm{d}}U$ of Equation (2) for the simple case of a spatially uniform and monochromatic flow velocity signal $u(x,t)=Ucos\omega t$ with angular frequency $\omega =2\pi /T$ and velocity amplitude U.

3. Solution Procedure

3.1. Model Development

Eliminating the flow velocity from Equations (4)–(6) gives a damped wave equation for the surface elevation, supplemented with nonhomogeneous boundary conditions:

$$\frac{{\partial }^2\zeta }{\partial t^2}+\frac{r\left(t\right)}{h}\frac{\partial \zeta }{\partial t}-gh\frac{{\partial }^2\zeta }{\partial x^2}=0,\zeta (0,t)=f\left(t\right),\frac{\partial \zeta }{\partial x}(l,t)=w\left(t\right),$$

(9)

where we introduced the dimensionless wind forcing $w\left(t\right)={\tau }_{\mathrm{w}}\left(t\right)/\left(\rho gh\right)$. To arrive at Equation (9), we first subtracted the space derivative of the momentum equation (multiplied by h) from the time derivative of the mass balance, and then used the mass balance again to express the friction term in $\zeta$. Here, it is essential that both $r\left(t\right)$ and ${\tau }_{\mathrm{w}}\left(t\right)$ do not depend on x. The boundary condition at $x=l$ follows from inserting $u(l,t)=0$ in the momentum equation. The problem in Equation (9) is linear in the unknown $\zeta$ and subject to two time-dependent forcing terms: the dimensionless wind forcing $w\left(t\right)$ and the free-surface elevation at the channel mouth, $f\left(t\right)$.

Alternatively, by eliminating $\zeta$ instead of u, one could also solve the problem in terms of the flow velocity u. In either case, the other quantity ($\zeta$ or u) will follow from applying one of the model equations.

3.2. Fourier Expansion

The governing equations in Equation (9) are solved in the frequency domain, leaning on the work by Chen et al. [14,15]. To this end, we expand all time-dependent variables in a discrete Fourier series according to

$$\left(\begin{array}{c}w\left(t\right)\\ f\left(t\right)\\ r\left(t\right)\\ \zeta (x,t)\\ u(x,t)\end{array}\right)=\sum _{m=-M}^M\left(\begin{array}{c}{W}_m\\ F_m\\ R_m\\ Z_m\left(x\right)\\ U_m\left(x\right)\end{array}\right)exp\left(\mathrm{i}{\omega }_{m}t\right),{\omega }_m=m{\omega }_1,$$

(10)

with truncation number M and angular frequency ${\omega }_m$. For each mode m, we introduce complex Fourier components $(W_m,F_m,R_m,Z_m\left(x\right),U_m\left(x\right))$ of the wind forcing, boundary elevation, friction coefficient, surface elevation and flow velocity, respectively. Importantly, the last two are functions of the along-channel coordinate x; they will be expressed as sums of analytical expressions. Since all physical quantities are real-valued, it follows that $W_{-m}=\overline{W_m}$, $F_{-m}=\overline{F_m}$, etc., with an overbar denoting complex conjugation.

The discrete Fourier series in Equation (10) implies that both forcing and response are actually periodic over a time span $T_{\mathrm{simul}}=2\pi /{\omega }_1$ associated with the lowest frequency ${\omega }_1$. By construction, the solution obtained here is directly in dynamic equilibrium with the periodic forcing; the transient propagation across the domain is not simulated. Choosing $T_{\mathrm{simul}}$ sufficiently large, here in the order of 10 days (Figure 3), prevents interference among subsequent storm events and ensures that the initial conditions of Section 2 are indeed satisfied.

3.3. Iterative Procedure

The problem is now solved following an iterative procedure consisting of six steps:

• Define an initial function of the friction coefficient $r\left(t\right)=r_{\mathrm{init}}$, e.g., zero or some other constant function, and determine the associated Fourier coefficients $R_m$.
• Find the surface elevation $\zeta (x,t)$ by solving Equation (9) in Fourier space. Since $r\left(t\right)$ and $R_m$ are known, this problem is linear in the unknown complex amplitude functions $Z_m\left(x\right)$. We solve it using matrix algebra (details in Appendix A.1).
• Obtain the flow velocity $u(x,t)$ by applying Equation (5). This provides the Fourier coefficients $U_m\left(x\right)$ (see Appendix A.2).
• Substitute $u(x,t)$ in Equation (7) to obtain an updated friction coefficient $\hat{r}\left(t\right)$. The integrals therein are evaluated on an equidistant grid with N points in the x-direction, using the trapezoidal rule. This step also produces new Fourier coefficients ${\hat{R}}_m$.
• Verify whether ${\hat{R}}_m$ and $R_m$ satisfy the convergence criterion

  $$\frac{1}{T_{\mathrm{simul}}}{\int }_0^{T_{\mathrm{simul}}}{\left[\hat{r}\left(t\right)-r\left(t\right)\right]}^2\mathrm{d}t=\sum _{m=-M}^M{|{\hat{R}}_m-R_m|}^2<{ϵ}^2,$$

  (11)

  with user-defined convergence tolerance $ϵ$. If so, the instantaneous power criterion is satisfied with sufficient precision over the entire simulation time, the current solution is converged in the sense of the second caveat on equivalent linearisations in Section 1, and the procedure is completed.
• If Equation (11) is not satisfied, apply relaxation to find the next guess for the friction coefficient: $r_{\mathrm{new}}\left(t\right)=\alpha \hat{r}\left(t\right)+(1-\alpha )r\left(t\right)$, with a relaxation parameter $\alpha$ that can be tuned to accelerate convergence. Then, repeat the procedure from step 2 onward using the new function $r\left(t\right)=r_{\mathrm{new}}\left(t\right)$ and the associated Fourier coefficients $R_m=R_{\mathrm{new},m}$.

The truncation number M, simulation period $T_{\mathrm{simul}}$, number N of spatial points used in the numerical integration of Equation (7), relaxation parameter $\alpha$ and convergence tolerance $ϵ$ are the user-defined parameters of our solution method: see

Table 1. Overview of parameter values (physical, forcing-related, numerical).

Parameter	Symbol	Value	Unit
Water depth	h	8	$\mathrm{m}$
Channel length	l	100	$\mathrm{km}$
Drag coefficient (bottom friction)	$c_{\mathrm{d}}$	$2.5\times {10}^{-3}$	-
Drag coefficient (wind stress)	$c_{\mathrm{w}}$	$2\times {10}^{-3}$	-
Gravitational acceleration	g	$9.81$	$\mathrm{m}{\mathrm{s}}^{-2}$
Water density	$\rho$	1000	$\mathrm{kg}{\mathrm{m}}^{-3}$
Air density	${\rho }_{\mathrm{air}}$	$1.225$	$\mathrm{kg}{\mathrm{m}}^{-3}$
Duration of storm event	$T_{\mathrm{event}}$	24	$\mathrm{h}$
Ramp-up time (equal to ramp-down time)	$T_{\mathrm{ramp}}$	12	$\mathrm{h}$
Peak value of wind stress	${\hat{\tau }}_{\mathrm{w}}$	1	$\mathrm{N}{\mathrm{m}}^{-2}$
Tidal period	$T_{\mathrm{tide}}$	12	$\mathrm{h}$
Tidal elevation amplitude	$\hat{F}$	1	$\mathrm{m}$
Truncation number of Fourier series in Equation (10)	M	512	-
Simulation period ${}^{*}$	$T_{\mathrm{simul}}$	240	$\mathrm{h}$
Relaxation parameter	$\alpha$	$\frac{2}{3}$	-
Number of points in space to evaluate power criterion	N	240	-
Convergence tolerance	$ϵ$	${10}^{-5}$	$\mathrm{m}{\mathrm{s}}^{-1}$
Space step ${}^{†}$	$\Delta x$	$1.68$	$\mathrm{km}$
Time step ${}^{†}$	$\Delta t$	60	$\mathrm{s}$

${}^{*}$ Period of the forcing and solution in our Fourier model (Figure 3); ${}^{†}$ Used in reference finite-difference solution with quadratic bed shear stress; see Section 3.5.

and Discussion in Section 5.

3.4. Procedure for the Steady Friction Coefficient

As mentioned in Section 1, we conducted simulations with a steady friction coefficient r for the purpose of comparison. Such a steady r-value is implemented by applying steps 1–3 of the procedure in Section 3.3 once. Unlike Lorentz’s original procedure, this formulation is not solely tidal and the steady r-value can be arbitrary, not necessarily tied to an energy criterion.

3.5. Reference Finite-Difference Solution with Quadratic Bed Shear Stress

A finite-difference solution with the conventional quadratic parametrisation of the bed shear stress was used to evaluate the goodness of the linear parametrisations with r and $r\left(t\right)$. To this end, we solved Equations (4) and (5) after replacing the bottom friction term $r\left(t\right)u/h$ with $c_{\mathrm{d}}\left|u\right|u/h$, recalling Equation (1).

We used a staggered uniform grid with spacing $\Delta x$ and time step $\Delta t$, consisting of elevation and velocity points that, relative to one another, were shifted over half a step in both space and time. In view of the boundary conditions in Equation (6), the grid was chosen such that $x=0$ coincided with an elevation point, and $x=l$ with a velocity point. We then applied a leap-frog scheme—i.e., central discretisations of the derivatives in Equations (4) and (5) evaluated at staggered points in space and time—and iteratively computed the velocity in the quadratic bottom friction term. The discretisation adopted is second-order accurate in space and time.

The initial conditions prescribe still water. Importantly, these conditions actually serve as a ‘cold start’ for the simulations with tidal fluctuations at the channel head. The model spins up for a first $T_{\mathrm{simul}}$-long span (see Figure 3) to simulate the transient and produces the desired solution in dynamic equilibrium with the periodic forcing in the second $T_{\mathrm{simul}}$-long span. This treatment enables a proper comparison with the Fourier solution, since in Section 3.2 we chose $T_{\mathrm{simul}}$ sufficiently long to fulfil this purpose.

The grid spacing and time step (Table 1) are such that the numerical solution is stable and sufficiently converged. Furthermore, the shortest oscillations resolved by the finite-difference solution ($T_{\min }^{\mathrm{FD}}=2\Delta t=120\mathrm{s}$) are much shorter than those associated with the highest frequency of the Fourier expansion (10) in our new model ($T_{\min }^{\mathrm{Four}}=T_{\mathrm{simul}}/M=1687.5\mathrm{s}$). Therefore, the finite-difference solution is indeed a trustworthy benchmark for the comparisons presented in Section 4.

4. Results

4.1. Wind Forcing Only: Synthetic Storm Event

First off, we tested the new model for the case of wind forcing only (Figure 4). Our reference channel (100 km long, 8 m deep, see Table 1) is subject to the synthetic wind stress signal as defined in Figure 3, with a peak value of ${\hat{\tau }}_{\mathrm{w}}=1\mathrm{N}{\mathrm{m}}^{-2}$, a duration $T_{\mathrm{event}}=24\mathrm{h}$ and a ramp-up time $T_{\mathrm{ramp}}=12\mathrm{h}$. The total simulation time is 10 days (Figure 4a). An a priori scale analysis (Appendix B) shows that, for these parameter values, bottom friction is of first order importance.

The iterative procedure converges to an unsteady friction coefficient $r\left(t\right)$ (Figure 4b), displaying pronounced peaks during ramp-up and ramp-down, as well as decaying oscillations at the tail of the stormThese oscillations are associated with the sloshing motion evident from both the elevation at the channel head (Figure 4c) and the velocity at the channel mouth (Figure 4d).

The results with the unsteady friction coefficient show excellent agreement with the reference solution retaining the quadratic bed shear stress (shown by the grey lines in Figure 4c,d). In contrast, such agreement breaks down for the two examples with steady friction also shown: one with the maximum value $r_{\max }$ and one with the mean value $r_{\mathrm{mean}}$ (averaged over the time span of $T_{\mathrm{event}}+T_{\mathrm{ramp}}=36\mathrm{h}$ with nonzero wind stress, shaded in pink in Figure 4a). The case with strong steady friction ($r=r_{\max }$) suppresses the sloshing of the surface elevation at the channel head and produces a deformed elevation curve underestimating and delaying the peak surge. Alternatively, the milder friction case ($r=r_{\mathrm{mean}}$) produces an early overestimated peak, a single undershoot after the event and no trailing oscillations. The same shortcomings occur in the velocity signal at the channel head (Figure 4f).

The scatter plots of Figure 5 contrast the results from all variants of the linear friction coefficient with the finite-difference quadratic reference. The bed shear stress at the channel mouth is also shown alongside the surface elevation at the head and the velocity at the mouth. The steady friction coefficients produce significantly poorer signals of those quantities than the unsteady formulation, regarding magnitude and/or phasing, which is confirmed by the RMSE-values in the insets. The same degree of goodness was also found at other locations in the channel (such as the centre, Figure 5d–f). Therefore, the new approach reproduces the flow solution caused by synthetic unsteady signals of the wind stress with great accuracy. Furthermore, the comparison with the reference solution confirms the merits of the unsteady linearisation of the friction coefficient over steady linearisations.

4.2. Wind and Tide Forcing

After the purely wind-driven flows, we addressed the tide-surge interactions induced by bottom stress. To this end, we considered the same channel geometry and the same synthetic wind stress signal ${\tau }_{\mathrm{w}}\left(t\right)$ as in Section 4.1 and Table 1, plus a tidal signal $f\left(t\right)$ for the surface elevation at the channel mouth (Figure 6a). The latter is an S2-signal with an amplitude $\hat{F}$, given by

$$f\left(t\right)=\hat{F}cos\left({\omega }_{\mathrm{tide}}t\right),{\omega }_{\mathrm{tide}}=\frac{2\pi }{T_{\mathrm{tide}}}.$$

(12)

Importantly, the total simulation period is an integer multiple of the tidal period, namely $T_{\mathrm{simul}}=20T_{\mathrm{tide}}$ with $T_{\mathrm{tide}}=12\mathrm{h}$.

As anticipated in the objectives of Section 1, the unsteady friction coefficient is subject to tidal fluctuations distorted by interference with the storm event (Figure 6). Unlike the wind-only case, the friction coefficient does not decay to zero because of the background tidal flow (Figure 6b).

The surface elevation at the channel head contains friction-induced tide-surge interactions (Figure 6c), because the combined wind-tide solution (orange line) differs from the sum of the wind-only and tide-only solutions (thin pink line). We note the good agreement between our new model (orange line) and the reference finite-difference solution (grey line) over the entire simulation period. The tide-surge interaction subsequently produces negative differences, i.e., a set-down overhead, during the onset of the storm; then produces positive differences, i.e., a set-up overhead, during the fall of the storm (black line); and vanishes after the storm event as expected.

4.3. Realistic Wind Forcing

Finally, we address the modelling of realistic, unsteady and irregular signals of the wind stress. To this end, we considered the same channel as in Section 4.1 and Table 1 and applied the wind stress signal of the Xavier or Sinterklaas storm in the Wadden Sea (Figure 1). Excluding tides for simplicity, we set the surface elevation at the channel mouth to zero, i.e., $f\left(t\right)=0$ in Equations (6) and (9).

In Figure 7, we compare the solution using the unsteady friction coefficient with the reference finite-difference solution. For the latter, as explained in Section 3.5, we imposed the periodic continuation of the wind stress signal in Figure 1 and simulated a period of $2T_{\mathrm{simul}}=480\mathrm{h}$, using the first $240\mathrm{h}$ for the spin-up and the second $240\mathrm{h}$ for the model results. These show that the unsteady friction coefficient follows the natural evolution of the storm (Figure 7a); and that the agreement with the reference finite-difference solution is excellent for elevation, velocity and bed shear stress (Figure 7b–d). Therefore, the new approach is also capable of accurately handling the irregular and unsteady signals of real-world storms.

5. Discussion

Denominator of bottom friction term. The bed shear stress parametrisation obviously pertains to the numerator of the bottom friction term in Equation (5). For the denominator, we assumed that $h+\zeta \approx h$ and thus relied on small values of $\left|\zeta \right|/h$. The peak water levels in the storm surge simulations presented in Section 4 are characterised by $\left|\zeta \right|/h\approx 0.25$ and this ratio is expected to increase for surges in shallower channels. Further research will address the significance and validity of the above approximation in the friction term for channels of arbitrary depth.

Power criterion. The unsteady linearisation of the bed shear stress circumvents the problem of finding an appropriate velocity scale for situations with episodic and irregular forcings, cf. U in Equation (2). In fact, this aspect is covered by the power criterion in Equation (7) and thus anchored in that physical statement. Not only do steady friction coefficients produce inaccurate wind-driven simulations (Section 4.1), but they are also arbitrary due to the lack of a physically sound basis. For example, why choose precisely $r=r_{\max }$ [16]? Or, when opting for a time-averaged value $r=r_{\mathrm{mean}}$, over what time span should one average? The choice of this time span should be independent of the artificial numerical parameter $T_{\mathrm{simul}}$.

Channel-wide bottom friction. An essential property of our approach is that the unsteady friction coefficient $r\left(t\right)$ is spatially uniform, i.e., it is representative of the entire channel. This is guaranteed for solutions in which the flow velocities $u(x,t)$ do not vary considerably over the domain. Otherwise, the spatial averaging in Equation (7) may imply a poor $r\left(t\right)$-estimate. One can test for and resolve this limitation by splitting the channel in several subchannels, each having an individual friction coefficient $r_j\left(t\right)$. This segmentation is actually a special case of the channel network discussed below. Furthermore, a constant drag coefficient $c_{\mathrm{d}}$, which effectively turns the general nonlinear stress law into a purely quadratic one, is by no means restrictive. Equation (7) can be readily adjusted to incorporate any law for the dependency of $c_{\mathrm{d}}$ on u and h.

Interpretation of the unsteady friction coefficient. Given the $r\left(t\right)$-signals shown in Figure 4b, Figure 6b and Figure 7a, we also inspect how the channel-wide unsteady friction coefficient relates to local and instantaneous flow velocities. To this end, Figure 8 shows the scatter plots of $r\left(t\right)$ versus $u(\hat{x},t)$ at the channel mouth, the channel centre and near the channel head, with each data point representing a time instant in a simulation. By construction, the magnitude of the instantaneous bed stress is proportional to the area of the rectangle having that point and the origin as opposite corners (see pink-shaded example). The grey V-shaped lines draw the values $r=c_{\mathrm{d}}\left|u\right|$ for which, recalling the bed stress parametrisation in Equation (1), the channel-wide linear formulation coincides with the local quadratic one. Please take note that this comparison is not a model verification, which was performed in Section 4 on the basis of time series of elevation, flow velocity, and bed stress. We carried out this analysis for the wind-only, tide-only, and combined wind-tide simulations presented in Section 4.1 and Section 4.2:

• For the wind-only simulations, $r\left(t\right)$ closely follows the quadratic parametrisation, as shown by the V-shaped arrangement of the data points in the top of Figure 8;
• For the tide-only simulations, $r\left(t\right)$ follows an oscillatory path and the continuous presence of nonzero currents inside the channel prevents it from approaching zero;
• Finally, for wind-tide simulations, the path of $r\left(t\right)$-values is similar to the combined path for wind and tide forcing but is additionally deformed by the tide-surge interactions.

Solution procedure and numerical settings. Our solution method in the frequency domain enables us to combine analytical solutions (Appendix A) in a way that is analogous to the classical approach with steady friction. The numerical parameters should be chosen with care. For example, as a consequence of imposing periodic solutions, the simulation period $T_{\mathrm{simul}}$ should be large enough to prevent the interference of subsequent events (as depicted in Figure 3). Furthermore, the truncation number M should be so large that the shortest fluctuations in the forcing and the solution, $T_{\mathrm{simul}}/M$ are resolved with sufficient accuracy. Finally, the iterative procedure appeared to be equally robust for the unsteady and steady friction coefficients. In all simulations carried out, it converged to a unique function $r\left(t\right)$, independent of the relaxation parameter $\alpha$ and the initial value $r_{\mathrm{init}}$ (see Section 3.3 and Table 1). The iterative procedure accounts for the dependency of $r\left(t\right)$ on $u(x,t)$ according to Equation (7) and emphasises the nonlinearity still present in the model (recall the second caveat on equivalent linearisations in Section 1).

Extension to channel network. The application of an unsteady friction coefficient to a single channel can be readily extended to a network of channels, as performed earlier for tides and storm surges—both separately [6] and jointly [16]. To this end, one should account for the orientation of subchannels with respect to wind direction. Furthermore, one must specify appropriate matching conditions at the internal nodes where subchannels meet (continuity of elevation and discharge) and boundary conditions at the terminal nodes (open or closed). The iterative procedure presented in Section 3 then concerns the friction coefficients $r_j\left(t\right)$ of the entire network. The extension to a channel network, implemented and tested by Pitzalis [17], is beyond the scope of the present study.

6. Conclusions

We devised and tested a novel time-dependent linearisation of bed shear stress extending Lorentz’s classical linearisation to the idealised modelling of storm surges in (tidal) channels. This approach consists of two elements: (i) the formulation of an instantaneous power equivalence criterion to specify an unsteady channel-averaged friction coefficient $r\left(t\right)$; and (ii) a straightforward iterative procedure, carried out in the frequency domain, to obtain both $r\left(t\right)$ and the corresponding solution of the linearised shallow-water equations. By construction, the unsteady friction coefficient $r\left(t\right)$ adjusts to the intensity of the flow inside the domain. Qualitatively, $r\left(t\right)$ is large whenever the flow is strong and small when it is weak. As a result, we circumvent the issue of specifying a single velocity scale for the bottom friction term during episodic and irregular wind forcing.

The verification with a reference solution in the time domain retaining the quadratic bed shear stress shows excellent agreement throughout the channel in the time series of elevation, flow velocity, and bed shear. Distinctive aspects of the surge dynamics influenced by bed shear stress, such as the timing and magnitude of peaks, sloshing motions and tide-surge interactions, are well captured qualitatively and quantitatively. This conclusion holds for both synthetic and realistic wind forcings. In contrast, steady friction coefficients fail to produce such an agreement.

The unsteady linearisation of bed stress has led to a process-based idealised model that is amenable for extensive sensitivity analyses of storm surges. The model produces time series of bed shear stress throughout the channel, which indicates a potential application in morphodynamics, e.g., in terms of event-averaged bed evolution.