- Email: [email protected]
Non-linear process based modelling of offshore sand waves
Continental Shelf Research 37 (2012) 26–35
Contents lists available at SciVerse ScienceDirect
Continental Shelf Research journal homepage: www.elsevier.com/locate/csr
Non-linear process based modelling of offshore sand waves J. van den Berg a, F. Sterlini b,n, S.J.M.H. Hulscher c, R. van Damme c a b c
MARIN, The Netherlands Water Engineering and Management, Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands University Twente, The Netherlands
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 August 2011 Received in revised form 19 January 2012 Accepted 23 January 2012 Available online 1 February 2012
Sand waves are more or less parallel ridges, making bed patterns in shallow seas. Due to their dimensions and migration rate, they can have a considerable effect on offshore human activities. In relation to such activities especially the variation in the sand wave characteristics (wave length, height and migration) are of importance. Due to the large spatial scale of sand waves and the different time scales involved for the tidal flow (hours) and the bed pattern changes (years to decades), modelling the sand wave behaviour from its initial stage up to its final equilibrium shape is a challenge, let alone to include the spatial variation within a sand wave field. In this paper a numerical solver is presented that describes sand waves from their initial state up to their final equilibrium. It is a 2DV idealized model, based on non-linear stability analysis. Both the model equations and the numerical setup are described and the model is validated against linear stability analysis. To investigate variations in a field of sand waves and the interaction between individual waves, simulations on large domains are presented and show that this solver is able to describe the development of a realistic sand wave field, with variations, from an initially flat sand bed with small random disturbances. The solver shows to be a promising tool, it is computational fast due to the efficient numerical algorithms and is easy to extend with other physical processes. Results show good agreement with analytical results in the linear regime and with field data. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Sand waves Modelling Non-linear
1. Introduction Sand waves form a bed pattern of more or less parallel ridges in shallow seas. For example, in the North Sea, their wave length (i.e. the distance between two adjacent crests) is about 500 m and the height is up to 10 m, which is a considerable amount of the total water depth (ca. 30 m). Sand waves can migrate at velocities if 10 m per year (McCave, 1971; Knaapen et al., 2001), or sometimes even up to speeds of 100 m per year (Buijsman and Ridderinkhof, 2008). The behaviour of sand waves is of both scientific and practical relevance. The latter becomes clear as they can affect the construction and safety of submarine cables and pipelines. On shipping lanes, they can present a danger to navigation. Therefore, these areas need to be monitored and, if necessary, dredged. A detailed overview of the practical relevance of sand waves can be found in Ne´meth et al. (2003).
n
Corresponding author. Tel.: þ31 534892585; fax: þ31 53 4895377. E-mail address: [email protected] (F. Sterlini).
0278-4343/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.csr.2012.01.012
In relation to offshore activities, not only the main characteristics (wave length, height and migration) of sand waves are important but especially their variation, maxima and minima, in both space and time. Sand waves are free instabilities, formed due to interactions between tidal flow and a sandy sea floor (Hulscher, 1996). The tidal flow can change the sea floor shape through sand transport; in turn the shape of the sea floor can then affect the tidal flow, creating a feedback mechanism. A model that can describe the physics of the evolution of sand waves correctly must incorporate descriptions of the tidal flow and the sand transport. Due to the difference in time scales of the tidal flow (hours) and of the bed pattern changes (years to decades) and the large spatial scale of sand waves, modelling the sand wave behaviour from its initial stage up to its final equilibrium shape is a challenge, let alone to include the spatial variation within a sand wave field. Using linear stability analysis on this kind of model, the linear growth of sand waves has been shown by Hulscher (1996), and further investigated by e.g. Gerkema (2000) (investigating the effect of parameters within the system), Besio et al. (2004) (migration and
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
tidal constituents) and Blondeaux and Vittori (2005) (inclusion of extra processes such as wind waves). Linear stability analysis is based on harmonic perturbations in space. What subsequently can be investigated is the effect of different physical conditions (e.g. flow velocity and water depth) on the initial wave length and migration rate. However, linear stability analysis gives no information on the amplitude behaviour or on the shape and migration in the final equilibrium state. Using the same model as Hulscher (1996), the non-linear evolution of a sand wave has been investigated by Ne´meth et al. (2006), using a numerical calculation method. Direct numerical simulation can be used to investigate non-linear amplitude behaviour. So far, it is not fast enough to use on large domains or to work with inclusion of other processes than the ones initially implemented (tidal flow and bed load sediment transport). In combination with the use of periodic boundary conditions, this makes the simulated sand wave length fully dependent on the domain size. In practice, the domain size is chosen according to the fastest growing mode under the circumstances investigated. Because of this, the non-linear behaviour of the wave length cannot be investigated, since the wave length is fixed to one special mode, or multiples of that mode. In addition, interaction between sand waves is impossible, because the simulated setting consists of only one sand wave. Since we want to investigate the variations in a field of sand waves, and the effect of various physical processes such as suspended sediment and surface waves we need to overcome these limitations. We aim to do numerical simulations on large domains. Using large domains, the full non-linear behaviour of a sand wave field, including the interaction between different sand waves, can be investigated. To achieve this goal, we develop an efficient numerical calculation approach. Using this approach, the development of a realistic sand wave field with variations can be achieved, originating from a flat sand bed with small random initial disturbances. We show that the presented solver is fast enough to include a realistic sand wave field and processes other than the basic tidal flow and bed load transport (e.g. suspended sediment transport and surface waves). In the next section we first describe the model that has been used for the numerical simulations. Since the choice of the numerical approach is crucial to this research, a detailed overview of the discretisation and solver is given in Section 3. Next, the validation of the new approach (Section 4) and results of long domain simulations (Section 5) are given. For validation, the numerical accuracy is tested and a comparison is made against linear stability analysis. Since the results from linear stability analysis have already been validated against field data (Hulscher and Van den Brink, 2001), coincidence with linear stability implies coincidence with this field data. Finally, in Sections 6 and 7 the discussion and conclusions are presented.
Fig. 1. One way coupling in the morphodynamic feedback.
Fig. 2. Definitions of H (average height of the water column), hðxÞ (bottom disturbance) and zðx,tÞ (surface elevation).
2.1. The 2D tidal flow For the tidal flow, we start with the equations for incompressible flow: the Navier Stokes equations. These describe the conservation of mass and momentum. In vector notation, they read @r ! ¼ r ðr u Þ @t ! @u ! ! rp ! þ nr2 u þ u ru ¼ r @t
ð1Þ
! where the flow velocity vector is denoted by u , p indicates the pressure scalar, r is the fluid density and n is the kinematic fluid viscosity. This general description can usually be simplified when it is applied to some specific setting. In our setting the following simplifications are made:
We define the z-direction in the line of gravity and x in the
2. Model description To simulate the behaviour of sand waves, it is necessary to model the tidal flow and the sediment transport due to this flow. The difference in time scales, between this tidal flow and the sediment transport, works in our advantage, since it justifies the use of an one-way coupling: during a tidal flow calculation the sea floor is assumed to be static (see Fig. 1). As a starting point, we use identical models as have been used by, for example Hulscher (1996) and Ne´meth et al. (2006). In this section, we first describe the precise class of model we use. In Section 3 we discuss our numerical setup for this problem.
27
principal horizontal direction. The variation in the y-direction is assumed to be negligible in this study. The corresponding flow velocities u and w are in the horizontal and vertical directions respectively (see Fig. 2). The density r is constant. Conservation of mass reduces to the condition that the velocity field should be divergence free: @u @w þ ¼0 @x @z
ð2Þ
Because the water is shallow compared with the characteristic length of a sand wave, the flow is nearly horizontal. Therefore, the vertical accelerations are small, which implies that the pressure can be assumed to be hydrostatic. That is, the pressure difference between any two points on the same
28
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
vertical line depends only on the weight of the fluid between those points (see for example Pedlosky, 1987). Using this hydrostatic assumption, rp in Eq. (1) reduces to gð@z=@xÞ, where g is the constant of gravity and z indicates the water surface elevation. Next, mathematical consistency of the model makes it necessary to drop the momentum equation in the vertical direction. The constant Av is the eddy viscosity in the vertical direction. In the horizontal direction, the turbulent viscosity is assumed negligible. An external body force F ðtÞ is added to mimic the tidal forcing. In Section 2.1.1 the choice of F ðtÞ is explained. Conservation of momentum now reads @u @u @u @2 u @z þu þw ¼ Av 2 g þF ðtÞ @t @x @z @x @z
ð3Þ
The sea floor topography changes over a timescale of years, whereas the tidal flow changes over a timescale of hours. It is therefore legitimate to assume the sea floor elevation hðxÞ is constant within the timeframe of the flow calculation. Next, we discuss the boundary conditions. On the bottom, we apply a partial slip condition @u Av ¼ Su ð4Þ @z z ¼ hðxÞ where Av represents the eddy viscosity and S the slip parameter. The limit S-1 would amount to a no-slip condition. As described by Soulsby (1990), partial slip in combination with a constant eddy viscosity results in a good approximation of the tidalcurrent boundary layer. There is no flow through the bottom, which results in @h ¼ 0 wu ð5Þ @x z ¼ hðxÞ The boundary conditions at the free surface are such that there is no friction at the free surface: @u ¼ 0 ð6Þ @z z ¼ H þ zðx,tÞ and no flow through the free surface: @z @z w¼ þu @t @x
unperturbed sea floor, it generates the tidal flow conditions we want to investigate. We first compute the analytic solution for a flat sea floor. For hðxÞ 0 the solution is independent of x, i.e. all derivatives with respect to x are zero and z 0. From mass conservation (2) we then know that @w ¼0 @z In combination with Eq. (5) this results in w ¼ 0. In a similar fashion, Eq. (3) reduces to an one-dimensional heat equation @uðz,tÞ @2 uðz,tÞ ¼ Av þ F ðtÞ @t @z2
ð8Þ
If we assume an M2 tide, with o corresponding to the tidal frequency Z 1 H u dz ¼ UðtÞ ¼ U 0 þU s sinðotÞ þ U c cosðotÞ H 0 and use the boundary conditions (4) and (6) the linear partial differential equation (8) can be solved analytically, resulting in the value of F ðtÞ: F ðtÞ ¼ F 0 þ F s sinðotÞ þF c cosðotÞ
ð9Þ
The subscripts 0, s, c indicate a unidirectional, sinusoidal and a cosine flow forcing respectively. For each UðtÞ a corresponding F ðtÞ can be found. This makes the flow description complete. 2.2. Sediment transport The evolution of the sea floor is dominated by bed load sediment transport (Blondeaux and Vittori, 2005). Conservation of mass for the sediment gives the evolution of the bottom disturbance h over time: @h @q ¼ @t @x
ð10Þ
where q is the bed load sediment flux. For unidirectional flow over a flat bed, the following empirical relation is commonly accepted: b
q ¼ a9tb 9 ð7Þ
z ¼ H þ zðx,tÞ
For the in and out flow conditions, on the vertical boundaries, we use periodic boundary conditions, as in Ne´meth et al. (2006). In our case this is valid, because we use computation domains that are large enough to limit the influence of this boundary condition. Both in the momentum equation and in the free surface boundary condition, only derivatives of z appear. The free surface is determined up to a constant. This can be imposed mathematically in different ways. The first option is to impose the free surface level at some point in space. However, for periodic boundary conditions it is more convenient to require Z 1 L z dx ¼ 0 L 0 where L is the horizontal length of the periodic domain. 2.1.1. Tidal forcing There are multiple ways to represent the tidal forcing. A common approach is to impose a gradient in the surface elevation. However, since we apply periodic boundary conditions this is not feasible. Instead we impose a forcing F ðtÞ, where t represents time. This forcing is chosen such that, in case of an
The parameters a and b are empirical. Typical values are a ¼ 0:3 and b ¼ 1:5 (see for example Van Rijn, 1993). The specific shear stress at the bottom tb is given by @u tb ¼ Av ð11Þ @z z ¼ hðxÞ Note that the unit of the specific shear stress is m2 =s2 , contrary to the usual unit for shear stress N=m2 . The difference stems from the fact that we divided by the density of the water. The tidal flow that generates sand waves is not unidirectional. Furthermore, when sand waves develop, the bed will no longer be flat. We therefore need to extend this basic formulation. First we acknowledge that sediment transport will, in general, be in the direction of the flow. Secondly, if the slope equals the angle of repose, no uphill transport will be possible. If we define l as the inverse of the tangent of the angle of repose g:
l
1 tan g
we can construct the following linear interpolation: " # tb @h b l q ¼ a9tb 9 @x 9tb 9
ð12Þ
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
The angle of repose g is in the range between 151 and 301, resulting in l 1:524:0. The shear stress should be taken tangential to the bottom. In summary, the flow is described by Eqs. (2) and (3) with boundary conditions (4)–(7) and the sea floor evolution by Eqs. (10) and (11).
29
Similarly, the boundary conditions at the free surface ((6) and (7)) lead to @u^ ¼ 0 ð18Þ ^ @z z^ ¼ H þ z^
^ ¼ w
1 @z^ @z^ þ u^ b @t^ @x^ ^ z ¼ H þ z^
ð19Þ
3. The numerical setup To obtain solutions to the morphodynamic model described in the previous section, a numerical solution procedure has been developed. While the flow description is two-dimensional, the sediment description is only one-dimensional. This implies that the bottleneck in terms of computational effort is in the flow calculation. The setup of the flow solver is as follows: first, a transformation is applied such that the domain with sea floor elevations is mapped to one with a flat sea floor (Section 3.1). Next, the transformed flow equations are discretised in space (Section 3.2) and time (Section 3.3). To solve the resulting discrete nonlinear set of equations, an iterative solver is used (Section 3.3). When the flow solution is obtained, the sea floor topology update is done using a semi-implicit finite difference method (Section 3.4).
3.1. Transformation In general numerical methods work much more efficiently if the mesh is uniform. Therefore, we try to map our continuously changing domain onto a square domain. This changes the bottom topology back to a flat sea floor, such that a rectangular structured grid can be used. Before discussing this transformation, we rewrite Eq. (3) into a conservative form using (2) @u @u2 @wu @2 u @z þ ¼ Av 2 g þ þF ðtÞ @t @z @x @x @z
ð13Þ
Next we use the following transformation: z^ ¼ H
x^ ¼ x,
zhðxÞ , HhðxÞ
t^ ¼ t
Using the following symbols a and b ^ z^ Þ aðx,
^ z^ H @hðxÞ , ^ HhðxÞ @x^
^ ¼ bðxÞ
H ^ HhðxÞ
we can define the following transformed variables: u^ ¼
u
b
,
^ ¼ w
a b
u þ w,
z^ ¼ bz
Together with Eq. (17) this results in Z z^ Z H @u^ @u^ ^ x,H ^ þ z^ ,tÞ ¼ dz dz wð ^ ^ H @x 0 @x Substituting this into Eq. (19) and using a Taylor expansion gives Z H ^ ^ ^ ^ ^ 2 @u^ ^ ¼ 1 @z þ @uðx,H, tÞz þ Oðz^ Þ ^ x,H, ^ tÞ dz^ ¼ wð ð20Þ ^ ^ ^ b @ x @ x @t 0 which is in conservative form. In our situation z^ is small (the 2 Froude number Fr ¼ u2 =gH 5 1), therefore the Oðz^ Þ can be neglected. The complete set of transformed equations is now formed by Eqs. (14)–(18) and (20).
3.2. Discretisation in space For the spatial discretisation, a mimetic method is used (see for example Bochev and Hyman, 2005; Hyman and Shashkov, 1999). The central idea behind mimetic discretisation is to locate the discrete quantities at logical places: fluxes on the cell boundaries and sources in the cell centers. The discrete control volumes are different for each conservation equation. This is unlike a finite volume approach, where all equations are evaluated over the same set of finite volumes. As with a finite volume approach, conservation of mass and momentum is guaranteed at the discrete level. The benefit of using a mimetic grid is that second-order accuracy can be achieved using a local stencil and without spurious oscillations in the solution. ^ þ 12ÞDz^ Þ and the discrete If we position the discrete u^ ij at ðiDx,ðj ^ ij at ðði þ 12ÞDx,j ^ Dz^ Þ (see Fig. 3) we obtain for mass conservation w (Eq. (14)): Z
^ (as proposed By dividing by the Jacobian determinant detðJÞ ¼ bðxÞ by Viviand, 1974; Vinokur, 1974), Eqs. (13) and (2) can be kept in conservative form under the transformation:
ðj þ 1ÞDz^
jDz^
Z
ði þ 1ÞDx^ iDx^
^ @u^ @w þ @x^ @z^
dx^ dz^ ¼ 0
Applying Gauss theorem then leads to ^ i,j w ^ i,j1 ÞDx^ ¼ 0 ðu^ i,j u^ i þ 1,j ÞDz^ þ ðw
^ @u^ @w þ ¼0 @x^ @z^
ð14Þ 2
g @z^
^ F ðtÞ
@u^ @ @ @ u^ 2 ^ uÞ ^ ¼ b2 Av þ ðbu^ Þ þ ðbw þ 2 b @x^ b @x^ @z^ @t^ @z^
ð15Þ
At the bottom, the boundary condition (4) now becomes @u^ ð16Þ Av b ¼ Su^ @z^ z^ ¼ 0 and (5) yields ^ ¼ 09z^ ¼ 0 w
To rewrite this boundary condition into a conservative form, we start by integrating mass conservation over the depth Z H þ z^ ^ @u^ @w þ dz ¼ 0 @x^ @z^ 0
ð17Þ
For the control volume, this is still without any approximation: mass conservation is guaranteed. For the conservation of momentum equation (15), a shifted control volume is used. In contrast to the conservation of mass equation, this equation does have source terms. The control volume is chosen such that these are in the center. We position ^ The source part s of the momentum the discrete z^ i at ðiþ 12ÞDx. equation becomes ! Z ðj þ 1ÞDz^ Z ði þ 1=2ÞDx^ 2^ ^ g @z^ F ðtÞ 2@ u þ Av b s¼ dx^ dz^ 2 b @x^ b jDz^ ði1=2ÞDx^ @z^
30
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
^ and the control volumes for mass Fig. 3. Locations of the discrete u^ and w, conservation (dot dashed) and momentum conservation (dashed).
Using central differences, and assuming the source terms to be constant over the control volume, this then becomes ! ^ ^ i,j þ 1 2u^ i,j þ u^ i,j1 z^ i z^ i1 F ðtÞ 2u s Av b þ g Dx^ Dz^ b Dx^ ðDz^ Þ2
Of course, the quadratic terms in all equations make that all modes are coupled. Combining the discretisation in space and time results in a non-linear set of discrete equations. To find solutions to this system of non-linear equations, a solver based on Newton iteration is used. The sediment transport changes the bottom profile h, resulting in a different tidal flow. Between two bed updates the tidal flows do not change much, since the sediment transport changes the bottom only gradually. Instead of re-calculating the complete flow each time, the Newton iteration is started from the previous solution. For each Newton iteration, the system of equations is solved numerically, either directly using a sparse LU decomposition, or iteratively using a sparse GMRES. As a preconditioner for the GMRES the LU of a previous calculation is used. In this way, GMRES converges within a small number of iterations. The amount of computational work is reduced even further by using incomplete Newton instead of full (Kelley, 1995). 3.4. Numerical treatment of the sand transport equation
For the advection terms a some averaging is needed: u^ i þ 1=2,j 12ðu^ i,j þ u^ i þ 1,j Þ u^ i,j þ 1=2 12ðu^ i,j þ u^ i,j þ 1 Þ ^ i1,j þ w ^ i,j Þ ^ i1=2,j 12ðw w Again using Gauss theorem, the integral formulation Z ðj þ 1ÞDz^ Z ði þ 1=2ÞDx^ @ @ 2 ^ uÞ ^ dx^ dz^ ðbu^ Þ þ ðbw a¼ @x^ @z^ jDz^ ði1=2ÞDx^ is approximated by 2
2
^ i1=2,j u^ i,j þ 1=2 w ^ i1=2,j1 u^ i,j1=2 ÞDx^ a ¼ ðbu^ i þ 1=2,j bu^ i1=2,j ÞDz^ þ bðw The resulting discretisation method is second-order accurate in space. This is validated in Section 4.1. 3.3. Discretisation in time The flow in the North sea is strongly harmonic in time: the physical driving force has a very dominant frequency of around 12.5 h and a weak double frequency. For the discretisation in time, a spectral method is used. We use a harmonic expansion in multiples of the dominant frequency. For 0 1 uðx,z,tÞ B wðx,z,tÞ C fðx,z,tÞ ¼ @ ð21Þ A zðx,tÞ
4. Validation
we write as an approximation:
fðx,z,tÞ f0 ðx,zÞ þ
N X
s
c
ðfk ðx,zÞ sinðkotÞ þ fk ðx,zÞ cosðkotÞÞ
k¼1
where the value of o corresponds with the tidal period. In practice, we truncate this expansion after the first term, i.e. N ¼ 1. In that situation, for every variable in Eq. (21) we have 2N þ 1 ¼ 3 equations and unknowns. As the equations are quadratic in the unknowns, we only have to reduce the product of any two variables. For a non-linear combination fi fj we obtain
fi fj f0i f0j þ 12ðfsi fsj þ fci fcj Þ þ ðf0i fsj þ f0j fsi Þ sinðotÞ 0
c
0
In the simulation of the morphodynamic feedback, the tidal flow is computed separately from the sediment transport. On the time scale of the tidal flow, the sea floor shape can be considered constant. After computing the tidal flow, a sea floor update is done. To update the sea floor shape, the sediment fluxes based on the tidal flow and the local sea floor slope are computed. A mass balance then gives the change of the sea floor shape. As the sea floor change within one tidal cycle is small, the tidal flow is not computed after every tidal cycle. The new sea floor is used with the old tidal cycle to calculate the next sea floor change. The time step for the tidal flow can be either fixed (e.g. 10 weeks) or variable (in which case it depends on the speed of the bed evolution). In case of a fixed time step, attention has to be paid to its length. When this is chosen to be too large, the evolution seems to proceed slower as the slow changes in the beginning seem to carry on, unrealistically, for longer. In the results presented in this paper, time steps of approximately 20 weeks gave reliable results (i.e. approximately the same as for shorter time steps). To compute the sea floor update, we use a semi-implicit finite difference scheme. The sea floor change between the old and the new time level is constructed using the tidal flow solution on the old time level and the sea floor shape on the new time level. For the space discretisation, second order accurate central differences are used. For the time integration we use a first order implicit Euler scheme.
In this section, we describe a number of tests done to validate the model results. The model equations we use are the same as those validated against field data by Hulscher and Van den Brink (2001). Therefore, we do not need to validate the model equations (the model) itself, but just the solution approach (the solver). For ease of comparison we use the same parameter settings as used by Hulscher (1996), listed in Table 1. First the convergence rate and accuracy is checked. Next the linear growth is compared with that found using linear stability analysis. 4.1. Theoretical tests
c
þ ðfi fj þ fj fi Þ cosðotÞ Higher harmonics are neglected. By grouping the sine, cosine and constant terms we obtain three equations for each spatial discrete equation.
The first check is on the order of accuracy in space. In a grid convergence study the results on different grid densities are systematically compared to check the convergence towards the analytical solution of the model equations, for decreasing grid
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
0.4
Table 1 Parameter settings used in a typical simulation run. Dimension
g Av S
9.81 0.03 0.01 1.4 10 4 0.3 1.5 2.5 30.0 0.0, 1.0, 0.0
m/s2 m2/s m/s 1/s s2/m
l H U0 , Us , Uc
m m/s
0.2 (de/dt)/e (1/year)
Value
a b
coarse 1x refined 2x refined analytic
0.3
Parameter
o
31
0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5
sizes. This confirms that the flow solver and the sea floor update are both second order accurate in space. The second test is to validate the accuracy in time. For the sea floor this is again done using a grid convergence study, confirming a first order accuracy in time. For the flow solver the accuracy in time is tested in a different way. The spectral discretisation in time is truncated at the first harmonic. Due to non-linear coupling it could be that the neglected higher harmonics become important. To make sure this is not the case, we also do a test run with a solver that does not use a spectral discretisation in time. Instead, it uses an implicit Euler time stepping. The setup for the test is as follows: the flow over a sine shaped sea floor with length 650 m and amplitude 5 m with an average water depth of 30 m is computed. The used time step is 5 min. After a settling time of several tides, we look at the resulting solution. For each spatial point and each variable we can now compute the spectrum in time over a tidal period. The highest found second harmonic, generated by the non-linear terms, has an amplitude less then 15% of the first harmonic, which we consider acceptable. 4.2. Linear growth In linear stability analysis the growth rate of small harmonic disturbances is investigated, starting from an equilibrium solution. The same analysis can be done numerically. Using the numerical code, the initial growth of a small disturbance can be simulated. From this the growth rate can be computed. By repeating this for a range of wave numbers, we can construct the linear growth curve. Numerical solvers compute approximate solutions to algebraic equations. Depending on the grid density, these solutions always have a truncation error. In Fig. 4 the growth curves for different grid densities are shown. As the grid density increases, the truncation error decreases and the analytic growth curve from linear stability analysis is more faithfully reproduced. Hulscher and Van den Brink (2001) already showed that this growth curve correlates well with sand waves in the North Sea, if the model parameters are tuned to the local conditions.
5. Sand wave development from random initial conditions In this section the results from a simulation on a large domain are presented. These show that a realistic sand wave field, with variations, can develop from a flat bed with small random initial disturbances (of grain size dimensions). In linear stability analysis as done by Hulscher (1996), Gerkema (2000) and Besio et al. (2003), the initial disturbance from which a sand wave develops is harmonic. The same holds for the numerical simulation as done by Ne´meth et al. (2006).
-0.6
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
k=2pi/L (1/m) Fig. 4. Comparison of linear growth found using the numerical approach and linear stability analysis. The numerical results are given for different grid densities. The coarse grid uses 30 cells in the horizontal direction and 12 cells in the vertical. The analytic curve is found using linear stability analysis.
Fig. 5. Given a small domain (600 m) all random initial disturbances evolve into the same sand wave. The results for two different initial random disturbances are plotted. The only real difference in the end results is a phase shift.
In reality the initial disturbance is unknown. Instead of a harmonic perturbation we use a random disturbance (for example the random elevations due to the sand grain shape). If we simulate the evolution of a random disturbance on a small domain (600 m wide), we always find the same sand wave shape as found by Ne´meth et al. (2006) (see Fig. 5, parameter settings in Table 1). The height, shape and length of the simulated sand wave are not dependent of the initial random perturbation. On an average desktop computer, simulating the growth during 40 years takes approximately 2 min (for the domain of 650 m, we used 30 grid points in the x-direction and 22 grid points in the y-direction). Now, we allow sand waves to form on a large domain (i.e. 8000 m), using the same parameter settings (Table 1). This gives a field of sand waves, with a variability in size, shape and growth rate. Again, the time that it takes an average desktop computer to run this simulation for 40 years is approximately 40 min (for this domain length, we used 350 grid points in the x-direction and 22 grid points in the y-direction). In Fig. 6 the sea floor shape after approximately 25 years of evolution is plotted, for different random initial conditions. Typically, waves with a wave length corresponding to the fastest growing mode show the largest amplitudes. In between these optimal waves, there are waves with smaller amplitudes that cannot grow (or shrink) to the optimal wave length because they are enclosed by stronger neighbours. Fig. 7 shows similar results, however with an
32
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
6 4
h (m)
2 0 -2 -4 -6
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000
8 6
h (m)
4 2 0 -2 -4 -6
10 8 6
h (m)
4 2 0 -2 -4 -6 -8
5000
6000
7000
8000
x (m) Fig. 6. From different random initial conditions, a different field of sand waves develops. The three plots each show a field of sand waves from a random disturbance with amplitude 0.30 m after approximately 30 years.
asymmetrical tide. A small unidirectional current is added (0.1 m/s), which results in asymmetric sand waves that develop slightly quicker than the sand waves under symmetric tidal conditions. This results in slightly higher sand waves after the same time period (30 years). Other characteristics are not significantly affected, although earlier research showed a big effect on the sand wave height, when an asymmetric tide was introduced (Sterlini et al., 2009). The same qualitative behaviour is seen in measured sand waves (see Fig. 8). In their raw form, the measurements show a lot of details that were not resolved in the numerical simulations. If we apply a Fourier low pass filter to the data, to remove these small details, the resemblance between the simulation results and the field data becomes even more striking. Fig. 9 shows the results of 10 simulations which only differ in their random initial conditions of the sea bed, together with the field data values. In Fig. 9 the left panel shows the sand wave length results. The shown dots are the minimum, mean and maximum per
simulation. The results are taken at the moment the growth in height significantly slows down (around 30 years of growth, left linear clusters of points in figure) and 10 years after that again (right linear clusters). The figure clearly shows the large variation that the random initial disturbances can create. Sand wave lengths between 320 m and 1300 m can occur, with a mean sand wave length between 500 and 600 m. The grey block represents the minimum and maximum value of the sand wave lengths in the field data. The black line and the open squares represent the mean values of the field data and the simulation results respectively. The results for the sand wave height are shown in the right hand panel of Fig. 9. Again there is a large variation, the height can vary between 2.2 m and 21 m, with a mean around 12 m, which coincides with both the shown model results (Fig. 6) while the field data values are in range with the found minimum heights. These results show that variations within sand wave fields can originate from small random disturbances in the initial conditions.
h (m)
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
12 10 8 6 4 2 0 −2 −4 −6 −8
0
1000
2000
3000
4000
33
5000
6000
7000
8000
x (m) 10 8 6 h (m)
4 2 0 −2 −4 −6
h (m)
−8
12 10 8 6 4 2 0 −2 −4 −6 −8
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
Fig. 7. From different random initial conditions, a different field of sand waves develops. The three plots each show a field of sand waves from a random disturbance with amplitude 0.30 m after approximately 30 years, for asymmetric tidal flow conditions.
8 original filtered
6
h (m)
4 2 0 -2 -4 -6
0
1000
2000
3000
4000
5000
6000
7000
8000
x (m) Fig. 8. Measured sand waves in the North sea. Data taken from the region Kort Verblijf, courtesy RWS Directie Noordzee. The dashed line depicts the same data after Fourier filtering.
34
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
1600
22 20
1400 1200
sand wave height (m)
sand wave length (m)
18
1000 800 600
16 14 12 10 8 6
400 4 200
min
mean max
values per simulation
2
min
mean max
values per simulation
Fig. 9. Statistical results for 10 long domain simulations starting from different randomly disturbed beds. Minimum, mean and maximum sand wave length and height are shown per simulation (open dots) with the mean of these three values (open squares). The results are taken after approximately 30 and 40 years (respectively left and right linear clusters of points for e.g. minimum sand wave length). The grey area represents the field data values (minimum to maximum with the mean represented in the black line).
currents, density gradients, sediment transport) and to describe them within broad classes of problems over different temporal and spatial scales. For example Tonnon et al. (2007) used a complex model to investigate which processes affected an existing sand lump with sand wave dimensions over a time span of 15 years. Due to the complex modelling they were unable to model on a morphological time scale and therefore do not give information on the formation and final state of natural sand waves. In general complex models have two main problems, when aiming for a long timescale and a long domain. Firstly, the computational time is far longer than for an idealized model, such as we use here. Secondly, the numerical diffusion can cause results that are not due to physical processes. Still they can be used for shorter timescales and investigation of specific processes (Borsje et al., in press). The strength of the here presented model is that it is able to model sand waves from initial small bed perturbations up to their full grown shape within a short computational time (order of minutes to tens of minutes). Furthermore the efficient numerical approach and the clear object oriented design processes can be added easily and do not lead to a problematic increase of computational time. It gives the opportunity to investigate the effect of separate processes on sand waves in more detail and to gain insight in leading effects on sand wave characteristics, both in the linear and the non-linear regime. The current aim is not to predict the exact changes in real time sand waves, but to gain insight in the effect of specific processes and correlated minimum and maximum values of e.g. sand wave height and migration.
6. Discussion This solver is already used to study sand waves near the Golden Gate, and showed reasonable agreement with the data when the asymmetry of the tidal current was taken into account (Sterlini et al., 2009). The model is easy to extend with other processes as shown by Van der Meer et al. (2008), Van der Meer et al. (2007), and Roos et al. (2007). Results show that the basic principles are presented correctly in the model and that especially the asymmetry in the tidal current can have a large impact on the sand wave characteristics. Still other processes (e.g. suspended sediment transport and surface waves) can significantly affect the sand waves, and are especially needed to improve the prediction of the final sand wave height. The solver is also used to study river dunes, where similar processes play a role (Hulscher and Dohmen-Janssen, 2005). The model is extended with a parameterization of the flow separation zone and shows good agreement with river dune data (Paarlberg et al., 2009). Due to the efficient numerical algorithm, the model can also be used for operational water management, e.g. for predicting the roughness during high water in rivers (Paarlberg et al., 2010). In shallow seas, the bed is rarely flat, even after dredging. Both larger bed forms (sand banks) and smaller bed forms (mega ripples) can coexist together with sand waves. Although the first can affect the water depth and the second the bed roughness, their influence is not directly included in the current model set up. The used random variation as initial bed is only a rough estimation of what could be found. Still, it can be used as a first indication on the (re)generation time of sand waves. Knaapen and Hulscher (2002) investigated the regeneration of sand waves and assumed an S-shaped growth curve. The results from the here presented model confirm this and could be used to estimate the regeneration of dredged sand waves. As dredging mainly affects the sand wave height and not the length, this process could be studied on a short domain (one wave length). Studying sand waves, one could argue to use a complex or ‘full’ process based models, instead of the idealized model we presented in this paper, as they are developed to include many different processes (e.g. tidal current, wind- and wave-driven
7. Conclusion Offshore sand waves in shallow seas show variability in size, shape and growth rate. Previous theoretical studies of the behaviour of sand waves could not predict this variability. Here, we have presented a solver for sand wave growth and dynamics from their initial small disturbances up to their final equilibrium shape, based on the tidal flow, bed load transport and the feedback mechanism between these processes. The model validation shows good agreement with the analytical results in the linear regime and reasonable accuracy. We use the same model equations as used for the linear stability studies and previous numerical simulations. With the new numerical approach the limitation on the domain size is lifted. The simulation results on large domains show that variations within sand wave fields can originate from small random disturbances in the initial conditions. The solver is able to describe these variations and results show comparable qualitative behaviour with measured field data. Using this solver, other studies show results that agree reasonably well with field and flume data. The model is easy to extent with other physical processes and helps to gain insight in the effect of these processes on sand wave and river dune evolution and final equilibrium. Overall the model is promising as a tool to investigate the effect of physical processes on both sand waves and river dunes. Furthermore, due to the efficient numerical algorithms the model is fast enough for operational use in river dune studies and in large domain sand wave studies. Acknowledgement The research described in this paper is funded by the Technology Foundation STW, applied science division of NWO and the technology program of the Ministry of Economic Affairs (project number TWO.5805: Modelling of spatial and temporal variations in offshore sand waves: process-oriented vs stochastic approach).
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
References Besio, G., Blondeaux, P., Frisina, P., 2003. A note on tidally generated sand waves. Journal of Fluid Dynamics 485, 171–190. Besio, G., Blondeaux, P., Brocchini, M., Vittori, G., 2004. On the modeling of sand wave migration. Journal of Geophysical Research 109, C04018. doi:10.1029/ 2002JC001622. Blondeaux, P., Vittori, G., 2005. Flow and sediment transport induces by tide propagation. 2. The wavy bottom case. Journal of Geophysical Research— Oceans 110, C08003. doi:10.1029/2004JC002545. Bochev, P.B., Hyman, J.M., 2005. Principles of Mimetic Discretizations of Differential Operators. Preprint: LA-UR-05-4244. Borsje, B.W., Roos, P.C., Kranenburg, W.M., Hulscher, S.J.M.H., 2011. Modeling sandwave formation in a numerical shallow water model. In: Shao, X., Wang, Z., Wang, G. (Eds.), RCEM 2011, Seventh IAHR Symposium on RiverCoastal and Estuarine Morphodynamics, Beijing, China. Buijsman, M.C., Ridderinkhof, H., 2008. Long-term evolution of sand waves in the Marsdiep inlet. I. High-resolution observations. Continental Shelf Research 28 (9), 1190–1201. doi:10.1016/j.csr.2007.10.011. Gerkema, T., 2000. A linear stability analysis of tidally generated sand waves. Journal of Fluid Mechanics 417, 303–322. Hulscher, S.J.M.H., 1996. Tidal-induced large-scale regular bed form patterns in a three-dimensional shallow water model. Journal of Geophysical Research 101, 727–744. Hulscher, S.J.M.H., Dohmen-Janssen, C.M., 2005. Introduction to special section on marine sand wave and river dune dynamics. Journal of Geophysical Research 110, F04S01. doi:10.1029/2005JF000404. Hulscher, S.J.M.H., Van den Brink, G.M., 2001. Comparison between predicted and observed sand waves and sand banks in the North Sea. Journal of Geophysical Research 106 (C5), 9327–9338. Hyman, J., Shashkov, M., 1999. The orthogonal decomposition theorems for mimetic finite difference methods. Society for Industrial and Applied Mathematics 36 (3), 788–818. Kelley, C.T., 1995. Iterative Methods for Linear and Nonlinear Equations. SIAM. Knaapen, M.A.F., Hulscher, S.J.M.H., 2002. Regeneration of sand waves after dredging. Coastal Engineering 46, 277–289. Knaapen, M.A.F., Hulscher, S.J.M.H., De Vriend, H.J., 2001. A new type of sea bed waves. Geophysical Research Letters 28 (7), 1323–1326. McCave, I.N., 1971. Sand waves in the North Sea off the coast of Holland. Marine Geology 10 (3), 199–225.
35
Ne´meth, A.A., Hulscher, S.J.M.H., De Vriend, H.J., 2003. Offshore sand wave dynamics, engineering problems and future solutions. Pipeline Gas Journal 230 (4), 67–69. Ne´meth, A.A., Hulscher, S.J.M.H., Van Damme, R.M.J., 2006. Simulating offshore sand waves. Coastal Engineering 53, 265–275. Paarlberg, A.J., Dohmen-Janssen, C.M., Hulscher, S.J.M.H., Termes, P., 2009. Modeling river dune evolution using a parameterization of flow separation. Journal of Geophysical Research—Earth Surface 114, F01014. Paarlberg, A.J., Dohmen-Janssen, C.M., Hulscher, S.J.M.H., Termes, P., Schielen, R., 2010. Modelling the effect of time-dependent river dune evolution on bed roughness and stage. Earth Surface Processes and Landforms 35, 1854–1866. doi:10.1002/esp.2074. Pedlosky, J., 1987. Geophysical Fluid Dynamics. Springer, New York. Roos, P.C., Hulscher, S.J.M.H., Meer, F.M.v.d., Dijk, T.A.G.P.v., Wientjes, I.G.M., Berg, J.v.d., 2007. Grain size sorting over offshore sand waves: observations and modelling. In: Dohmen-Janssen, C.M., Hulscher, S.J.M.H. (Eds.), River, Coastal and Estuarine Morphodynamics: RCEM 2007. Taylor & Francis Group, London. Soulsby, R.L., 1990. Tidal-current boundary layers. In: le Mehaute, B., Hanes, D.M. (Eds.), The Sea, Ocean Engineering Science. John Wiley and Sons, New York, pp. 523–566. Sterlini, F., Hulscher, S.J.H.M., Hanes, D.M., 2009. Simulating and understanding sand wave variation: a case study of the Golden Gate sand waves. Journal of Geophysical Research 114, F02007. doi:10.1029/2008JF000999. Tonnon, P.K., Van-Rijn, L.C., Walstra, D.J.R., 2007. The morphodynamic modelling of tidal sand waves on the shoreface. Coastal Engineering 54, 279–296. Van der Meer, F., Hulscher, S.J.M.H., van den Berg, J., 2007. On the influence of suspended sediment transport on the generation of offshore sand waves. In: B.J.G., et al. (Eds.), Particle Laden Flow: From Geophysical to Kolmogorov Scales. Springer, pp. 29–41. Van der Meer, F., Hulscher, S.J.M.H., Dodd, N., 2008. Modelling surface wave effects on evolved offshore sand waves. In: Souza, Alejandro J., Proctor, R. (Eds.), PECS 2008: Physics of Estuaries and Coastal Seas, Liverpool, UK. Van Rijn, L.C., 1993. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas, vol. III. Aqua Publications, Amsterdam. Vinokur, M., 1974. Conservation equations of gas-dynamics in curvilinear coordinate systems. Journal of Computational Physics 14, 105–125. Viviand, H., 1974. Conservative Forms of Gas Dynamic Equations, vol. 1. Rech. Aerosp., pp. 65–68.
Contents lists available at SciVerse ScienceDirect
Continental Shelf Research journal homepage: www.elsevier.com/locate/csr
Non-linear process based modelling of offshore sand waves J. van den Berg a, F. Sterlini b,n, S.J.M.H. Hulscher c, R. van Damme c a b c
MARIN, The Netherlands Water Engineering and Management, Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands University Twente, The Netherlands
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 August 2011 Received in revised form 19 January 2012 Accepted 23 January 2012 Available online 1 February 2012
Sand waves are more or less parallel ridges, making bed patterns in shallow seas. Due to their dimensions and migration rate, they can have a considerable effect on offshore human activities. In relation to such activities especially the variation in the sand wave characteristics (wave length, height and migration) are of importance. Due to the large spatial scale of sand waves and the different time scales involved for the tidal flow (hours) and the bed pattern changes (years to decades), modelling the sand wave behaviour from its initial stage up to its final equilibrium shape is a challenge, let alone to include the spatial variation within a sand wave field. In this paper a numerical solver is presented that describes sand waves from their initial state up to their final equilibrium. It is a 2DV idealized model, based on non-linear stability analysis. Both the model equations and the numerical setup are described and the model is validated against linear stability analysis. To investigate variations in a field of sand waves and the interaction between individual waves, simulations on large domains are presented and show that this solver is able to describe the development of a realistic sand wave field, with variations, from an initially flat sand bed with small random disturbances. The solver shows to be a promising tool, it is computational fast due to the efficient numerical algorithms and is easy to extend with other physical processes. Results show good agreement with analytical results in the linear regime and with field data. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Sand waves Modelling Non-linear
1. Introduction Sand waves form a bed pattern of more or less parallel ridges in shallow seas. For example, in the North Sea, their wave length (i.e. the distance between two adjacent crests) is about 500 m and the height is up to 10 m, which is a considerable amount of the total water depth (ca. 30 m). Sand waves can migrate at velocities if 10 m per year (McCave, 1971; Knaapen et al., 2001), or sometimes even up to speeds of 100 m per year (Buijsman and Ridderinkhof, 2008). The behaviour of sand waves is of both scientific and practical relevance. The latter becomes clear as they can affect the construction and safety of submarine cables and pipelines. On shipping lanes, they can present a danger to navigation. Therefore, these areas need to be monitored and, if necessary, dredged. A detailed overview of the practical relevance of sand waves can be found in Ne´meth et al. (2003).
n
Corresponding author. Tel.: þ31 534892585; fax: þ31 53 4895377. E-mail address: [email protected] (F. Sterlini).
0278-4343/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.csr.2012.01.012
In relation to offshore activities, not only the main characteristics (wave length, height and migration) of sand waves are important but especially their variation, maxima and minima, in both space and time. Sand waves are free instabilities, formed due to interactions between tidal flow and a sandy sea floor (Hulscher, 1996). The tidal flow can change the sea floor shape through sand transport; in turn the shape of the sea floor can then affect the tidal flow, creating a feedback mechanism. A model that can describe the physics of the evolution of sand waves correctly must incorporate descriptions of the tidal flow and the sand transport. Due to the difference in time scales of the tidal flow (hours) and of the bed pattern changes (years to decades) and the large spatial scale of sand waves, modelling the sand wave behaviour from its initial stage up to its final equilibrium shape is a challenge, let alone to include the spatial variation within a sand wave field. Using linear stability analysis on this kind of model, the linear growth of sand waves has been shown by Hulscher (1996), and further investigated by e.g. Gerkema (2000) (investigating the effect of parameters within the system), Besio et al. (2004) (migration and
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
tidal constituents) and Blondeaux and Vittori (2005) (inclusion of extra processes such as wind waves). Linear stability analysis is based on harmonic perturbations in space. What subsequently can be investigated is the effect of different physical conditions (e.g. flow velocity and water depth) on the initial wave length and migration rate. However, linear stability analysis gives no information on the amplitude behaviour or on the shape and migration in the final equilibrium state. Using the same model as Hulscher (1996), the non-linear evolution of a sand wave has been investigated by Ne´meth et al. (2006), using a numerical calculation method. Direct numerical simulation can be used to investigate non-linear amplitude behaviour. So far, it is not fast enough to use on large domains or to work with inclusion of other processes than the ones initially implemented (tidal flow and bed load sediment transport). In combination with the use of periodic boundary conditions, this makes the simulated sand wave length fully dependent on the domain size. In practice, the domain size is chosen according to the fastest growing mode under the circumstances investigated. Because of this, the non-linear behaviour of the wave length cannot be investigated, since the wave length is fixed to one special mode, or multiples of that mode. In addition, interaction between sand waves is impossible, because the simulated setting consists of only one sand wave. Since we want to investigate the variations in a field of sand waves, and the effect of various physical processes such as suspended sediment and surface waves we need to overcome these limitations. We aim to do numerical simulations on large domains. Using large domains, the full non-linear behaviour of a sand wave field, including the interaction between different sand waves, can be investigated. To achieve this goal, we develop an efficient numerical calculation approach. Using this approach, the development of a realistic sand wave field with variations can be achieved, originating from a flat sand bed with small random initial disturbances. We show that the presented solver is fast enough to include a realistic sand wave field and processes other than the basic tidal flow and bed load transport (e.g. suspended sediment transport and surface waves). In the next section we first describe the model that has been used for the numerical simulations. Since the choice of the numerical approach is crucial to this research, a detailed overview of the discretisation and solver is given in Section 3. Next, the validation of the new approach (Section 4) and results of long domain simulations (Section 5) are given. For validation, the numerical accuracy is tested and a comparison is made against linear stability analysis. Since the results from linear stability analysis have already been validated against field data (Hulscher and Van den Brink, 2001), coincidence with linear stability implies coincidence with this field data. Finally, in Sections 6 and 7 the discussion and conclusions are presented.
Fig. 1. One way coupling in the morphodynamic feedback.
Fig. 2. Definitions of H (average height of the water column), hðxÞ (bottom disturbance) and zðx,tÞ (surface elevation).
2.1. The 2D tidal flow For the tidal flow, we start with the equations for incompressible flow: the Navier Stokes equations. These describe the conservation of mass and momentum. In vector notation, they read @r ! ¼ r ðr u Þ @t ! @u ! ! rp ! þ nr2 u þ u ru ¼ r @t
ð1Þ
! where the flow velocity vector is denoted by u , p indicates the pressure scalar, r is the fluid density and n is the kinematic fluid viscosity. This general description can usually be simplified when it is applied to some specific setting. In our setting the following simplifications are made:
We define the z-direction in the line of gravity and x in the
2. Model description To simulate the behaviour of sand waves, it is necessary to model the tidal flow and the sediment transport due to this flow. The difference in time scales, between this tidal flow and the sediment transport, works in our advantage, since it justifies the use of an one-way coupling: during a tidal flow calculation the sea floor is assumed to be static (see Fig. 1). As a starting point, we use identical models as have been used by, for example Hulscher (1996) and Ne´meth et al. (2006). In this section, we first describe the precise class of model we use. In Section 3 we discuss our numerical setup for this problem.
27
principal horizontal direction. The variation in the y-direction is assumed to be negligible in this study. The corresponding flow velocities u and w are in the horizontal and vertical directions respectively (see Fig. 2). The density r is constant. Conservation of mass reduces to the condition that the velocity field should be divergence free: @u @w þ ¼0 @x @z
ð2Þ
Because the water is shallow compared with the characteristic length of a sand wave, the flow is nearly horizontal. Therefore, the vertical accelerations are small, which implies that the pressure can be assumed to be hydrostatic. That is, the pressure difference between any two points on the same
28
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
vertical line depends only on the weight of the fluid between those points (see for example Pedlosky, 1987). Using this hydrostatic assumption, rp in Eq. (1) reduces to gð@z=@xÞ, where g is the constant of gravity and z indicates the water surface elevation. Next, mathematical consistency of the model makes it necessary to drop the momentum equation in the vertical direction. The constant Av is the eddy viscosity in the vertical direction. In the horizontal direction, the turbulent viscosity is assumed negligible. An external body force F ðtÞ is added to mimic the tidal forcing. In Section 2.1.1 the choice of F ðtÞ is explained. Conservation of momentum now reads @u @u @u @2 u @z þu þw ¼ Av 2 g þF ðtÞ @t @x @z @x @z
ð3Þ
The sea floor topography changes over a timescale of years, whereas the tidal flow changes over a timescale of hours. It is therefore legitimate to assume the sea floor elevation hðxÞ is constant within the timeframe of the flow calculation. Next, we discuss the boundary conditions. On the bottom, we apply a partial slip condition @u Av ¼ Su ð4Þ @z z ¼ hðxÞ where Av represents the eddy viscosity and S the slip parameter. The limit S-1 would amount to a no-slip condition. As described by Soulsby (1990), partial slip in combination with a constant eddy viscosity results in a good approximation of the tidalcurrent boundary layer. There is no flow through the bottom, which results in @h ¼ 0 wu ð5Þ @x z ¼ hðxÞ The boundary conditions at the free surface are such that there is no friction at the free surface: @u ¼ 0 ð6Þ @z z ¼ H þ zðx,tÞ and no flow through the free surface: @z @z w¼ þu @t @x
unperturbed sea floor, it generates the tidal flow conditions we want to investigate. We first compute the analytic solution for a flat sea floor. For hðxÞ 0 the solution is independent of x, i.e. all derivatives with respect to x are zero and z 0. From mass conservation (2) we then know that @w ¼0 @z In combination with Eq. (5) this results in w ¼ 0. In a similar fashion, Eq. (3) reduces to an one-dimensional heat equation @uðz,tÞ @2 uðz,tÞ ¼ Av þ F ðtÞ @t @z2
ð8Þ
If we assume an M2 tide, with o corresponding to the tidal frequency Z 1 H u dz ¼ UðtÞ ¼ U 0 þU s sinðotÞ þ U c cosðotÞ H 0 and use the boundary conditions (4) and (6) the linear partial differential equation (8) can be solved analytically, resulting in the value of F ðtÞ: F ðtÞ ¼ F 0 þ F s sinðotÞ þF c cosðotÞ
ð9Þ
The subscripts 0, s, c indicate a unidirectional, sinusoidal and a cosine flow forcing respectively. For each UðtÞ a corresponding F ðtÞ can be found. This makes the flow description complete. 2.2. Sediment transport The evolution of the sea floor is dominated by bed load sediment transport (Blondeaux and Vittori, 2005). Conservation of mass for the sediment gives the evolution of the bottom disturbance h over time: @h @q ¼ @t @x
ð10Þ
where q is the bed load sediment flux. For unidirectional flow over a flat bed, the following empirical relation is commonly accepted: b
q ¼ a9tb 9 ð7Þ
z ¼ H þ zðx,tÞ
For the in and out flow conditions, on the vertical boundaries, we use periodic boundary conditions, as in Ne´meth et al. (2006). In our case this is valid, because we use computation domains that are large enough to limit the influence of this boundary condition. Both in the momentum equation and in the free surface boundary condition, only derivatives of z appear. The free surface is determined up to a constant. This can be imposed mathematically in different ways. The first option is to impose the free surface level at some point in space. However, for periodic boundary conditions it is more convenient to require Z 1 L z dx ¼ 0 L 0 where L is the horizontal length of the periodic domain. 2.1.1. Tidal forcing There are multiple ways to represent the tidal forcing. A common approach is to impose a gradient in the surface elevation. However, since we apply periodic boundary conditions this is not feasible. Instead we impose a forcing F ðtÞ, where t represents time. This forcing is chosen such that, in case of an
The parameters a and b are empirical. Typical values are a ¼ 0:3 and b ¼ 1:5 (see for example Van Rijn, 1993). The specific shear stress at the bottom tb is given by @u tb ¼ Av ð11Þ @z z ¼ hðxÞ Note that the unit of the specific shear stress is m2 =s2 , contrary to the usual unit for shear stress N=m2 . The difference stems from the fact that we divided by the density of the water. The tidal flow that generates sand waves is not unidirectional. Furthermore, when sand waves develop, the bed will no longer be flat. We therefore need to extend this basic formulation. First we acknowledge that sediment transport will, in general, be in the direction of the flow. Secondly, if the slope equals the angle of repose, no uphill transport will be possible. If we define l as the inverse of the tangent of the angle of repose g:
l
1 tan g
we can construct the following linear interpolation: " # tb @h b l q ¼ a9tb 9 @x 9tb 9
ð12Þ
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
The angle of repose g is in the range between 151 and 301, resulting in l 1:524:0. The shear stress should be taken tangential to the bottom. In summary, the flow is described by Eqs. (2) and (3) with boundary conditions (4)–(7) and the sea floor evolution by Eqs. (10) and (11).
29
Similarly, the boundary conditions at the free surface ((6) and (7)) lead to @u^ ¼ 0 ð18Þ ^ @z z^ ¼ H þ z^
^ ¼ w
1 @z^ @z^ þ u^ b @t^ @x^ ^ z ¼ H þ z^
ð19Þ
3. The numerical setup To obtain solutions to the morphodynamic model described in the previous section, a numerical solution procedure has been developed. While the flow description is two-dimensional, the sediment description is only one-dimensional. This implies that the bottleneck in terms of computational effort is in the flow calculation. The setup of the flow solver is as follows: first, a transformation is applied such that the domain with sea floor elevations is mapped to one with a flat sea floor (Section 3.1). Next, the transformed flow equations are discretised in space (Section 3.2) and time (Section 3.3). To solve the resulting discrete nonlinear set of equations, an iterative solver is used (Section 3.3). When the flow solution is obtained, the sea floor topology update is done using a semi-implicit finite difference method (Section 3.4).
3.1. Transformation In general numerical methods work much more efficiently if the mesh is uniform. Therefore, we try to map our continuously changing domain onto a square domain. This changes the bottom topology back to a flat sea floor, such that a rectangular structured grid can be used. Before discussing this transformation, we rewrite Eq. (3) into a conservative form using (2) @u @u2 @wu @2 u @z þ ¼ Av 2 g þ þF ðtÞ @t @z @x @x @z
ð13Þ
Next we use the following transformation: z^ ¼ H
x^ ¼ x,
zhðxÞ , HhðxÞ
t^ ¼ t
Using the following symbols a and b ^ z^ Þ aðx,
^ z^ H @hðxÞ , ^ HhðxÞ @x^
^ ¼ bðxÞ
H ^ HhðxÞ
we can define the following transformed variables: u^ ¼
u
b
,
^ ¼ w
a b
u þ w,
z^ ¼ bz
Together with Eq. (17) this results in Z z^ Z H @u^ @u^ ^ x,H ^ þ z^ ,tÞ ¼ dz dz wð ^ ^ H @x 0 @x Substituting this into Eq. (19) and using a Taylor expansion gives Z H ^ ^ ^ ^ ^ 2 @u^ ^ ¼ 1 @z þ @uðx,H, tÞz þ Oðz^ Þ ^ x,H, ^ tÞ dz^ ¼ wð ð20Þ ^ ^ ^ b @ x @ x @t 0 which is in conservative form. In our situation z^ is small (the 2 Froude number Fr ¼ u2 =gH 5 1), therefore the Oðz^ Þ can be neglected. The complete set of transformed equations is now formed by Eqs. (14)–(18) and (20).
3.2. Discretisation in space For the spatial discretisation, a mimetic method is used (see for example Bochev and Hyman, 2005; Hyman and Shashkov, 1999). The central idea behind mimetic discretisation is to locate the discrete quantities at logical places: fluxes on the cell boundaries and sources in the cell centers. The discrete control volumes are different for each conservation equation. This is unlike a finite volume approach, where all equations are evaluated over the same set of finite volumes. As with a finite volume approach, conservation of mass and momentum is guaranteed at the discrete level. The benefit of using a mimetic grid is that second-order accuracy can be achieved using a local stencil and without spurious oscillations in the solution. ^ þ 12ÞDz^ Þ and the discrete If we position the discrete u^ ij at ðiDx,ðj ^ ij at ðði þ 12ÞDx,j ^ Dz^ Þ (see Fig. 3) we obtain for mass conservation w (Eq. (14)): Z
^ (as proposed By dividing by the Jacobian determinant detðJÞ ¼ bðxÞ by Viviand, 1974; Vinokur, 1974), Eqs. (13) and (2) can be kept in conservative form under the transformation:
ðj þ 1ÞDz^
jDz^
Z
ði þ 1ÞDx^ iDx^
^ @u^ @w þ @x^ @z^
dx^ dz^ ¼ 0
Applying Gauss theorem then leads to ^ i,j w ^ i,j1 ÞDx^ ¼ 0 ðu^ i,j u^ i þ 1,j ÞDz^ þ ðw
^ @u^ @w þ ¼0 @x^ @z^
ð14Þ 2
g @z^
^ F ðtÞ
@u^ @ @ @ u^ 2 ^ uÞ ^ ¼ b2 Av þ ðbu^ Þ þ ðbw þ 2 b @x^ b @x^ @z^ @t^ @z^
ð15Þ
At the bottom, the boundary condition (4) now becomes @u^ ð16Þ Av b ¼ Su^ @z^ z^ ¼ 0 and (5) yields ^ ¼ 09z^ ¼ 0 w
To rewrite this boundary condition into a conservative form, we start by integrating mass conservation over the depth Z H þ z^ ^ @u^ @w þ dz ¼ 0 @x^ @z^ 0
ð17Þ
For the control volume, this is still without any approximation: mass conservation is guaranteed. For the conservation of momentum equation (15), a shifted control volume is used. In contrast to the conservation of mass equation, this equation does have source terms. The control volume is chosen such that these are in the center. We position ^ The source part s of the momentum the discrete z^ i at ðiþ 12ÞDx. equation becomes ! Z ðj þ 1ÞDz^ Z ði þ 1=2ÞDx^ 2^ ^ g @z^ F ðtÞ 2@ u þ Av b s¼ dx^ dz^ 2 b @x^ b jDz^ ði1=2ÞDx^ @z^
30
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
^ and the control volumes for mass Fig. 3. Locations of the discrete u^ and w, conservation (dot dashed) and momentum conservation (dashed).
Using central differences, and assuming the source terms to be constant over the control volume, this then becomes ! ^ ^ i,j þ 1 2u^ i,j þ u^ i,j1 z^ i z^ i1 F ðtÞ 2u s Av b þ g Dx^ Dz^ b Dx^ ðDz^ Þ2
Of course, the quadratic terms in all equations make that all modes are coupled. Combining the discretisation in space and time results in a non-linear set of discrete equations. To find solutions to this system of non-linear equations, a solver based on Newton iteration is used. The sediment transport changes the bottom profile h, resulting in a different tidal flow. Between two bed updates the tidal flows do not change much, since the sediment transport changes the bottom only gradually. Instead of re-calculating the complete flow each time, the Newton iteration is started from the previous solution. For each Newton iteration, the system of equations is solved numerically, either directly using a sparse LU decomposition, or iteratively using a sparse GMRES. As a preconditioner for the GMRES the LU of a previous calculation is used. In this way, GMRES converges within a small number of iterations. The amount of computational work is reduced even further by using incomplete Newton instead of full (Kelley, 1995). 3.4. Numerical treatment of the sand transport equation
For the advection terms a some averaging is needed: u^ i þ 1=2,j 12ðu^ i,j þ u^ i þ 1,j Þ u^ i,j þ 1=2 12ðu^ i,j þ u^ i,j þ 1 Þ ^ i1,j þ w ^ i,j Þ ^ i1=2,j 12ðw w Again using Gauss theorem, the integral formulation Z ðj þ 1ÞDz^ Z ði þ 1=2ÞDx^ @ @ 2 ^ uÞ ^ dx^ dz^ ðbu^ Þ þ ðbw a¼ @x^ @z^ jDz^ ði1=2ÞDx^ is approximated by 2
2
^ i1=2,j u^ i,j þ 1=2 w ^ i1=2,j1 u^ i,j1=2 ÞDx^ a ¼ ðbu^ i þ 1=2,j bu^ i1=2,j ÞDz^ þ bðw The resulting discretisation method is second-order accurate in space. This is validated in Section 4.1. 3.3. Discretisation in time The flow in the North sea is strongly harmonic in time: the physical driving force has a very dominant frequency of around 12.5 h and a weak double frequency. For the discretisation in time, a spectral method is used. We use a harmonic expansion in multiples of the dominant frequency. For 0 1 uðx,z,tÞ B wðx,z,tÞ C fðx,z,tÞ ¼ @ ð21Þ A zðx,tÞ
4. Validation
we write as an approximation:
fðx,z,tÞ f0 ðx,zÞ þ
N X
s
c
ðfk ðx,zÞ sinðkotÞ þ fk ðx,zÞ cosðkotÞÞ
k¼1
where the value of o corresponds with the tidal period. In practice, we truncate this expansion after the first term, i.e. N ¼ 1. In that situation, for every variable in Eq. (21) we have 2N þ 1 ¼ 3 equations and unknowns. As the equations are quadratic in the unknowns, we only have to reduce the product of any two variables. For a non-linear combination fi fj we obtain
fi fj f0i f0j þ 12ðfsi fsj þ fci fcj Þ þ ðf0i fsj þ f0j fsi Þ sinðotÞ 0
c
0
In the simulation of the morphodynamic feedback, the tidal flow is computed separately from the sediment transport. On the time scale of the tidal flow, the sea floor shape can be considered constant. After computing the tidal flow, a sea floor update is done. To update the sea floor shape, the sediment fluxes based on the tidal flow and the local sea floor slope are computed. A mass balance then gives the change of the sea floor shape. As the sea floor change within one tidal cycle is small, the tidal flow is not computed after every tidal cycle. The new sea floor is used with the old tidal cycle to calculate the next sea floor change. The time step for the tidal flow can be either fixed (e.g. 10 weeks) or variable (in which case it depends on the speed of the bed evolution). In case of a fixed time step, attention has to be paid to its length. When this is chosen to be too large, the evolution seems to proceed slower as the slow changes in the beginning seem to carry on, unrealistically, for longer. In the results presented in this paper, time steps of approximately 20 weeks gave reliable results (i.e. approximately the same as for shorter time steps). To compute the sea floor update, we use a semi-implicit finite difference scheme. The sea floor change between the old and the new time level is constructed using the tidal flow solution on the old time level and the sea floor shape on the new time level. For the space discretisation, second order accurate central differences are used. For the time integration we use a first order implicit Euler scheme.
In this section, we describe a number of tests done to validate the model results. The model equations we use are the same as those validated against field data by Hulscher and Van den Brink (2001). Therefore, we do not need to validate the model equations (the model) itself, but just the solution approach (the solver). For ease of comparison we use the same parameter settings as used by Hulscher (1996), listed in Table 1. First the convergence rate and accuracy is checked. Next the linear growth is compared with that found using linear stability analysis. 4.1. Theoretical tests
c
þ ðfi fj þ fj fi Þ cosðotÞ Higher harmonics are neglected. By grouping the sine, cosine and constant terms we obtain three equations for each spatial discrete equation.
The first check is on the order of accuracy in space. In a grid convergence study the results on different grid densities are systematically compared to check the convergence towards the analytical solution of the model equations, for decreasing grid
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
0.4
Table 1 Parameter settings used in a typical simulation run. Dimension
g Av S
9.81 0.03 0.01 1.4 10 4 0.3 1.5 2.5 30.0 0.0, 1.0, 0.0
m/s2 m2/s m/s 1/s s2/m
l H U0 , Us , Uc
m m/s
0.2 (de/dt)/e (1/year)
Value
a b
coarse 1x refined 2x refined analytic
0.3
Parameter
o
31
0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5
sizes. This confirms that the flow solver and the sea floor update are both second order accurate in space. The second test is to validate the accuracy in time. For the sea floor this is again done using a grid convergence study, confirming a first order accuracy in time. For the flow solver the accuracy in time is tested in a different way. The spectral discretisation in time is truncated at the first harmonic. Due to non-linear coupling it could be that the neglected higher harmonics become important. To make sure this is not the case, we also do a test run with a solver that does not use a spectral discretisation in time. Instead, it uses an implicit Euler time stepping. The setup for the test is as follows: the flow over a sine shaped sea floor with length 650 m and amplitude 5 m with an average water depth of 30 m is computed. The used time step is 5 min. After a settling time of several tides, we look at the resulting solution. For each spatial point and each variable we can now compute the spectrum in time over a tidal period. The highest found second harmonic, generated by the non-linear terms, has an amplitude less then 15% of the first harmonic, which we consider acceptable. 4.2. Linear growth In linear stability analysis the growth rate of small harmonic disturbances is investigated, starting from an equilibrium solution. The same analysis can be done numerically. Using the numerical code, the initial growth of a small disturbance can be simulated. From this the growth rate can be computed. By repeating this for a range of wave numbers, we can construct the linear growth curve. Numerical solvers compute approximate solutions to algebraic equations. Depending on the grid density, these solutions always have a truncation error. In Fig. 4 the growth curves for different grid densities are shown. As the grid density increases, the truncation error decreases and the analytic growth curve from linear stability analysis is more faithfully reproduced. Hulscher and Van den Brink (2001) already showed that this growth curve correlates well with sand waves in the North Sea, if the model parameters are tuned to the local conditions.
5. Sand wave development from random initial conditions In this section the results from a simulation on a large domain are presented. These show that a realistic sand wave field, with variations, can develop from a flat bed with small random initial disturbances (of grain size dimensions). In linear stability analysis as done by Hulscher (1996), Gerkema (2000) and Besio et al. (2003), the initial disturbance from which a sand wave develops is harmonic. The same holds for the numerical simulation as done by Ne´meth et al. (2006).
-0.6
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
k=2pi/L (1/m) Fig. 4. Comparison of linear growth found using the numerical approach and linear stability analysis. The numerical results are given for different grid densities. The coarse grid uses 30 cells in the horizontal direction and 12 cells in the vertical. The analytic curve is found using linear stability analysis.
Fig. 5. Given a small domain (600 m) all random initial disturbances evolve into the same sand wave. The results for two different initial random disturbances are plotted. The only real difference in the end results is a phase shift.
In reality the initial disturbance is unknown. Instead of a harmonic perturbation we use a random disturbance (for example the random elevations due to the sand grain shape). If we simulate the evolution of a random disturbance on a small domain (600 m wide), we always find the same sand wave shape as found by Ne´meth et al. (2006) (see Fig. 5, parameter settings in Table 1). The height, shape and length of the simulated sand wave are not dependent of the initial random perturbation. On an average desktop computer, simulating the growth during 40 years takes approximately 2 min (for the domain of 650 m, we used 30 grid points in the x-direction and 22 grid points in the y-direction). Now, we allow sand waves to form on a large domain (i.e. 8000 m), using the same parameter settings (Table 1). This gives a field of sand waves, with a variability in size, shape and growth rate. Again, the time that it takes an average desktop computer to run this simulation for 40 years is approximately 40 min (for this domain length, we used 350 grid points in the x-direction and 22 grid points in the y-direction). In Fig. 6 the sea floor shape after approximately 25 years of evolution is plotted, for different random initial conditions. Typically, waves with a wave length corresponding to the fastest growing mode show the largest amplitudes. In between these optimal waves, there are waves with smaller amplitudes that cannot grow (or shrink) to the optimal wave length because they are enclosed by stronger neighbours. Fig. 7 shows similar results, however with an
32
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
6 4
h (m)
2 0 -2 -4 -6
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000
8 6
h (m)
4 2 0 -2 -4 -6
10 8 6
h (m)
4 2 0 -2 -4 -6 -8
5000
6000
7000
8000
x (m) Fig. 6. From different random initial conditions, a different field of sand waves develops. The three plots each show a field of sand waves from a random disturbance with amplitude 0.30 m after approximately 30 years.
asymmetrical tide. A small unidirectional current is added (0.1 m/s), which results in asymmetric sand waves that develop slightly quicker than the sand waves under symmetric tidal conditions. This results in slightly higher sand waves after the same time period (30 years). Other characteristics are not significantly affected, although earlier research showed a big effect on the sand wave height, when an asymmetric tide was introduced (Sterlini et al., 2009). The same qualitative behaviour is seen in measured sand waves (see Fig. 8). In their raw form, the measurements show a lot of details that were not resolved in the numerical simulations. If we apply a Fourier low pass filter to the data, to remove these small details, the resemblance between the simulation results and the field data becomes even more striking. Fig. 9 shows the results of 10 simulations which only differ in their random initial conditions of the sea bed, together with the field data values. In Fig. 9 the left panel shows the sand wave length results. The shown dots are the minimum, mean and maximum per
simulation. The results are taken at the moment the growth in height significantly slows down (around 30 years of growth, left linear clusters of points in figure) and 10 years after that again (right linear clusters). The figure clearly shows the large variation that the random initial disturbances can create. Sand wave lengths between 320 m and 1300 m can occur, with a mean sand wave length between 500 and 600 m. The grey block represents the minimum and maximum value of the sand wave lengths in the field data. The black line and the open squares represent the mean values of the field data and the simulation results respectively. The results for the sand wave height are shown in the right hand panel of Fig. 9. Again there is a large variation, the height can vary between 2.2 m and 21 m, with a mean around 12 m, which coincides with both the shown model results (Fig. 6) while the field data values are in range with the found minimum heights. These results show that variations within sand wave fields can originate from small random disturbances in the initial conditions.
h (m)
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
12 10 8 6 4 2 0 −2 −4 −6 −8
0
1000
2000
3000
4000
33
5000
6000
7000
8000
x (m) 10 8 6 h (m)
4 2 0 −2 −4 −6
h (m)
−8
12 10 8 6 4 2 0 −2 −4 −6 −8
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
0
1000
2000
3000
4000 x (m)
5000
6000
7000
8000
Fig. 7. From different random initial conditions, a different field of sand waves develops. The three plots each show a field of sand waves from a random disturbance with amplitude 0.30 m after approximately 30 years, for asymmetric tidal flow conditions.
8 original filtered
6
h (m)
4 2 0 -2 -4 -6
0
1000
2000
3000
4000
5000
6000
7000
8000
x (m) Fig. 8. Measured sand waves in the North sea. Data taken from the region Kort Verblijf, courtesy RWS Directie Noordzee. The dashed line depicts the same data after Fourier filtering.
34
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
1600
22 20
1400 1200
sand wave height (m)
sand wave length (m)
18
1000 800 600
16 14 12 10 8 6
400 4 200
min
mean max
values per simulation
2
min
mean max
values per simulation
Fig. 9. Statistical results for 10 long domain simulations starting from different randomly disturbed beds. Minimum, mean and maximum sand wave length and height are shown per simulation (open dots) with the mean of these three values (open squares). The results are taken after approximately 30 and 40 years (respectively left and right linear clusters of points for e.g. minimum sand wave length). The grey area represents the field data values (minimum to maximum with the mean represented in the black line).
currents, density gradients, sediment transport) and to describe them within broad classes of problems over different temporal and spatial scales. For example Tonnon et al. (2007) used a complex model to investigate which processes affected an existing sand lump with sand wave dimensions over a time span of 15 years. Due to the complex modelling they were unable to model on a morphological time scale and therefore do not give information on the formation and final state of natural sand waves. In general complex models have two main problems, when aiming for a long timescale and a long domain. Firstly, the computational time is far longer than for an idealized model, such as we use here. Secondly, the numerical diffusion can cause results that are not due to physical processes. Still they can be used for shorter timescales and investigation of specific processes (Borsje et al., in press). The strength of the here presented model is that it is able to model sand waves from initial small bed perturbations up to their full grown shape within a short computational time (order of minutes to tens of minutes). Furthermore the efficient numerical approach and the clear object oriented design processes can be added easily and do not lead to a problematic increase of computational time. It gives the opportunity to investigate the effect of separate processes on sand waves in more detail and to gain insight in leading effects on sand wave characteristics, both in the linear and the non-linear regime. The current aim is not to predict the exact changes in real time sand waves, but to gain insight in the effect of specific processes and correlated minimum and maximum values of e.g. sand wave height and migration.
6. Discussion This solver is already used to study sand waves near the Golden Gate, and showed reasonable agreement with the data when the asymmetry of the tidal current was taken into account (Sterlini et al., 2009). The model is easy to extend with other processes as shown by Van der Meer et al. (2008), Van der Meer et al. (2007), and Roos et al. (2007). Results show that the basic principles are presented correctly in the model and that especially the asymmetry in the tidal current can have a large impact on the sand wave characteristics. Still other processes (e.g. suspended sediment transport and surface waves) can significantly affect the sand waves, and are especially needed to improve the prediction of the final sand wave height. The solver is also used to study river dunes, where similar processes play a role (Hulscher and Dohmen-Janssen, 2005). The model is extended with a parameterization of the flow separation zone and shows good agreement with river dune data (Paarlberg et al., 2009). Due to the efficient numerical algorithm, the model can also be used for operational water management, e.g. for predicting the roughness during high water in rivers (Paarlberg et al., 2010). In shallow seas, the bed is rarely flat, even after dredging. Both larger bed forms (sand banks) and smaller bed forms (mega ripples) can coexist together with sand waves. Although the first can affect the water depth and the second the bed roughness, their influence is not directly included in the current model set up. The used random variation as initial bed is only a rough estimation of what could be found. Still, it can be used as a first indication on the (re)generation time of sand waves. Knaapen and Hulscher (2002) investigated the regeneration of sand waves and assumed an S-shaped growth curve. The results from the here presented model confirm this and could be used to estimate the regeneration of dredged sand waves. As dredging mainly affects the sand wave height and not the length, this process could be studied on a short domain (one wave length). Studying sand waves, one could argue to use a complex or ‘full’ process based models, instead of the idealized model we presented in this paper, as they are developed to include many different processes (e.g. tidal current, wind- and wave-driven
7. Conclusion Offshore sand waves in shallow seas show variability in size, shape and growth rate. Previous theoretical studies of the behaviour of sand waves could not predict this variability. Here, we have presented a solver for sand wave growth and dynamics from their initial small disturbances up to their final equilibrium shape, based on the tidal flow, bed load transport and the feedback mechanism between these processes. The model validation shows good agreement with the analytical results in the linear regime and reasonable accuracy. We use the same model equations as used for the linear stability studies and previous numerical simulations. With the new numerical approach the limitation on the domain size is lifted. The simulation results on large domains show that variations within sand wave fields can originate from small random disturbances in the initial conditions. The solver is able to describe these variations and results show comparable qualitative behaviour with measured field data. Using this solver, other studies show results that agree reasonably well with field and flume data. The model is easy to extent with other physical processes and helps to gain insight in the effect of these processes on sand wave and river dune evolution and final equilibrium. Overall the model is promising as a tool to investigate the effect of physical processes on both sand waves and river dunes. Furthermore, due to the efficient numerical algorithms the model is fast enough for operational use in river dune studies and in large domain sand wave studies. Acknowledgement The research described in this paper is funded by the Technology Foundation STW, applied science division of NWO and the technology program of the Ministry of Economic Affairs (project number TWO.5805: Modelling of spatial and temporal variations in offshore sand waves: process-oriented vs stochastic approach).
J. van den Berg et al. / Continental Shelf Research 37 (2012) 26–35
References Besio, G., Blondeaux, P., Frisina, P., 2003. A note on tidally generated sand waves. Journal of Fluid Dynamics 485, 171–190. Besio, G., Blondeaux, P., Brocchini, M., Vittori, G., 2004. On the modeling of sand wave migration. Journal of Geophysical Research 109, C04018. doi:10.1029/ 2002JC001622. Blondeaux, P., Vittori, G., 2005. Flow and sediment transport induces by tide propagation. 2. The wavy bottom case. Journal of Geophysical Research— Oceans 110, C08003. doi:10.1029/2004JC002545. Bochev, P.B., Hyman, J.M., 2005. Principles of Mimetic Discretizations of Differential Operators. Preprint: LA-UR-05-4244. Borsje, B.W., Roos, P.C., Kranenburg, W.M., Hulscher, S.J.M.H., 2011. Modeling sandwave formation in a numerical shallow water model. In: Shao, X., Wang, Z., Wang, G. (Eds.), RCEM 2011, Seventh IAHR Symposium on RiverCoastal and Estuarine Morphodynamics, Beijing, China. Buijsman, M.C., Ridderinkhof, H., 2008. Long-term evolution of sand waves in the Marsdiep inlet. I. High-resolution observations. Continental Shelf Research 28 (9), 1190–1201. doi:10.1016/j.csr.2007.10.011. Gerkema, T., 2000. A linear stability analysis of tidally generated sand waves. Journal of Fluid Mechanics 417, 303–322. Hulscher, S.J.M.H., 1996. Tidal-induced large-scale regular bed form patterns in a three-dimensional shallow water model. Journal of Geophysical Research 101, 727–744. Hulscher, S.J.M.H., Dohmen-Janssen, C.M., 2005. Introduction to special section on marine sand wave and river dune dynamics. Journal of Geophysical Research 110, F04S01. doi:10.1029/2005JF000404. Hulscher, S.J.M.H., Van den Brink, G.M., 2001. Comparison between predicted and observed sand waves and sand banks in the North Sea. Journal of Geophysical Research 106 (C5), 9327–9338. Hyman, J., Shashkov, M., 1999. The orthogonal decomposition theorems for mimetic finite difference methods. Society for Industrial and Applied Mathematics 36 (3), 788–818. Kelley, C.T., 1995. Iterative Methods for Linear and Nonlinear Equations. SIAM. Knaapen, M.A.F., Hulscher, S.J.M.H., 2002. Regeneration of sand waves after dredging. Coastal Engineering 46, 277–289. Knaapen, M.A.F., Hulscher, S.J.M.H., De Vriend, H.J., 2001. A new type of sea bed waves. Geophysical Research Letters 28 (7), 1323–1326. McCave, I.N., 1971. Sand waves in the North Sea off the coast of Holland. Marine Geology 10 (3), 199–225.
35
Ne´meth, A.A., Hulscher, S.J.M.H., De Vriend, H.J., 2003. Offshore sand wave dynamics, engineering problems and future solutions. Pipeline Gas Journal 230 (4), 67–69. Ne´meth, A.A., Hulscher, S.J.M.H., Van Damme, R.M.J., 2006. Simulating offshore sand waves. Coastal Engineering 53, 265–275. Paarlberg, A.J., Dohmen-Janssen, C.M., Hulscher, S.J.M.H., Termes, P., 2009. Modeling river dune evolution using a parameterization of flow separation. Journal of Geophysical Research—Earth Surface 114, F01014. Paarlberg, A.J., Dohmen-Janssen, C.M., Hulscher, S.J.M.H., Termes, P., Schielen, R., 2010. Modelling the effect of time-dependent river dune evolution on bed roughness and stage. Earth Surface Processes and Landforms 35, 1854–1866. doi:10.1002/esp.2074. Pedlosky, J., 1987. Geophysical Fluid Dynamics. Springer, New York. Roos, P.C., Hulscher, S.J.M.H., Meer, F.M.v.d., Dijk, T.A.G.P.v., Wientjes, I.G.M., Berg, J.v.d., 2007. Grain size sorting over offshore sand waves: observations and modelling. In: Dohmen-Janssen, C.M., Hulscher, S.J.M.H. (Eds.), River, Coastal and Estuarine Morphodynamics: RCEM 2007. Taylor & Francis Group, London. Soulsby, R.L., 1990. Tidal-current boundary layers. In: le Mehaute, B., Hanes, D.M. (Eds.), The Sea, Ocean Engineering Science. John Wiley and Sons, New York, pp. 523–566. Sterlini, F., Hulscher, S.J.H.M., Hanes, D.M., 2009. Simulating and understanding sand wave variation: a case study of the Golden Gate sand waves. Journal of Geophysical Research 114, F02007. doi:10.1029/2008JF000999. Tonnon, P.K., Van-Rijn, L.C., Walstra, D.J.R., 2007. The morphodynamic modelling of tidal sand waves on the shoreface. Coastal Engineering 54, 279–296. Van der Meer, F., Hulscher, S.J.M.H., van den Berg, J., 2007. On the influence of suspended sediment transport on the generation of offshore sand waves. In: B.J.G., et al. (Eds.), Particle Laden Flow: From Geophysical to Kolmogorov Scales. Springer, pp. 29–41. Van der Meer, F., Hulscher, S.J.M.H., Dodd, N., 2008. Modelling surface wave effects on evolved offshore sand waves. In: Souza, Alejandro J., Proctor, R. (Eds.), PECS 2008: Physics of Estuaries and Coastal Seas, Liverpool, UK. Van Rijn, L.C., 1993. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas, vol. III. Aqua Publications, Amsterdam. Vinokur, M., 1974. Conservation equations of gas-dynamics in curvilinear coordinate systems. Journal of Computational Physics 14, 105–125. Viviand, H., 1974. Conservative Forms of Gas Dynamic Equations, vol. 1. Rech. Aerosp., pp. 65–68.