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
, the linearisation can be written as
where
is the water density,
is the dimensionless drag coefficient of the conventional quadratic stress parametrisation, here assumed to be constant, and
r is the linear friction coefficient (in
).
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,
, is equal for both parametrisations in Equation (
1). Assuming a monochromatic tidal flow signal
with angular frequency
and (yet unknown) amplitude
U, this gives
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
, 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
allowing
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 ; (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
, the ‘mouth’, and a closed boundary at
, the ‘head’. Let
and
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
and adopting the linearisation posed by Equations (
1) and (
3), the cross-sectionally averaged mass and momentum balances in linearised form become
Here,
g is the gravitational acceleration,
is the water density, and
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
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
prescribe a surface elevation signal
at the channel mouth and zero velocity at the channel head. As the initial conditions (
), 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
is specified by equating the
instantaneous power dissipated by both parametrisations in Equation (
1),
integrated over the entire channel:
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
Equation (
8) reduces to Lorentz’s expression
of Equation (
2) for the simple case of a spatially uniform and monochromatic flow velocity signal
with angular frequency
and velocity amplitude
U.
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
and thus relied on small values of
. The peak water levels in the storm surge simulations presented in
Section 4 are characterised by
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
[
16]? Or, when opting for a time-averaged value
, over what time span should one average? The choice of this time span should be independent of the artificial numerical parameter
.
Channel-wide bottom friction. An essential property of our approach is that the unsteady friction coefficient
is spatially uniform, i.e., it is representative of the entire channel. This is guaranteed for solutions in which the flow velocities
do not vary considerably over the domain. Otherwise, the spatial averaging in Equation (
7) may imply a poor
-estimate. One can test for and resolve this limitation by splitting the channel in several subchannels, each having an individual friction coefficient
. This segmentation is actually a special case of the channel network discussed below. Furthermore, a constant drag coefficient
, 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
on
u and
h.
Interpretation of the unsteady friction coefficient. Given the
-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
versus
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
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,
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, 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 -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
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,
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
, independent of the relaxation parameter
and the initial value
(see
Section 3.3 and
Table 1). The iterative procedure accounts for the dependency of
on
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
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 ; and (ii) a straightforward iterative procedure, carried out in the frequency domain, to obtain both and the corresponding solution of the linearised shallow-water equations. By construction, the unsteady friction coefficient adjusts to the intensity of the flow inside the domain. Qualitatively, 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.