- Email: [email protected]
Shoreline motion in nonlinear shallow water coastal models
Coastal Engineering 56 (2009) 495–505
Contents lists available at ScienceDirect
Coastal Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / c o a s t a l e n g
Shoreline motion in nonlinear shallow water coastal models Riccardo Briganti ⁎, Nicholas Dodd Civil Engineering Dept., Faculty of Engineering, University of Nottingham, Nottingham NG7 2RD, UK
a r t i c l e
i n f o
Article history: Received 16 May 2008 Received in revised form 22 October 2008 Accepted 28 October 2008 Available online 19 December 2008 Keywords: Bore-driven swash Swash zone modeling Shoreline boundary conditions Shock capturing methods
a b s t r a c t Different shoreline boundary conditions for numerical models of the Non-Linear Shallow Water Equations based on Godunov-type schemes are compared. The study focuses on the Peregrine and Williams [Peregrine, D.H., Williams, S.M., 2001. Swash overtopping a truncated plane beach. Journal of Fluid Mechanics 440, 391– 399.] problem of a single bore collapsing on a slope. This is considered the best test to assess performances of the shoreline boundary treatments in terms of all the parameters of interest in swash zone modelling. Emphasis is given to the shoreline trajectory and flow velocity modelling. A mismatch of the velocity at the early stage of the motion is highlighted. Most of the tested techniques perform similarly in terms of maximum run-up, the backwash phase is critical in all cases. Starting from the Brocchini et al. [Brocchini, M., Bernetti, R., Mancinelli, A., Albertini, G., 2001. An efficient solver for nearshore flows based on the WAF method. Coastal Engineering 43(2), 105–129.] shoreline boundary treatment, a simple technique that improves the accuracy of velocity predictions is also developed. A sensitivity analysis of the domain resolution and the threshold value of the water depth that defines a wet cell is also presented. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Numerical modelling of shoreline motion plays a key role in a number of coastal engineering applications, ranging from morphodynamics in the swash zone to wave run-up, coastal flooding and tsunami propagation, significantly influencing the accuracy of the prediction of the phenomenon in hand. In the evaluation of the impact of storms or tsunami events on a coastal region, the estimation of the maximum runup distance is generally the main issue. On the other hand, flow velocity and the duration of a swash event have a primary importance in wavestructure interaction problems, as well as in the prediction of morphodynamical changes. This latter aspect has recently been analyzed in the study of swash zone bedforms (see Dodd et al., 2008). The increasing attention given to these topics has led to sustained efforts by the scientific community to describe swash motion from an analytical and a numerical point of view. The starting point, in both cases, is often the solution of the Non-Linear Shallow Water Equations (NLSWE). Here we focus on the numerical modelling of this set of equations; the reader is referred to Brocchini and Dodd (2008) for a review of numerical techniques. In developing a numerical code it is necessary to use a suitable shoreline boundary condition (SBC) in order to model accurately the swash motion. Several approaches have been proposed in different numerical schemes. Hibberd and Peregrine (1979) used a Lax– Wendroff scheme, implementing an extrapolation technique to
⁎ Corresponding author. E-mail address: [email protected] (R. Briganti). 0378-3839/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.coastaleng.2008.10.008
allow the model to handle movable boundaries. Titov and Synolakis (1995) used a Godunov scheme and proposed a simpler treatment for the shoreline. Extrapolation techniques are also popular in Boussinesqtype models (e.g. Lynett et al., 2002). Other techniques developed in this framework avoid dealing with a wet/dry interface by modelling the beach as porous or slotted (e.g. Kennedy et al., 2000). Among the numerical approaches that have gained popularity in the last decade, finite volume techniques based on upwind Godunovtype schemes (see Toro, 1999, 2001 for a detailed explanation of the related theory) have been extensively used, both in research and engineering practice. These schemes solve the conservative form of the NLSWE and their shock-capturing nature makes them suitable for simulating wave propagation in shallow water, where broken waves behave as bores. One of their main advantages is that they may treat interfaces between wet and dry regions using the Riemann problem analytical solution for the wave propagation speeds. This implies that the treatment of a SBC can be relatively simple and yet accurate. Although different studies tackle this issue, no comparative paper that assesses the different SBCs for these kind of models in the swash exists. In order to fill this gap we have here implemented different SBCs in a state of the art NLSWE solver based on the WAF technique. We also choose a suitable benchmark problem. Previous numerical models used, as a benchmark, analytical solutions of run-up and run-down like those in Carrier and Greenspan (1958) and in Thacker (1981). These solutions do not involve the collapse of a bore on a slope, so they are not suitable to test the performances of SBCs in modelling swash events generated by breaking waves. Furthermore, many previous studies (e.g. Hu et al., 2001; Brocchini et al., 2001) use as a benchmark the analytical solution for a dam break problem on a flat bottom. Although very severe for
496
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Eq. (3) is discretized in space into a set of N control volumes (cells of equal width Δx). Each ith cell (i = 1,N) has its centre located at xi. Hence, at a generic time level n, Wni = W(iΔx,nΔtn). In the following we assume piecewise constant states in each cell. If a generic pair of adjacent cells is considered, they can be referred to as left and right states: (WL,WR) ≡ (Wi,Wi + 1). These define the initial conditions of a subgrid Riemann problem at the inter-cell boundary xi + 1, at each time step. The numerical 2 solution of this system is based on the solution of the relationship n + 1
Wi
n + 12 + 12
Fi
n
= Wi −
i Δt h n + 12 n + 1 Fi + 1 − Fi − 12 + ΔtS: Δx 2 2
is the intercell flux at x = xi
+
1 2
ð4Þ
. We use a Strang approach
(Strang, 1968) to split Eq. (3) into a partial differential equation (PDE) and an ordinary differential equation (ODE): Fig. 1. Sketch of the flow variables and initial conditions of the PW01 problem.
testing wet/dry problems, it allows testing only of the wetting process, without giving any indication of important quantities such as the maximum run-up, shoreline retreat, the elevation and velocity in the backwash and the swash period. Shen and Meyer (1963) described analytically the swash from a uniform bore collapsing on a plane beach. This solution was obtained solving the NLWSE as a boundary value problem. More recently Peregrine and Williams (2001) (PW01, hereinafter) found a simpler and elegant interpretation of this solution, solving an initial value problem for the NLSWE. It follows naturally that this represents a very valuable benchmark for numerical models aimed at the study of swash motions in the presence of breaking waves. It is even more valuable for those applications in which a single event has to be modeled; this may occur, for example, in tsunami simulations. In this paper we model the PW01 problem and compare several SBC treatments. An emphasis is given to the shoreline motion and velocity, an aspect often overlooked in swash zone modelling. The paper is organized as follows: Section 1 is the Introduction, and Section 2 briefly describes the numerical solver. Section 3 reviews the SBCs proposed in the literature and introduces those tested. Section 4 describes and discusses the numerical results. Conclusions are given in Section 5 and an Appendix A gives details on the numerical scheme implemented. 2. The numerical model We solve the NLSWE in conservative form, using the frame of reference shown in Fig. 1. Here x is the slope parallel abscissa, U denotes the depth-averaged slope parallel velocity, h is the local water depth (normal to the slope): Ah AhU + =0 At Ax
AhU + At
A hU 2 +
ð1Þ
2 1 2 gh
Ax
cos α
= − gh sin α
ð2Þ
where g is the gravity acceleration and α is slope angle. In order to derive the scheme we rewrite the NLSWE in a more convenient vector form: AW AFðWÞ + =S At Ax
ð3Þ
T where W=[h,hU]T is the vector of unknowns, F = hU; hU 2 + 12 gh2 cosα T is the flux vector and S = [0,−ghsin α] is the vector of the source terms.
AW AFðWÞ + =0 At Ax
ð5Þ
dW =S dt
ð6Þ
The following sequence of time split operators assures a second order accuracy in time: ODE :
dW =S: dt
PDE :
AW AFðWÞ + =0: At Ax
W⁎
ODE :
dW =S: dt
W⁎⁎
Δt =2
Wn
Y Wn⁎
n + 1
Y
n + 1
Y Wn + 1
Δt
+ 1
n + 1
W⁎⁎
ð7Þ
Δt =2
n+1 are the intermediate solutions at the end of the W ⁎n + 1 and W⁎⁎ first and second step respectively. Note that we use the same sequence of operators used in Hu et al. (2001), while in Brocchini et al. (2001) only two steps are used. The solution of the ODE is carried out at each stage with a RungeKutta fourth order scheme with stepsize control (see Press et al., 1992). The flux in Eq. (4) is evaluated making use of a WAF technique with an HLL approximate Riemann solver. Since the implementation of this technique is not the main focus of the paper, its key features are summarized in the Appendix A.
3. Treatments for the moving shoreline In numerical implementations, the shoreline boundary conditions approximate the theoretically correct ones: hðxs Þ = 0
U ðxs Þ =
dxs dt
ð8Þ
In Godunov-type schemes a shoreline is defined at the interface between a wet and a dry cell. This is, inevitably, an approximation of the actual one, since the zero depth condition is never met. Fig. 2 shows an example of the mismatch between the numerical shoreline and the actual one. It also shows that the cell width is an important factor in shoreline motion modelling. Moreover, the fact that NLSWE are solved in conservative form introduces numerical errors when the water depth becomes very small (Uni + 1 = (Uh)ni + 1 / hni + 1), which requires us to approximate the state of these very shallow cells. A tolerance depth (here hmin) is therefore used to identify those cells. Note that hmin can also be used to locate the shoreline itself. Fig. 3 depicts an advancing shoreline in these schemes. The nominal shoreline can be taken initially at xI + 1 (t =tn), where i =I identifies the 2
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
497
Fig. 2. Conceptual sketch of the numerical and actual shoreline positions at different grid resolutions.
last cell for which hni Nhmin at time tn (top of Fig. 3); I is thus variable in time and follows the shoreline. After integration the solution variables are updated; in particular, cell I(tn) +1 is wetted. However, hI(tnn)++11 bhmin, so that I(tn + 1) =I(tn), and different strategies are used in different SBCs to cope with this, some of which discard water at I(tn) +1. In the final step, in Fig. 3, a further integration yields hI(tnn)++21 Nhmin; therefore at the end of this time step I(tn + 2) =I(tn)+ 1. Note that had we set hI(tnn)++11 =0, we must then have hI(tnn)++21 =0 by a CFL condition. Different SBCs again apply different strategies for where the shoreline is located (e.g. at an interface hL N 0, hR =0, or hL N hmin, hR b hmin, etc.) and for assigning velocities. In the backwash phase, water is progressively drained from the cells until their depth falls below hmin, causing the nominal shoreline to retreat. Note that it is possible that this happens for more than one cell at a time. Different approaches have been used by researchers to implement these steps; a brief literature review is therefore useful. Watson et al. (1992), implemented a WAF scheme based on an exact Riemann solver to study wave propagation, and the geometrical source term is handled using an accelerated reference frame. They used an hmin
as a tolerance value below which the velocity is no longer updated using the momentum equation. Depths such that h b hmin are retained, but the velocity is extrapolated: UI + 1 =UI (see Fig. 3). This extrapolation is similar to that proposed by Titov and Synolakis (1995). Kim et al. (2004) proposed setting UI + 1 = 0. According to Watson et al. (1992), the wave speeds are computed using Eq. (A.1) everywhere except when a dry bed Riemann problem is found (hL N 0 and hR = 0), in which case they are: SL = UL − CL
SR = UL + 2CL
ð9Þ
(see also Toro, 1999). The inter-cell flux is computed everywhere using the WAF expression (Eq. (A.5)) and an exact Riemann solver for the inter-cell problem. Hu et al. (2001) used a MUSCL-Hancock scheme and an HLL approximate Riemann solver. They use a ‘minimum wet depth’, denoted by hmin. The essence of their approach is to adopt a minimum depth to define a wet/dry Riemann problem. Cells with h b hmin are considered dry, i.e. their depth and velocity is set to zero. In the original formulation,
Fig. 3. Numerical shoreline advancing. Cell i = I indicates the position of the last wet cell.
498
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Hu et al. (2001) used the wave estimates given by Eq. (A.1). Note that inter-cell flux at x N xI + 1 is disregarded. 2 Another method was introduced by Brocchini et al. (2001). They used a WAF technique and an exact Riemann solver as in Watson et al. (1992). The SBC treatment relies again on the definition of a tolerance depth hmin to locate the wet/dry interface. The wave speeds are given, again, by Eq. (9). Brocchini et al. (2001) also provided a different estimate of flux at the wet/dry boundary. Unlike the previous methods, which consider all the cells as part of the computational domain, Brocchini et al. (2001) explicitly excludes the dry cells from it. This means that the numerical grid is defined at each time steps by means of a wetting/drying mechanism. Making reference to Fig. 3, Brocchini et al. (2001) consider the cell i =I at time step n to be the last cell to be solved for. Flow variables at i =I + 1 are computed in two different ways, depending on whether UI N 0 (wetting phase) or UI b 0 (drying phase). In the wetting phase the flow variables in cell I + 1 at time step n +1 1 are computed using Eq. (4) in which WIn+ 1 = 0 and the n + WAF flux FI + 12 = 0. The ‘depth history’ of the cell is stored so that the volume of water computed using Eq. (4) is added to that already present. When h N hmin the cell is activated and its velocity computed using Eq. (4). In the drying phase (UI b 0) the cell I + 1 is simply removed from the computational domain, i.e. considered dry when h b hmin. The inter-cell flux at xI + 1 is computed according to: 2
n + 12 + 12
FI
=
n n 1 n 1 F WI + ðc − cI Þ Ffan − F WI 2 2 s
ð10Þ
(b) Water remains in the cell at i = I + 1, which is updated using the continuity equation and the Watson et al. (1992)/Titov and Synolakis (1995) approach to estimate U. The wave speeds are still given by Eq. (A.1). (c) As for (b), but UI + 1 is estimated using the Lynett et al. (2002) linear extrapolation method. (d) As for (b), but the state of the cell I + 1 is estimated by solving the non-conservative form of the NLSWE. (e) This is divided into two sub-options: a) As for (a) but wave speeds estimates, given by Eq. (9), are used. b) As for (b) Eq. (9) is used also when the water depth hI + 1 is computed using the continuity equation and UI + 1 = UI is used. Option 3 Original Brocchini et al. (2001) formulation. Option 4 Brocchini et al. (2001) approach without the ‘depth history’ mechanism and the same wetting/drying technique as in Option 2ea. Different choices are tested for the expression of the flux Ffan: (a) The one proposed in Brocchini et al. (2001), i.e. Eq. (11), h iT ˆ = hˆU; ˆ hˆ U ˆ2 + 1g h ˆ 2 cosα : (b) Ffan = F W 2
Option 5
As for Option 4, except that, instead of the heuristic approximation of the variables at the centre of the rarefaction fan, we use the analytical solution for WI + 1 given in Toro (1999). 2 This, in the case of hL N 0 and hR = 0, reads:
which is a more compact form of Eq. (4.37) in Brocchini et al. (2001). R Δt Here cs = SΔx . Note that Eq. (10) represents a weighted average as in Eq. (A.5), although here no WAF limiter is used. Ffan is the flux in the rarefaction fan. Here SR =Us is the velocity of the shoreline. The flux Ffan is estimated as a function of the left and right depth and velocity, and the ˆ = [hˆ , Uˆ ]T. The authors variables at the centre of the rarefaction fan W 1 ˆ ˆ = 1 ðUI + Us Þ. The provide the heuristic approximation h = 2 hI and U 2 Brocchini et al. (2001) flux estimation reads: n
Ffan =
1 n 1 ˆ n 1 n F WI + F W + F WI + 1 : 4 2 4
8 > < WI
3.1. The tested SBC treatments The aforementioned approaches are included in the numerical solver described in Section 2 and compared to the PW01 solution. Within each approach it is also possible to test the effect of some modifications. Thus, the following alternatives have been tested: Option 1 The original Watson et al. (1992) approach. Option 2 The Hu et al. (2001) approach. After each integration step, the cell with index i=I+1 (Fig. 3), having hbhmin, can be considered dry. Alternatively, water is not removed and the velocity is approximated. The fluxes of cells i=I+k (k=2,…, N−I) are neglected. Different choices for this approximation and for the evaluation of the wave speeds give birth to different alternatives: (a) Cell I + 1 is dry. We use the Hu et al. (2001) choice for wave speeds, i.e. those for a wet/wet problem given in Eq. (A.1).
1 2
=
1 2 ½U + 2CI 9g I
> : U = 1 ½U + 2C I 3 I
ð12Þ
A similar approach has been applied in Pan et al. (2007) within a model based on a Godunov-MUSCL scheme. Again two choices for the flux in the rarefaction fan are tested:
ð11Þ
However, WI + 1 = [0, Us]T, i.e. F(WI + 1) = 0. Finally, Hubbard and Dodd (2002) take an approach more commonly found in fluvial applications. A dry cell (i.e. h b hmin) in which the bed level is lower than the water level in any adjacent wet cells is deemed to be liable to wet in the next integration step. In this case this cell is primed with an initial h =hmin, which allows the interface xI + 1 to be treated as 2 a wet/wet one. This avoids the well-known problem of the wet/dry interface in the Roe solver.
+
h=
(a) The one proposed in Brocchini et al. (2001), i.e. Eq. (11), (b) Ffan = F WI + 1 : 2
Option 6
Hubbard and Dodd (2002) approach.
4. Numerical test results We model the problem described in PW01. Its initial condition is a uniform bore at rest on an infinite slope, see Fig. 1. The solution is given in terms of water depth and velocity: hðx; t Þ =
pffiffiffiffiffiffi 2 1 2 sin α − 2x 4t gA −gt 36gt 2 cosα
ð13Þ
U ðx; t Þ =
2 pffiffiffiffiffiffi 2 4t gA − gt sin α + x : 3t
ð14Þ
2A is the maximum run-up height (see Fig. 1). The bore height is A therefore h0 = cosα . The shoreline trajectory (xs) and velocity (Us) are: xs ðt Þ = 2t
pffiffiffiffiffiffi 1 2 gA − gt sin α 2
pffiffiffiffiffiffi Us ðt Þ = 2 gA − gt sin α
ð15Þ
ð16Þ
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Eq. p ffiffiffi (16) allows computation of the swash event duration Tswash = A 4 pffiffigsinα . The collapse of the bore is always located at x = 0. The numerical tests discussed in this section are carried out using m = tan α = 0.1. The initial bore height is h0 = 1m, and the Courant number Cr = 0.8, in all the tests. A convenient way to express the computational resolution in dimensionless form is to use the number of computational cells used to resolve the theoretical maximum run-up distance, i.e. the ratio Nswash = x s,max / Δx. Here xs,max denotes the maximum run-up predicted by the PW01 solution. We use Nswash = 220, with Δx =0.09 m, which is reasonable for practical computations. In the next section we also confirm that the results are consistent for any Nswash. hmin = 10− 6 m is used in the tests, unless stated otherwise. This value has been chosen since it is the same used in Hu et al. (2001) and within the range of the values suggested in other publications (e.g. Kim et al., 2004). Note that the formula proposed by Toro (2001) yields larger values, in this case hmin = 0.1 tan αΔx = 0.9 × 10− 3 m. It is useful to express the quantities in Eqs. (14)–(17) in dimensionless form. The same scaling arguments in PW01 are used here. Dimensionless quantities pffiffiffiffiffiffiHence x′ =x sin pffiffiffiffiffiffiffiffiffiffi are indicated using a prime. α/A, t V= t sin α g = A, h′ =h cos α/A, U V= U = gA. The maximum run up distance and swash event duration are therefore x′s,max = 2 and T′swash = 4. We note that the PW01 solution, when non-dimensionalised is self similar for any slope and bore height, which implies that all values of A and α are represented by the non-dimensional simulation. First we compare the solution in terms of h′ and U′ resulting from the different tested options. The cell with i =I identifies the numerical shoreline. Fig. 4 shows the prediction of the free surface given by four among the tested options (Option 2a, Option 2ea, Option 3 and Option 4b), chosen as representative of all the techniques used, along with the PM01 solution. The free surface is very well reproduced by the present model (panel (a)) as is the velocity (panel (b)). The differences among the various SBC treatments are concentrated in the tip region. It is evident that the numerical shoreline is located offshore of the analytical one, so
Fig. 4. Numerical results for the PW01 problem at different time steps. (a) h, (b) U. Option 2a (green lines), Option 2ea (red lines), Option 3 (magenta lines) and Option 4 (blue lines) compared with PW01 (black solid lines, run-up phase; black dashed lines, backwash).
499
that the run-up distance is always underestimated. Because of this underestimation of x′s, U′s is always underestimated (overestimated) in the run-up (run-down) phase. In the proximity of the numerical tip the numerical solution does not behave linearly. Fig. 5 is a close-up of the tip region at different times; note the different scales used in the panels. The figure shows four snapshots of h′ (upper panels) and the associated U′ (lower panels). As noted before, most of the differences are seen in the velocity profile. Note that Option 3 shows a significant departure from the analytical curve at the tip during the run-up (t′ = 0.625, t′ = 1.87). Later it will be shown that U′s also shows a noisier behaviour during the backwash. Note that due to the very small steepness of the free surface near the tip, the mismatch between the theoretical and the numerical shoreline is quite significant along the x′ axis. This feature is also present at different resolutions; i.e. it is a ‘self similar’ shape. The different options are now compared in terms of x′s and U′s. First we analyze the performances of Option 1 and Option 2a (see Fig. 6). x′s is normalized with the maximum run-up distance x′s,max. In panel (c) of the figure the phase diagram is also shown, this relates the two variables directly. Results in terms of non-dimensional predicted maximum runup (x′s,max,num) and non-dimensional swash period T′swash are shown in Table 1 for all the options tested. The two options displayed give similar x′s and U′s. U′s is underestimated immediately after the bore collapse, but later this mismatch reduces significantly. It is also interesting to note that the velocity curve is not smooth and U′s grows with a sawtooth-like shape. These behaviours are common to most of the methods tested here. Near the end of the swash cycle the numerical solution is affected by some noise. This problem is emphasized in the phase diagram. Here ‘spikes’ occur in the backwash phase, characterized by a sudden change in U′s while x′s remains constant. During the drying process the nominal shoreline, sometimes, does not retreat for several time steps; the subsequent retreat of the shoreline (often by more than one cell in a time step), causes Us to appear irregular. All the methods tested also underestimate the swash period because of the underprediction of x′s. Fig. 7 shows a comparison of two of the options belonging to the Option 2 family. Only the options that performed best are shown. Option 2c produces a significant underestimation of x′s, and is therefore no longer considered. Also we notice that Option 2d did not produce any benefit compared to Option 2ea and Option 2b, therefore this is not shown as well. Options 2b and 2d, and Option 2eb are not significantly different from Options 2a and 2ea respectively. Furthermore, Options 2b and Options 2eb are found to be impractical. Results generally show that the shoreline trajectory is insensitive to the wave speed estimates used in the SBC. This is evident when Option 2a and Option 2ea are compared directly. If UI + 1 = UI is used (this is done in Options 2b and 2eb) no improvement is seen in the run-up estimation. However we use this approach only in the run-up phase, since it produces very poor results if applied during the backwash as well. This limits the usefulness of this method. Overall all the variants of Option 2 give very similar results. A better performance in the backwash and the capability of preserving the symmetry of the x′s curve make Option 2ea preferable. Therefore, in the following we focus only on this approach within the Option 2 family. Option 3 results differ from those from other methods, mainly in terms of velocity. Fig. 8 shows a comparison between the options that use the basic idea of Brocchini et al. (2001) with modifications; Option 2ea is shown for reference. Note that the shoreline parabola seems to lose symmetry. The maximum of the run-up slightly anticipates the analytical one and the backwash phase lasts longer than the run-up one. To understand this we note that the cells are flooded or dried according to the sign of UI (see Fig. 3). In particular, if UI b 0, the cells with hbhmin are removed from the computational domain. In the other methods the criterion used to dry the ‘shallow’ cells is independent of the value of UI. This sign reversal at x′s =x′s,max induces faster initial shoreline retreat than for the other methods. Later the x′s retreat is slower.
500
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Fig. 5. Numerical results for the PW01 problem. Details at the tip of the swash lens at different time steps. Upper panels: Water depth (h), lower panels: velocity (U). Option 2a (green lines), Option 2ea (red lines), Option 3 (magenta lines) and Option 4 (blue lines) compared with PW01 (black solid lines, run-up phase, black dashed lines, backwash).
The most evident feature of Option 3 is the particularly ‘noisy’ velocity curve in the run-up phase. Again the wetting/drying technique plays an important role. Hence, in Options 4a and 4b the ‘depth history’ mechanism is removed and replaced with the same approach as in Option 2ea. Option 4a (not shown) gives very similar results to the Option 2 ‘family’ for x′s and U′s. However, Table 1 shows poorer results for x′s,max. Therefore it is no longer considered.
As seen in Fig. 8, Option 4b, on the other hand, improves significantly the prediction of the velocity after bore collapse; thereafter differences with the other methods reduce. The x′s curve is also better behaved, (i.e. symmetric and smooth). Note that for Option 3 and 4 U′s decreases with time and then increases to reconnect to the same curve followed by the other methods at around t′ = 0.625. A different approach for the flux computation is that implemented in Option 5a and b. As explained above the two alternatives differ only for
Fig. 6. Comparison between PW01 solution and numerical results. Panel (a): x′s(t)/x′s,max. Panel (b): U′s(t). Panel (c): phase diagram (x′s vs U′s). PW01 solution (black lines) Option 1 (blue lines) and Option 2a (red lines).
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505 Table 1 Maximum run up distance x′s.max,num and swash period T′swash for the tested SBC techniques. Model
x′s.max,num
T′swash
hmin [m]
PW01 Option Option Option Option Option Option Option Option Option Option Option Option Option
2.00 1.767 1.749 1.749 1.749 1.749 1.731 1.767 1.740 1.722 1.749 1.722 1.731 1.731
4.00 3.936 3.873 3.878 3.875 3.875 3.804 3.875 3.900 3.856 3.865 3.91 3.852 3.807
10− 6 10− 6 10− 6 10− 6 10− 6 10− 5 10− 6 10− 6 10− 6 10− 6 10− 6 10− 6 10− 5
1 2a 2b 2d 2ea 2ea 2eb 3 4a 4b 5a 5b 6
Nswash = 220 in all tests.
expressions for the flux in the rarefaction fan. Option 5b gives slightly better results, and so is shown here. However, Fig. 8 shows that Option 5b still performs similarly to Option 2ea. Note that Option 4b and Option 5b differ only by the choice of Ffan. The use of the heuristic approach in Brocchini et al. (2001) gives a higher U′s in the early run-up. Note that in the first time steps the velocity is substantially overestimated. As shown in Fig. 8 this results also in a slightly better prediction of the shoreline trajectory immediately after bore collapse. In applying Option 6 we use hmin =10− 5 m, it being the minimum value that gave stable results in the backwash. Fig. 9 shows the phase diagram for the results of Option 6 compared to Option 2ea using hmin = 10− 5 m (blue line) and hmin = 10− 6 m (green line). There is no significant difference between Option 6 and 2ea but note the difference in U′s in the backwash, due to the larger hmin value for Option 6. The fact that the velocity in the run-up phase is less sensitive to hmin is a feature of all the tested techniques and is a relevant aspect for practical applications. The numerical SBC approximates Eq. (8), which therefore is not satisfied exactly. It is instructive to assess the degree to which Eq. (8) is satisfied in each SBC. In order to do so we introduce the error: SBC =
R xs ðt Þ − 0t Us ðτÞdτ xs;max
ð17Þ
which shows the accumulated error in shoreline modelling. This is plotted in Fig. 10. All results show the underestimation of xs seen in earlier plots (see Fig. 6–9). The monotonic increasing curves show that this grows throughout the swash cycle. Note, however, that for Option 2ea, 4b and 5b a significant portion of the mid to late run-up phase
501
shows no increase in error (i.e. good modelling), after an initially rapidly increasing error. The backwash phase marks the start of further underprediction: see Fig. 10. Of further note are the significantly poorer performance of Option 3 and the corresponding better performance of Option 4b, for which the improved prediction of U′s immediately after bore collapse is apparent. All the alternatives of the Option 2 ‘family’ show a very similar magnitude of εSBC as do Option 5b and 6 (these latter two are not shown here). Thanks to the lower error in the velocity prediction Option 4b performs better in terms of SBC balance. However, the results in the backwash do not differ significantly from those of the other methods. 4.1. On grid spacing and hmin Both the domain resolution and the value of hmin are critical aspects in modelling of the swash zone. In Fig. 11 x′s and U′s are plotted for different mesh sizes (i.e. for different Nswash =xs,max / Δx) for Option 2ea (Panels (a) and (b)), chosen as representative of the ones derived from Hu et al. (2001), and for Option 4b (Panels (c) and (d)); which improves the velocity description at bore collapse. Note that the mismatch of the velocity at the early run-up stages reduces with the increasing resolution (Panels (b) and (d)), Option 4b improves the prediction of U′s at any resolution in the run-up phase. The initial U′s is the same for any resolution. This is because the solution of the Riemann problem at t = 0 is the same at any resolution. Since the highest gradients occur at the beginning of the event, the resolution of the steep solution is crucial. Hence, the decreasing cell dimension reduces the time interval during which the initial conditions influence the accuracy. The influence of hmin and Nswash is here studied by analyzing the errors in the maximum run-up prediction and in U′s as a function of the ratio h′min /Δx′=hmin /Δx tan α) and of Nswash. For the run-up we define: run−up =
xs;max;num − xs;max × 100 xs;max
ð18Þ
To evaluate the accuracy of U′s for the whole of the swash cycle we use the root mean square error (RMSE) defined as:
U sV
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uN ts uX ðU Vns −U Vs ðtn ÞÞ2 =t : Nts n=1
ð19Þ
Nts is the number of time steps in the swash cycle, U′ns is the numerical prediction of U′s at the n-th time step and U′s(tn) is the corresponding analytical value. In Fig. 12 contours of εrun-up (upper panel) and εU′s (lower panel) are shown as functions of h′min / Δx′ and Nswash for Option
Fig. 7. Comparison between PW01 solution and numerical results. Panel (a): x′s(t) / x′s,max. Panel (b): U′s(t). Panel (c): phase diagram (x′s(t) / x′s,max vs U′s). PW01 solution (black lines) Option 2a (blue lines) and Option 2ea (red lines).
502
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Fig. 8. Comparison between PW01 solution and numerical results. Panels (a) and (d): x′s(t)/x′s,max. Panels (b) and (e): U′s(t). Panels (c) and (f): phase diagram (x′s(t)/x′s,max vs U′s). Panels (a), (b), (c): PW01 solution (black lines) Option 2ea (blue lines) and Option 3 (red lines). Panels (d), (e), (f): PW01 solution (black lines) Option 4b (blue lines) and Option 5b (red lines).
2ea. In the plots h′min / Δx′ varies from 1 × 10− 1 to 1 × 10− 5, whilst the maximum Nswash is 1550. Tests performed with smaller h′min / Δx′ and larger Nswash resulted in poor performances in the backwash, so they are not shown. The first conclusion that can be drawn from the figure is that εrun-up monotonically decreases as either h′min / Δx′ decreases or Nswash increases. For εU′s the reduction of h′min / Δx′ is not always beneficial. For h′min / Δx′ b 1 × 10− 3, εU′s does not necessarily indicate convergence for
Fig. 9. Phase diagram (x′s(t)/x′s,max vs U′s(t)) for Option 6 using hmin = 10− 5 m, (red line) compared with Option 2ea using hmin = 10− 6 m (green line) and hmin = 10− 5 m (blue line). PW01 solution is also plotted (black line).
increasing Nswash. This is due to the increasing error in the backwash phase for very small h′min / Δx′. The contour plot of εrun-up shows that, for h′min / Δx′ N 1 × 10− 3, decreasing h′min / Δx′ for constant Nswash produces a greater benefit than increasing resolution alone. For h′min / Δx′ b 1 × 10− 3 resolution becomes increasingly more important in reducing εrun-up. Obviously, while decreasing h′min is not computationally expensive increasing the resolution is. Therefore, ideally we choose a suitable value of h′min / Δx′ that preserves accuracy whilst saving on resolution.
Fig. 10. SBC balance for Options 2ea, 3, 4b and 5b.
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
503
Fig. 11. x′s(t) and U′s(t) for different values of Nswash. Panels (a) and (b): Option 2ea. Nswash = 125 (magenta line), Nswash = 250 (green line), Nswash = 500 (blue line), Nswash = 1000 (red line). Panels (c) and (d): Option 4b. Nswash = 125 (magenta line), Nswash = 250 (green line), Nswash = 500 (blue line), Nswash = 1000 (red line). PW01 solution is always plotted with a black line.
We conclude that h′min / Δx′ = 1 × 10− 3 is an optimal value since it guarantees that both εrun-up and εU′s indicate convergence for increasing Nswash; furthermore smaller values of the ratio produce
moderate benefits for the run-up but a poorer performance for the velocity of the backwash, causing the increase in εU′s. Option 4b produces very similar results to the ones shown in Fig. 12, with εrun-up
Fig. 12. Upper panel: contour plot of εrun-up for Option 2ea for different values of h′min / Δx′ and Nswash. Lower panel: contour plot of εU′s for Option 2ea for different values of h′min / Δx′ and Nswash. Note that max(U′s) = 2.
504
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
slightly larger at all resolutions. However, we found that εU′s is larger when h′min / Δx′ N 1 × 10− 3 as a result of worse performance in the backwash. Note that h′min /Δx′ =hmin/(tan αΔx) = 1 × 10− 3 is smaller than the ratio suggested by Watson et al. (1992) (i.e. h′min /(tan αΔx) = 1 × 10− 1). Furthermore, note that transformation from an x-axis running along the bed to a horizontal axis, leaves the ratio hmin /Δx unchanged. 5. Conclusions A comparative study of different treatments of the shoreline boundary conditions in Godunov-type schemes has been presented in this paper. We focus on swash motion and use the PW01 problem as a demanding benchmark for a state of the art WAF solver. Most of the tested SBCs show comparable performances in terms of x′s and U′s. The modelling of the bulk of the flow, on the other hand, is not significantly affected by the choice of the SBC. A critical aspect when modelling the shoreline motion with the solvers used here is the accuracy of the velocity immediately after bore collapse. Most of the techniques tested exhibit a very similar behaviour, with the velocity at the tip increasing in time and reconnecting with the theoretical curve at a later stage. The duration of this ‘growing curve’ depends strongly on the numerical resolution. A practical solution to this problem is the use of a modified version of the Brocchini et al. (2001) technique, referred to here as Option 4b. This improves the prediction of the velocity at early stage and gives better results than the original Brocchini et al. (2001) formulation. Among the approaches derived from the Hu et al. (2001) ‘minimum wet depth’ technique, Option 2ea performs slightly better and proves to be robust at any resolution both in the run-up and the backwash phase. The sensitivity analysis carried out shows that the improvement in U′s in Option 4b (compared to Option 2ea) is more significant at resolutions Nswash b 500. At this resolution, the initial velocity is still better modeled by Option 4b, but the maximum run-up is slightly better predicted by Option 2ea. Hence we suggest to use either one or the other according to the grid size required in the problem under consideration. In all the methods tested the prediction of U′s in the backwash phase seems to be a demanding aspect. It is shown that the drying process is usually not smooth, especially at lower resolutions, affecting, in turn, the U′s prediction. The fact that in the PW01 solution the free surface is always tangential to the bottom at the tip, implies that in the backwash a thin film of water is present, causing the numerical round-off errors to be important in the solution. This effect is less relevant when friction is included. We also analysed the relative role of hmin and Nswash in the accuracy of the maximum run-up and U′s prediction. The outcome of this analysis is that h′min / Δx′ ≈ 1 × 10− 3 can be used as a predimensioning rule for testing a code against the PW01 solution. When Nswash = 1500 and the suggested ratio is used, εrun-up = 5% a value that is acceptable for engineering purposes. Slightly better performances can be obtained by reducing h′min / Δx′, but with possibly poorer performances for U′s. Although this resolution may seem quite high, it is important to remember that the (a) methods presented are relatively computationally inexpensive, and parallel computing can reduce the computation time for two dimensional simulations; and (b) for practical engineering problems this kind of resolution is unlikely to be necessary because of the presence of friction, and pure dam-break problem such as that of the PW01 solution are rarely encountered in practical applications. Other numerical models, such as those based on the method of characteristics, achieve higher accuracy than those used here. However, those methods are usually significantly more demanding to implement, especially in two dimensions, and are computationally more expensive, and so they rarely find practical engineering application. Much of the mismatch discussed here stems from the fact that the analytical solution is tangential to the bottom at the shoreline, a con-
dition that is not met by the numerical prediction. This makes the PW01 solution a very demanding test for swash zone applications. When friction is present in the system this condition no longer pertains, and in turn results are expected to be more accurate. Similarly, for other swash events (see e.g. Hibberd and Peregrine, 1979) bore height will diminish when encountering the shore, thus alleviating the problem with estimating the initial U′s that we illustrated here. Nevertheless, good modelling of the PW01 solution will certainly produce good performances for swash flows with friction or in less demanding cases. The implications of the above results for both swash zone morphodynamic modelling and wave overtopping are currently being studied and represent the natural continuation of this work. Acknowledgements This study was supported by the EPSRC through the research project “Experimental and Numerical Modelling Study of Swash Zone Hydrodynamics and Sediment Transport” (EP/EO10407/1). Their support is gratefully acknowledged. David Kelly of the University of Nottingham, Civil Engineering Department is also acknowledged for the many productive discussions during this research. Appendix A. The WAF technique The solution of the homogeneous problem is carried out using a WAF approach for the computation of the inter-cell flux appearing in Eq. (4). This requires the solution of a subgrid Riemann problem. At this purpose we chose to use a HLL (Harten et al., 1983) approximate Riemann solver. This is widely used and performs well in cases of critical rarefactions and dry fronts. The solver assumes the same wave structure for any inter-cell problem solution, i.e. two waves that separate three regions of constant states. The propagation speed estimates for the two waves (SL, SR) are (Toro, 1999): SL = UL − CL qL
SR = UR + CR qR
ðA:1Þ
where the subscripts L and R indicate the left and right state atffi each pffiffiffiffiffi inter-cell boundary. C is the wave celerity in shallow water ( gh) and qK (K = L,R) is given by:
qK =
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8v 3 u 2 ⁎ ⁎ > u > u < t1 4 h + hK h 5 h2K
2
> > :
1
if if
h⁎ N hK ð22Þ
ðA:2Þ
⁎
h V hK
where h⁎ is an estimate for the solution in the so called ’star region’, i.e. the intermediate region between the two waves. This estimate is provided by a Two-Rarefaction approximate Riemann solver: ⁎
h =
2 1 1 1 ðCL + CR Þ + ðUL − UR Þ : g 2 4
ðA:3Þ
Following Toro (1999) is possible to obtain the HLL estimate of the flux in the star region: = FHLL ⁎
SR FL − SL FR + SR SL ðWR − WL Þ : SR − SL
ðA:4Þ
Once the flux in the star region is known, it is possible to compute the inter-cell flux according to the WAF technique. In the TVD (total variation diminishing) version of the WAF, its expression is given by the relationship: n + 12 + 12
Fi
=
2 n 1 n 1X ðlÞ F Wi + F W i + 1 + signðcl Þl ΔFi + 1 2 2 2 l=1
ðA:5Þ
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
where F(Wni ) = FnL is the flux vector computed considering the left state and F(Wni + 1) = FnR is the flux related to the right state. Note that they are both computed at the time level n. ϕl (l = 1,2) is a WAF limiter function. In this study we use the limiter SUPERA (see Toro, 1999). ðlÞ Finally, ΔFi + 1 is the flux jump across each zone identified in the 2 wave structure of the solution. In Eq. (A.5), cl is the Courant number defined for each l-th wave of speed Sl, i.e.: cl =
Sl Δt Δx
ðA:6Þ
References Brocchini, M., Dodd, N., 2008. Nonlinear shallow water equations modeling for coastal engineering. Journal of Waterway, Port, Coastal and Ocean Engineering 134 (2) (6), 104–120. Brocchini, M., Bernetti, R., Mancinelli, A., Albertini, G., 2001. An efficient solver for nearshore flows based on the WAF method. Coastal Engineering 43 (2), 105–129. Carrier, G., Greenspan, H., 1958. Water waves of finite amplitude on a sloping beach. J. Fluid Mech. 4, 97–109. Dodd, N., Stoker, A.M., Calvete, D., Sriariyawat, A., 2008. On the formation of beach cusps. Journal of Fluid Mechanics 597, 145–169. Harten, A., Lax, P., van Leer, B., 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25 (1), 35–61. Hibberd, S., Peregrine, D.H., 1979. Surf and run-up on a beach. Journal of Fluid Mechanics 95, 323–345. Hubbard, M.E., Dodd, N., 2002. A 2D numerical model of wave run-up and overtopping. Coastal Engineering 47 (1), 1–26. Hu, K., Mingham, C.G., Causon, D.M., 2001. Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations. Coastal Engineering 41 (4), 433–465.
505
Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000. Boussinesq modeling of wave transformation, breaking and runup. I: 1D. Journal of Waterways Port Coastal and Ocean Engineering-ASCE 126, 39–47. Kim, A., Cho, Y.S., Kim, W.G., 2004. Weighted averaged flux-type scheme for shallowwater equations with fractional step methods. Journal of Engineering Mechanics 130 (2), 152–160. Lynett, P., Wu, T.-R., Liu, P.L.-F., 2002. Modeling wave runup with depth-integrated equations. Coastal Engineering 46 (2), 89–107. Pan, C.-H., Lin, B.-Y., Mao, X.-Z., 2007. Case study: numerical modeling of the tidal bore on Qiantang River, China. Journal of Hydraulic Engineering 133 (2), 130–138. Peregrine, D.H., Williams, S.M., 2001. Swash overtopping a truncated plane beach. Journal of Fluid Mechanics 440, 391–399. Press, W., Flannery, B., Teukolsky, S., Vetterling, W., 1992. Numerical Recipes: The Art of Scientific Computing. Cambridge Univ. Press, New York. 992 pp. Shen, M.C., Meyer, R.E., 1963. Climb of a bore on a beach. Part 3. Runup. Journal of Fluid Mechanics 16, 113–125. Strang, G., 1968. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 124, 506–517. Thacker, W.C., 1981. Some exact solutions to the nonlinear shallow water wave equations. Journal of Fluid Mechanics 107, 499–508. Titov, V.V., Synolakis, C.E., 1995. Modeling of breaking and nonbreaking long wave evolution and runup using VTCS-2. Journal of Waterway, Port, Coastal and Ocean Engineering 121 (6), 308–316. Toro, E.F., 1999. Riemann solvers and numerical methods for fluid dynamics. Springer, New York. Toro, E.F., 2001. Shock-capturing methods for free-surface shallow flows. Springer, New York. Watson, G., Peregrine, D.H., Toro, E.F., 1992. Numerical solution of the shallow water equations on a beach using the weighted average flux method. In: Hirsch, P., Kordulla (Eds.), Proceedings of the First European Computational Fluid Dynamics Conference. Elsevier, Brussels, pp. 495–502. September 7th to 11th, 1992.
Contents lists available at ScienceDirect
Coastal Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / c o a s t a l e n g
Shoreline motion in nonlinear shallow water coastal models Riccardo Briganti ⁎, Nicholas Dodd Civil Engineering Dept., Faculty of Engineering, University of Nottingham, Nottingham NG7 2RD, UK
a r t i c l e
i n f o
Article history: Received 16 May 2008 Received in revised form 22 October 2008 Accepted 28 October 2008 Available online 19 December 2008 Keywords: Bore-driven swash Swash zone modeling Shoreline boundary conditions Shock capturing methods
a b s t r a c t Different shoreline boundary conditions for numerical models of the Non-Linear Shallow Water Equations based on Godunov-type schemes are compared. The study focuses on the Peregrine and Williams [Peregrine, D.H., Williams, S.M., 2001. Swash overtopping a truncated plane beach. Journal of Fluid Mechanics 440, 391– 399.] problem of a single bore collapsing on a slope. This is considered the best test to assess performances of the shoreline boundary treatments in terms of all the parameters of interest in swash zone modelling. Emphasis is given to the shoreline trajectory and flow velocity modelling. A mismatch of the velocity at the early stage of the motion is highlighted. Most of the tested techniques perform similarly in terms of maximum run-up, the backwash phase is critical in all cases. Starting from the Brocchini et al. [Brocchini, M., Bernetti, R., Mancinelli, A., Albertini, G., 2001. An efficient solver for nearshore flows based on the WAF method. Coastal Engineering 43(2), 105–129.] shoreline boundary treatment, a simple technique that improves the accuracy of velocity predictions is also developed. A sensitivity analysis of the domain resolution and the threshold value of the water depth that defines a wet cell is also presented. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Numerical modelling of shoreline motion plays a key role in a number of coastal engineering applications, ranging from morphodynamics in the swash zone to wave run-up, coastal flooding and tsunami propagation, significantly influencing the accuracy of the prediction of the phenomenon in hand. In the evaluation of the impact of storms or tsunami events on a coastal region, the estimation of the maximum runup distance is generally the main issue. On the other hand, flow velocity and the duration of a swash event have a primary importance in wavestructure interaction problems, as well as in the prediction of morphodynamical changes. This latter aspect has recently been analyzed in the study of swash zone bedforms (see Dodd et al., 2008). The increasing attention given to these topics has led to sustained efforts by the scientific community to describe swash motion from an analytical and a numerical point of view. The starting point, in both cases, is often the solution of the Non-Linear Shallow Water Equations (NLSWE). Here we focus on the numerical modelling of this set of equations; the reader is referred to Brocchini and Dodd (2008) for a review of numerical techniques. In developing a numerical code it is necessary to use a suitable shoreline boundary condition (SBC) in order to model accurately the swash motion. Several approaches have been proposed in different numerical schemes. Hibberd and Peregrine (1979) used a Lax– Wendroff scheme, implementing an extrapolation technique to
⁎ Corresponding author. E-mail address: [email protected] (R. Briganti). 0378-3839/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.coastaleng.2008.10.008
allow the model to handle movable boundaries. Titov and Synolakis (1995) used a Godunov scheme and proposed a simpler treatment for the shoreline. Extrapolation techniques are also popular in Boussinesqtype models (e.g. Lynett et al., 2002). Other techniques developed in this framework avoid dealing with a wet/dry interface by modelling the beach as porous or slotted (e.g. Kennedy et al., 2000). Among the numerical approaches that have gained popularity in the last decade, finite volume techniques based on upwind Godunovtype schemes (see Toro, 1999, 2001 for a detailed explanation of the related theory) have been extensively used, both in research and engineering practice. These schemes solve the conservative form of the NLSWE and their shock-capturing nature makes them suitable for simulating wave propagation in shallow water, where broken waves behave as bores. One of their main advantages is that they may treat interfaces between wet and dry regions using the Riemann problem analytical solution for the wave propagation speeds. This implies that the treatment of a SBC can be relatively simple and yet accurate. Although different studies tackle this issue, no comparative paper that assesses the different SBCs for these kind of models in the swash exists. In order to fill this gap we have here implemented different SBCs in a state of the art NLSWE solver based on the WAF technique. We also choose a suitable benchmark problem. Previous numerical models used, as a benchmark, analytical solutions of run-up and run-down like those in Carrier and Greenspan (1958) and in Thacker (1981). These solutions do not involve the collapse of a bore on a slope, so they are not suitable to test the performances of SBCs in modelling swash events generated by breaking waves. Furthermore, many previous studies (e.g. Hu et al., 2001; Brocchini et al., 2001) use as a benchmark the analytical solution for a dam break problem on a flat bottom. Although very severe for
496
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Eq. (3) is discretized in space into a set of N control volumes (cells of equal width Δx). Each ith cell (i = 1,N) has its centre located at xi. Hence, at a generic time level n, Wni = W(iΔx,nΔtn). In the following we assume piecewise constant states in each cell. If a generic pair of adjacent cells is considered, they can be referred to as left and right states: (WL,WR) ≡ (Wi,Wi + 1). These define the initial conditions of a subgrid Riemann problem at the inter-cell boundary xi + 1, at each time step. The numerical 2 solution of this system is based on the solution of the relationship n + 1
Wi
n + 12 + 12
Fi
n
= Wi −
i Δt h n + 12 n + 1 Fi + 1 − Fi − 12 + ΔtS: Δx 2 2
is the intercell flux at x = xi
+
1 2
ð4Þ
. We use a Strang approach
(Strang, 1968) to split Eq. (3) into a partial differential equation (PDE) and an ordinary differential equation (ODE): Fig. 1. Sketch of the flow variables and initial conditions of the PW01 problem.
testing wet/dry problems, it allows testing only of the wetting process, without giving any indication of important quantities such as the maximum run-up, shoreline retreat, the elevation and velocity in the backwash and the swash period. Shen and Meyer (1963) described analytically the swash from a uniform bore collapsing on a plane beach. This solution was obtained solving the NLWSE as a boundary value problem. More recently Peregrine and Williams (2001) (PW01, hereinafter) found a simpler and elegant interpretation of this solution, solving an initial value problem for the NLSWE. It follows naturally that this represents a very valuable benchmark for numerical models aimed at the study of swash motions in the presence of breaking waves. It is even more valuable for those applications in which a single event has to be modeled; this may occur, for example, in tsunami simulations. In this paper we model the PW01 problem and compare several SBC treatments. An emphasis is given to the shoreline motion and velocity, an aspect often overlooked in swash zone modelling. The paper is organized as follows: Section 1 is the Introduction, and Section 2 briefly describes the numerical solver. Section 3 reviews the SBCs proposed in the literature and introduces those tested. Section 4 describes and discusses the numerical results. Conclusions are given in Section 5 and an Appendix A gives details on the numerical scheme implemented. 2. The numerical model We solve the NLSWE in conservative form, using the frame of reference shown in Fig. 1. Here x is the slope parallel abscissa, U denotes the depth-averaged slope parallel velocity, h is the local water depth (normal to the slope): Ah AhU + =0 At Ax
AhU + At
A hU 2 +
ð1Þ
2 1 2 gh
Ax
cos α
= − gh sin α
ð2Þ
where g is the gravity acceleration and α is slope angle. In order to derive the scheme we rewrite the NLSWE in a more convenient vector form: AW AFðWÞ + =S At Ax
ð3Þ
T where W=[h,hU]T is the vector of unknowns, F = hU; hU 2 + 12 gh2 cosα T is the flux vector and S = [0,−ghsin α] is the vector of the source terms.
AW AFðWÞ + =0 At Ax
ð5Þ
dW =S dt
ð6Þ
The following sequence of time split operators assures a second order accuracy in time: ODE :
dW =S: dt
PDE :
AW AFðWÞ + =0: At Ax
W⁎
ODE :
dW =S: dt
W⁎⁎
Δt =2
Wn
Y Wn⁎
n + 1
Y
n + 1
Y Wn + 1
Δt
+ 1
n + 1
W⁎⁎
ð7Þ
Δt =2
n+1 are the intermediate solutions at the end of the W ⁎n + 1 and W⁎⁎ first and second step respectively. Note that we use the same sequence of operators used in Hu et al. (2001), while in Brocchini et al. (2001) only two steps are used. The solution of the ODE is carried out at each stage with a RungeKutta fourth order scheme with stepsize control (see Press et al., 1992). The flux in Eq. (4) is evaluated making use of a WAF technique with an HLL approximate Riemann solver. Since the implementation of this technique is not the main focus of the paper, its key features are summarized in the Appendix A.
3. Treatments for the moving shoreline In numerical implementations, the shoreline boundary conditions approximate the theoretically correct ones: hðxs Þ = 0
U ðxs Þ =
dxs dt
ð8Þ
In Godunov-type schemes a shoreline is defined at the interface between a wet and a dry cell. This is, inevitably, an approximation of the actual one, since the zero depth condition is never met. Fig. 2 shows an example of the mismatch between the numerical shoreline and the actual one. It also shows that the cell width is an important factor in shoreline motion modelling. Moreover, the fact that NLSWE are solved in conservative form introduces numerical errors when the water depth becomes very small (Uni + 1 = (Uh)ni + 1 / hni + 1), which requires us to approximate the state of these very shallow cells. A tolerance depth (here hmin) is therefore used to identify those cells. Note that hmin can also be used to locate the shoreline itself. Fig. 3 depicts an advancing shoreline in these schemes. The nominal shoreline can be taken initially at xI + 1 (t =tn), where i =I identifies the 2
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
497
Fig. 2. Conceptual sketch of the numerical and actual shoreline positions at different grid resolutions.
last cell for which hni Nhmin at time tn (top of Fig. 3); I is thus variable in time and follows the shoreline. After integration the solution variables are updated; in particular, cell I(tn) +1 is wetted. However, hI(tnn)++11 bhmin, so that I(tn + 1) =I(tn), and different strategies are used in different SBCs to cope with this, some of which discard water at I(tn) +1. In the final step, in Fig. 3, a further integration yields hI(tnn)++21 Nhmin; therefore at the end of this time step I(tn + 2) =I(tn)+ 1. Note that had we set hI(tnn)++11 =0, we must then have hI(tnn)++21 =0 by a CFL condition. Different SBCs again apply different strategies for where the shoreline is located (e.g. at an interface hL N 0, hR =0, or hL N hmin, hR b hmin, etc.) and for assigning velocities. In the backwash phase, water is progressively drained from the cells until their depth falls below hmin, causing the nominal shoreline to retreat. Note that it is possible that this happens for more than one cell at a time. Different approaches have been used by researchers to implement these steps; a brief literature review is therefore useful. Watson et al. (1992), implemented a WAF scheme based on an exact Riemann solver to study wave propagation, and the geometrical source term is handled using an accelerated reference frame. They used an hmin
as a tolerance value below which the velocity is no longer updated using the momentum equation. Depths such that h b hmin are retained, but the velocity is extrapolated: UI + 1 =UI (see Fig. 3). This extrapolation is similar to that proposed by Titov and Synolakis (1995). Kim et al. (2004) proposed setting UI + 1 = 0. According to Watson et al. (1992), the wave speeds are computed using Eq. (A.1) everywhere except when a dry bed Riemann problem is found (hL N 0 and hR = 0), in which case they are: SL = UL − CL
SR = UL + 2CL
ð9Þ
(see also Toro, 1999). The inter-cell flux is computed everywhere using the WAF expression (Eq. (A.5)) and an exact Riemann solver for the inter-cell problem. Hu et al. (2001) used a MUSCL-Hancock scheme and an HLL approximate Riemann solver. They use a ‘minimum wet depth’, denoted by hmin. The essence of their approach is to adopt a minimum depth to define a wet/dry Riemann problem. Cells with h b hmin are considered dry, i.e. their depth and velocity is set to zero. In the original formulation,
Fig. 3. Numerical shoreline advancing. Cell i = I indicates the position of the last wet cell.
498
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Hu et al. (2001) used the wave estimates given by Eq. (A.1). Note that inter-cell flux at x N xI + 1 is disregarded. 2 Another method was introduced by Brocchini et al. (2001). They used a WAF technique and an exact Riemann solver as in Watson et al. (1992). The SBC treatment relies again on the definition of a tolerance depth hmin to locate the wet/dry interface. The wave speeds are given, again, by Eq. (9). Brocchini et al. (2001) also provided a different estimate of flux at the wet/dry boundary. Unlike the previous methods, which consider all the cells as part of the computational domain, Brocchini et al. (2001) explicitly excludes the dry cells from it. This means that the numerical grid is defined at each time steps by means of a wetting/drying mechanism. Making reference to Fig. 3, Brocchini et al. (2001) consider the cell i =I at time step n to be the last cell to be solved for. Flow variables at i =I + 1 are computed in two different ways, depending on whether UI N 0 (wetting phase) or UI b 0 (drying phase). In the wetting phase the flow variables in cell I + 1 at time step n +1 1 are computed using Eq. (4) in which WIn+ 1 = 0 and the n + WAF flux FI + 12 = 0. The ‘depth history’ of the cell is stored so that the volume of water computed using Eq. (4) is added to that already present. When h N hmin the cell is activated and its velocity computed using Eq. (4). In the drying phase (UI b 0) the cell I + 1 is simply removed from the computational domain, i.e. considered dry when h b hmin. The inter-cell flux at xI + 1 is computed according to: 2
n + 12 + 12
FI
=
n n 1 n 1 F WI + ðc − cI Þ Ffan − F WI 2 2 s
ð10Þ
(b) Water remains in the cell at i = I + 1, which is updated using the continuity equation and the Watson et al. (1992)/Titov and Synolakis (1995) approach to estimate U. The wave speeds are still given by Eq. (A.1). (c) As for (b), but UI + 1 is estimated using the Lynett et al. (2002) linear extrapolation method. (d) As for (b), but the state of the cell I + 1 is estimated by solving the non-conservative form of the NLSWE. (e) This is divided into two sub-options: a) As for (a) but wave speeds estimates, given by Eq. (9), are used. b) As for (b) Eq. (9) is used also when the water depth hI + 1 is computed using the continuity equation and UI + 1 = UI is used. Option 3 Original Brocchini et al. (2001) formulation. Option 4 Brocchini et al. (2001) approach without the ‘depth history’ mechanism and the same wetting/drying technique as in Option 2ea. Different choices are tested for the expression of the flux Ffan: (a) The one proposed in Brocchini et al. (2001), i.e. Eq. (11), h iT ˆ = hˆU; ˆ hˆ U ˆ2 + 1g h ˆ 2 cosα : (b) Ffan = F W 2
Option 5
As for Option 4, except that, instead of the heuristic approximation of the variables at the centre of the rarefaction fan, we use the analytical solution for WI + 1 given in Toro (1999). 2 This, in the case of hL N 0 and hR = 0, reads:
which is a more compact form of Eq. (4.37) in Brocchini et al. (2001). R Δt Here cs = SΔx . Note that Eq. (10) represents a weighted average as in Eq. (A.5), although here no WAF limiter is used. Ffan is the flux in the rarefaction fan. Here SR =Us is the velocity of the shoreline. The flux Ffan is estimated as a function of the left and right depth and velocity, and the ˆ = [hˆ , Uˆ ]T. The authors variables at the centre of the rarefaction fan W 1 ˆ ˆ = 1 ðUI + Us Þ. The provide the heuristic approximation h = 2 hI and U 2 Brocchini et al. (2001) flux estimation reads: n
Ffan =
1 n 1 ˆ n 1 n F WI + F W + F WI + 1 : 4 2 4
8 > < WI
3.1. The tested SBC treatments The aforementioned approaches are included in the numerical solver described in Section 2 and compared to the PW01 solution. Within each approach it is also possible to test the effect of some modifications. Thus, the following alternatives have been tested: Option 1 The original Watson et al. (1992) approach. Option 2 The Hu et al. (2001) approach. After each integration step, the cell with index i=I+1 (Fig. 3), having hbhmin, can be considered dry. Alternatively, water is not removed and the velocity is approximated. The fluxes of cells i=I+k (k=2,…, N−I) are neglected. Different choices for this approximation and for the evaluation of the wave speeds give birth to different alternatives: (a) Cell I + 1 is dry. We use the Hu et al. (2001) choice for wave speeds, i.e. those for a wet/wet problem given in Eq. (A.1).
1 2
=
1 2 ½U + 2CI 9g I
> : U = 1 ½U + 2C I 3 I
ð12Þ
A similar approach has been applied in Pan et al. (2007) within a model based on a Godunov-MUSCL scheme. Again two choices for the flux in the rarefaction fan are tested:
ð11Þ
However, WI + 1 = [0, Us]T, i.e. F(WI + 1) = 0. Finally, Hubbard and Dodd (2002) take an approach more commonly found in fluvial applications. A dry cell (i.e. h b hmin) in which the bed level is lower than the water level in any adjacent wet cells is deemed to be liable to wet in the next integration step. In this case this cell is primed with an initial h =hmin, which allows the interface xI + 1 to be treated as 2 a wet/wet one. This avoids the well-known problem of the wet/dry interface in the Roe solver.
+
h=
(a) The one proposed in Brocchini et al. (2001), i.e. Eq. (11), (b) Ffan = F WI + 1 : 2
Option 6
Hubbard and Dodd (2002) approach.
4. Numerical test results We model the problem described in PW01. Its initial condition is a uniform bore at rest on an infinite slope, see Fig. 1. The solution is given in terms of water depth and velocity: hðx; t Þ =
pffiffiffiffiffiffi 2 1 2 sin α − 2x 4t gA −gt 36gt 2 cosα
ð13Þ
U ðx; t Þ =
2 pffiffiffiffiffiffi 2 4t gA − gt sin α + x : 3t
ð14Þ
2A is the maximum run-up height (see Fig. 1). The bore height is A therefore h0 = cosα . The shoreline trajectory (xs) and velocity (Us) are: xs ðt Þ = 2t
pffiffiffiffiffiffi 1 2 gA − gt sin α 2
pffiffiffiffiffiffi Us ðt Þ = 2 gA − gt sin α
ð15Þ
ð16Þ
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Eq. p ffiffiffi (16) allows computation of the swash event duration Tswash = A 4 pffiffigsinα . The collapse of the bore is always located at x = 0. The numerical tests discussed in this section are carried out using m = tan α = 0.1. The initial bore height is h0 = 1m, and the Courant number Cr = 0.8, in all the tests. A convenient way to express the computational resolution in dimensionless form is to use the number of computational cells used to resolve the theoretical maximum run-up distance, i.e. the ratio Nswash = x s,max / Δx. Here xs,max denotes the maximum run-up predicted by the PW01 solution. We use Nswash = 220, with Δx =0.09 m, which is reasonable for practical computations. In the next section we also confirm that the results are consistent for any Nswash. hmin = 10− 6 m is used in the tests, unless stated otherwise. This value has been chosen since it is the same used in Hu et al. (2001) and within the range of the values suggested in other publications (e.g. Kim et al., 2004). Note that the formula proposed by Toro (2001) yields larger values, in this case hmin = 0.1 tan αΔx = 0.9 × 10− 3 m. It is useful to express the quantities in Eqs. (14)–(17) in dimensionless form. The same scaling arguments in PW01 are used here. Dimensionless quantities pffiffiffiffiffiffiHence x′ =x sin pffiffiffiffiffiffiffiffiffiffi are indicated using a prime. α/A, t V= t sin α g = A, h′ =h cos α/A, U V= U = gA. The maximum run up distance and swash event duration are therefore x′s,max = 2 and T′swash = 4. We note that the PW01 solution, when non-dimensionalised is self similar for any slope and bore height, which implies that all values of A and α are represented by the non-dimensional simulation. First we compare the solution in terms of h′ and U′ resulting from the different tested options. The cell with i =I identifies the numerical shoreline. Fig. 4 shows the prediction of the free surface given by four among the tested options (Option 2a, Option 2ea, Option 3 and Option 4b), chosen as representative of all the techniques used, along with the PM01 solution. The free surface is very well reproduced by the present model (panel (a)) as is the velocity (panel (b)). The differences among the various SBC treatments are concentrated in the tip region. It is evident that the numerical shoreline is located offshore of the analytical one, so
Fig. 4. Numerical results for the PW01 problem at different time steps. (a) h, (b) U. Option 2a (green lines), Option 2ea (red lines), Option 3 (magenta lines) and Option 4 (blue lines) compared with PW01 (black solid lines, run-up phase; black dashed lines, backwash).
499
that the run-up distance is always underestimated. Because of this underestimation of x′s, U′s is always underestimated (overestimated) in the run-up (run-down) phase. In the proximity of the numerical tip the numerical solution does not behave linearly. Fig. 5 is a close-up of the tip region at different times; note the different scales used in the panels. The figure shows four snapshots of h′ (upper panels) and the associated U′ (lower panels). As noted before, most of the differences are seen in the velocity profile. Note that Option 3 shows a significant departure from the analytical curve at the tip during the run-up (t′ = 0.625, t′ = 1.87). Later it will be shown that U′s also shows a noisier behaviour during the backwash. Note that due to the very small steepness of the free surface near the tip, the mismatch between the theoretical and the numerical shoreline is quite significant along the x′ axis. This feature is also present at different resolutions; i.e. it is a ‘self similar’ shape. The different options are now compared in terms of x′s and U′s. First we analyze the performances of Option 1 and Option 2a (see Fig. 6). x′s is normalized with the maximum run-up distance x′s,max. In panel (c) of the figure the phase diagram is also shown, this relates the two variables directly. Results in terms of non-dimensional predicted maximum runup (x′s,max,num) and non-dimensional swash period T′swash are shown in Table 1 for all the options tested. The two options displayed give similar x′s and U′s. U′s is underestimated immediately after the bore collapse, but later this mismatch reduces significantly. It is also interesting to note that the velocity curve is not smooth and U′s grows with a sawtooth-like shape. These behaviours are common to most of the methods tested here. Near the end of the swash cycle the numerical solution is affected by some noise. This problem is emphasized in the phase diagram. Here ‘spikes’ occur in the backwash phase, characterized by a sudden change in U′s while x′s remains constant. During the drying process the nominal shoreline, sometimes, does not retreat for several time steps; the subsequent retreat of the shoreline (often by more than one cell in a time step), causes Us to appear irregular. All the methods tested also underestimate the swash period because of the underprediction of x′s. Fig. 7 shows a comparison of two of the options belonging to the Option 2 family. Only the options that performed best are shown. Option 2c produces a significant underestimation of x′s, and is therefore no longer considered. Also we notice that Option 2d did not produce any benefit compared to Option 2ea and Option 2b, therefore this is not shown as well. Options 2b and 2d, and Option 2eb are not significantly different from Options 2a and 2ea respectively. Furthermore, Options 2b and Options 2eb are found to be impractical. Results generally show that the shoreline trajectory is insensitive to the wave speed estimates used in the SBC. This is evident when Option 2a and Option 2ea are compared directly. If UI + 1 = UI is used (this is done in Options 2b and 2eb) no improvement is seen in the run-up estimation. However we use this approach only in the run-up phase, since it produces very poor results if applied during the backwash as well. This limits the usefulness of this method. Overall all the variants of Option 2 give very similar results. A better performance in the backwash and the capability of preserving the symmetry of the x′s curve make Option 2ea preferable. Therefore, in the following we focus only on this approach within the Option 2 family. Option 3 results differ from those from other methods, mainly in terms of velocity. Fig. 8 shows a comparison between the options that use the basic idea of Brocchini et al. (2001) with modifications; Option 2ea is shown for reference. Note that the shoreline parabola seems to lose symmetry. The maximum of the run-up slightly anticipates the analytical one and the backwash phase lasts longer than the run-up one. To understand this we note that the cells are flooded or dried according to the sign of UI (see Fig. 3). In particular, if UI b 0, the cells with hbhmin are removed from the computational domain. In the other methods the criterion used to dry the ‘shallow’ cells is independent of the value of UI. This sign reversal at x′s =x′s,max induces faster initial shoreline retreat than for the other methods. Later the x′s retreat is slower.
500
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Fig. 5. Numerical results for the PW01 problem. Details at the tip of the swash lens at different time steps. Upper panels: Water depth (h), lower panels: velocity (U). Option 2a (green lines), Option 2ea (red lines), Option 3 (magenta lines) and Option 4 (blue lines) compared with PW01 (black solid lines, run-up phase, black dashed lines, backwash).
The most evident feature of Option 3 is the particularly ‘noisy’ velocity curve in the run-up phase. Again the wetting/drying technique plays an important role. Hence, in Options 4a and 4b the ‘depth history’ mechanism is removed and replaced with the same approach as in Option 2ea. Option 4a (not shown) gives very similar results to the Option 2 ‘family’ for x′s and U′s. However, Table 1 shows poorer results for x′s,max. Therefore it is no longer considered.
As seen in Fig. 8, Option 4b, on the other hand, improves significantly the prediction of the velocity after bore collapse; thereafter differences with the other methods reduce. The x′s curve is also better behaved, (i.e. symmetric and smooth). Note that for Option 3 and 4 U′s decreases with time and then increases to reconnect to the same curve followed by the other methods at around t′ = 0.625. A different approach for the flux computation is that implemented in Option 5a and b. As explained above the two alternatives differ only for
Fig. 6. Comparison between PW01 solution and numerical results. Panel (a): x′s(t)/x′s,max. Panel (b): U′s(t). Panel (c): phase diagram (x′s vs U′s). PW01 solution (black lines) Option 1 (blue lines) and Option 2a (red lines).
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505 Table 1 Maximum run up distance x′s.max,num and swash period T′swash for the tested SBC techniques. Model
x′s.max,num
T′swash
hmin [m]
PW01 Option Option Option Option Option Option Option Option Option Option Option Option Option
2.00 1.767 1.749 1.749 1.749 1.749 1.731 1.767 1.740 1.722 1.749 1.722 1.731 1.731
4.00 3.936 3.873 3.878 3.875 3.875 3.804 3.875 3.900 3.856 3.865 3.91 3.852 3.807
10− 6 10− 6 10− 6 10− 6 10− 6 10− 5 10− 6 10− 6 10− 6 10− 6 10− 6 10− 6 10− 5
1 2a 2b 2d 2ea 2ea 2eb 3 4a 4b 5a 5b 6
Nswash = 220 in all tests.
expressions for the flux in the rarefaction fan. Option 5b gives slightly better results, and so is shown here. However, Fig. 8 shows that Option 5b still performs similarly to Option 2ea. Note that Option 4b and Option 5b differ only by the choice of Ffan. The use of the heuristic approach in Brocchini et al. (2001) gives a higher U′s in the early run-up. Note that in the first time steps the velocity is substantially overestimated. As shown in Fig. 8 this results also in a slightly better prediction of the shoreline trajectory immediately after bore collapse. In applying Option 6 we use hmin =10− 5 m, it being the minimum value that gave stable results in the backwash. Fig. 9 shows the phase diagram for the results of Option 6 compared to Option 2ea using hmin = 10− 5 m (blue line) and hmin = 10− 6 m (green line). There is no significant difference between Option 6 and 2ea but note the difference in U′s in the backwash, due to the larger hmin value for Option 6. The fact that the velocity in the run-up phase is less sensitive to hmin is a feature of all the tested techniques and is a relevant aspect for practical applications. The numerical SBC approximates Eq. (8), which therefore is not satisfied exactly. It is instructive to assess the degree to which Eq. (8) is satisfied in each SBC. In order to do so we introduce the error: SBC =
R xs ðt Þ − 0t Us ðτÞdτ xs;max
ð17Þ
which shows the accumulated error in shoreline modelling. This is plotted in Fig. 10. All results show the underestimation of xs seen in earlier plots (see Fig. 6–9). The monotonic increasing curves show that this grows throughout the swash cycle. Note, however, that for Option 2ea, 4b and 5b a significant portion of the mid to late run-up phase
501
shows no increase in error (i.e. good modelling), after an initially rapidly increasing error. The backwash phase marks the start of further underprediction: see Fig. 10. Of further note are the significantly poorer performance of Option 3 and the corresponding better performance of Option 4b, for which the improved prediction of U′s immediately after bore collapse is apparent. All the alternatives of the Option 2 ‘family’ show a very similar magnitude of εSBC as do Option 5b and 6 (these latter two are not shown here). Thanks to the lower error in the velocity prediction Option 4b performs better in terms of SBC balance. However, the results in the backwash do not differ significantly from those of the other methods. 4.1. On grid spacing and hmin Both the domain resolution and the value of hmin are critical aspects in modelling of the swash zone. In Fig. 11 x′s and U′s are plotted for different mesh sizes (i.e. for different Nswash =xs,max / Δx) for Option 2ea (Panels (a) and (b)), chosen as representative of the ones derived from Hu et al. (2001), and for Option 4b (Panels (c) and (d)); which improves the velocity description at bore collapse. Note that the mismatch of the velocity at the early run-up stages reduces with the increasing resolution (Panels (b) and (d)), Option 4b improves the prediction of U′s at any resolution in the run-up phase. The initial U′s is the same for any resolution. This is because the solution of the Riemann problem at t = 0 is the same at any resolution. Since the highest gradients occur at the beginning of the event, the resolution of the steep solution is crucial. Hence, the decreasing cell dimension reduces the time interval during which the initial conditions influence the accuracy. The influence of hmin and Nswash is here studied by analyzing the errors in the maximum run-up prediction and in U′s as a function of the ratio h′min /Δx′=hmin /Δx tan α) and of Nswash. For the run-up we define: run−up =
xs;max;num − xs;max × 100 xs;max
ð18Þ
To evaluate the accuracy of U′s for the whole of the swash cycle we use the root mean square error (RMSE) defined as:
U sV
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uN ts uX ðU Vns −U Vs ðtn ÞÞ2 =t : Nts n=1
ð19Þ
Nts is the number of time steps in the swash cycle, U′ns is the numerical prediction of U′s at the n-th time step and U′s(tn) is the corresponding analytical value. In Fig. 12 contours of εrun-up (upper panel) and εU′s (lower panel) are shown as functions of h′min / Δx′ and Nswash for Option
Fig. 7. Comparison between PW01 solution and numerical results. Panel (a): x′s(t) / x′s,max. Panel (b): U′s(t). Panel (c): phase diagram (x′s(t) / x′s,max vs U′s). PW01 solution (black lines) Option 2a (blue lines) and Option 2ea (red lines).
502
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
Fig. 8. Comparison between PW01 solution and numerical results. Panels (a) and (d): x′s(t)/x′s,max. Panels (b) and (e): U′s(t). Panels (c) and (f): phase diagram (x′s(t)/x′s,max vs U′s). Panels (a), (b), (c): PW01 solution (black lines) Option 2ea (blue lines) and Option 3 (red lines). Panels (d), (e), (f): PW01 solution (black lines) Option 4b (blue lines) and Option 5b (red lines).
2ea. In the plots h′min / Δx′ varies from 1 × 10− 1 to 1 × 10− 5, whilst the maximum Nswash is 1550. Tests performed with smaller h′min / Δx′ and larger Nswash resulted in poor performances in the backwash, so they are not shown. The first conclusion that can be drawn from the figure is that εrun-up monotonically decreases as either h′min / Δx′ decreases or Nswash increases. For εU′s the reduction of h′min / Δx′ is not always beneficial. For h′min / Δx′ b 1 × 10− 3, εU′s does not necessarily indicate convergence for
Fig. 9. Phase diagram (x′s(t)/x′s,max vs U′s(t)) for Option 6 using hmin = 10− 5 m, (red line) compared with Option 2ea using hmin = 10− 6 m (green line) and hmin = 10− 5 m (blue line). PW01 solution is also plotted (black line).
increasing Nswash. This is due to the increasing error in the backwash phase for very small h′min / Δx′. The contour plot of εrun-up shows that, for h′min / Δx′ N 1 × 10− 3, decreasing h′min / Δx′ for constant Nswash produces a greater benefit than increasing resolution alone. For h′min / Δx′ b 1 × 10− 3 resolution becomes increasingly more important in reducing εrun-up. Obviously, while decreasing h′min is not computationally expensive increasing the resolution is. Therefore, ideally we choose a suitable value of h′min / Δx′ that preserves accuracy whilst saving on resolution.
Fig. 10. SBC balance for Options 2ea, 3, 4b and 5b.
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
503
Fig. 11. x′s(t) and U′s(t) for different values of Nswash. Panels (a) and (b): Option 2ea. Nswash = 125 (magenta line), Nswash = 250 (green line), Nswash = 500 (blue line), Nswash = 1000 (red line). Panels (c) and (d): Option 4b. Nswash = 125 (magenta line), Nswash = 250 (green line), Nswash = 500 (blue line), Nswash = 1000 (red line). PW01 solution is always plotted with a black line.
We conclude that h′min / Δx′ = 1 × 10− 3 is an optimal value since it guarantees that both εrun-up and εU′s indicate convergence for increasing Nswash; furthermore smaller values of the ratio produce
moderate benefits for the run-up but a poorer performance for the velocity of the backwash, causing the increase in εU′s. Option 4b produces very similar results to the ones shown in Fig. 12, with εrun-up
Fig. 12. Upper panel: contour plot of εrun-up for Option 2ea for different values of h′min / Δx′ and Nswash. Lower panel: contour plot of εU′s for Option 2ea for different values of h′min / Δx′ and Nswash. Note that max(U′s) = 2.
504
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
slightly larger at all resolutions. However, we found that εU′s is larger when h′min / Δx′ N 1 × 10− 3 as a result of worse performance in the backwash. Note that h′min /Δx′ =hmin/(tan αΔx) = 1 × 10− 3 is smaller than the ratio suggested by Watson et al. (1992) (i.e. h′min /(tan αΔx) = 1 × 10− 1). Furthermore, note that transformation from an x-axis running along the bed to a horizontal axis, leaves the ratio hmin /Δx unchanged. 5. Conclusions A comparative study of different treatments of the shoreline boundary conditions in Godunov-type schemes has been presented in this paper. We focus on swash motion and use the PW01 problem as a demanding benchmark for a state of the art WAF solver. Most of the tested SBCs show comparable performances in terms of x′s and U′s. The modelling of the bulk of the flow, on the other hand, is not significantly affected by the choice of the SBC. A critical aspect when modelling the shoreline motion with the solvers used here is the accuracy of the velocity immediately after bore collapse. Most of the techniques tested exhibit a very similar behaviour, with the velocity at the tip increasing in time and reconnecting with the theoretical curve at a later stage. The duration of this ‘growing curve’ depends strongly on the numerical resolution. A practical solution to this problem is the use of a modified version of the Brocchini et al. (2001) technique, referred to here as Option 4b. This improves the prediction of the velocity at early stage and gives better results than the original Brocchini et al. (2001) formulation. Among the approaches derived from the Hu et al. (2001) ‘minimum wet depth’ technique, Option 2ea performs slightly better and proves to be robust at any resolution both in the run-up and the backwash phase. The sensitivity analysis carried out shows that the improvement in U′s in Option 4b (compared to Option 2ea) is more significant at resolutions Nswash b 500. At this resolution, the initial velocity is still better modeled by Option 4b, but the maximum run-up is slightly better predicted by Option 2ea. Hence we suggest to use either one or the other according to the grid size required in the problem under consideration. In all the methods tested the prediction of U′s in the backwash phase seems to be a demanding aspect. It is shown that the drying process is usually not smooth, especially at lower resolutions, affecting, in turn, the U′s prediction. The fact that in the PW01 solution the free surface is always tangential to the bottom at the tip, implies that in the backwash a thin film of water is present, causing the numerical round-off errors to be important in the solution. This effect is less relevant when friction is included. We also analysed the relative role of hmin and Nswash in the accuracy of the maximum run-up and U′s prediction. The outcome of this analysis is that h′min / Δx′ ≈ 1 × 10− 3 can be used as a predimensioning rule for testing a code against the PW01 solution. When Nswash = 1500 and the suggested ratio is used, εrun-up = 5% a value that is acceptable for engineering purposes. Slightly better performances can be obtained by reducing h′min / Δx′, but with possibly poorer performances for U′s. Although this resolution may seem quite high, it is important to remember that the (a) methods presented are relatively computationally inexpensive, and parallel computing can reduce the computation time for two dimensional simulations; and (b) for practical engineering problems this kind of resolution is unlikely to be necessary because of the presence of friction, and pure dam-break problem such as that of the PW01 solution are rarely encountered in practical applications. Other numerical models, such as those based on the method of characteristics, achieve higher accuracy than those used here. However, those methods are usually significantly more demanding to implement, especially in two dimensions, and are computationally more expensive, and so they rarely find practical engineering application. Much of the mismatch discussed here stems from the fact that the analytical solution is tangential to the bottom at the shoreline, a con-
dition that is not met by the numerical prediction. This makes the PW01 solution a very demanding test for swash zone applications. When friction is present in the system this condition no longer pertains, and in turn results are expected to be more accurate. Similarly, for other swash events (see e.g. Hibberd and Peregrine, 1979) bore height will diminish when encountering the shore, thus alleviating the problem with estimating the initial U′s that we illustrated here. Nevertheless, good modelling of the PW01 solution will certainly produce good performances for swash flows with friction or in less demanding cases. The implications of the above results for both swash zone morphodynamic modelling and wave overtopping are currently being studied and represent the natural continuation of this work. Acknowledgements This study was supported by the EPSRC through the research project “Experimental and Numerical Modelling Study of Swash Zone Hydrodynamics and Sediment Transport” (EP/EO10407/1). Their support is gratefully acknowledged. David Kelly of the University of Nottingham, Civil Engineering Department is also acknowledged for the many productive discussions during this research. Appendix A. The WAF technique The solution of the homogeneous problem is carried out using a WAF approach for the computation of the inter-cell flux appearing in Eq. (4). This requires the solution of a subgrid Riemann problem. At this purpose we chose to use a HLL (Harten et al., 1983) approximate Riemann solver. This is widely used and performs well in cases of critical rarefactions and dry fronts. The solver assumes the same wave structure for any inter-cell problem solution, i.e. two waves that separate three regions of constant states. The propagation speed estimates for the two waves (SL, SR) are (Toro, 1999): SL = UL − CL qL
SR = UR + CR qR
ðA:1Þ
where the subscripts L and R indicate the left and right state atffi each pffiffiffiffiffi inter-cell boundary. C is the wave celerity in shallow water ( gh) and qK (K = L,R) is given by:
qK =
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8v 3 u 2 ⁎ ⁎ > u > u < t1 4 h + hK h 5 h2K
2
> > :
1
if if
h⁎ N hK ð22Þ
ðA:2Þ
⁎
h V hK
where h⁎ is an estimate for the solution in the so called ’star region’, i.e. the intermediate region between the two waves. This estimate is provided by a Two-Rarefaction approximate Riemann solver: ⁎
h =
2 1 1 1 ðCL + CR Þ + ðUL − UR Þ : g 2 4
ðA:3Þ
Following Toro (1999) is possible to obtain the HLL estimate of the flux in the star region: = FHLL ⁎
SR FL − SL FR + SR SL ðWR − WL Þ : SR − SL
ðA:4Þ
Once the flux in the star region is known, it is possible to compute the inter-cell flux according to the WAF technique. In the TVD (total variation diminishing) version of the WAF, its expression is given by the relationship: n + 12 + 12
Fi
=
2 n 1 n 1X ðlÞ F Wi + F W i + 1 + signðcl Þl ΔFi + 1 2 2 2 l=1
ðA:5Þ
R. Briganti, N. Dodd / Coastal Engineering 56 (2009) 495–505
where F(Wni ) = FnL is the flux vector computed considering the left state and F(Wni + 1) = FnR is the flux related to the right state. Note that they are both computed at the time level n. ϕl (l = 1,2) is a WAF limiter function. In this study we use the limiter SUPERA (see Toro, 1999). ðlÞ Finally, ΔFi + 1 is the flux jump across each zone identified in the 2 wave structure of the solution. In Eq. (A.5), cl is the Courant number defined for each l-th wave of speed Sl, i.e.: cl =
Sl Δt Δx
ðA:6Þ
References Brocchini, M., Dodd, N., 2008. Nonlinear shallow water equations modeling for coastal engineering. Journal of Waterway, Port, Coastal and Ocean Engineering 134 (2) (6), 104–120. Brocchini, M., Bernetti, R., Mancinelli, A., Albertini, G., 2001. An efficient solver for nearshore flows based on the WAF method. Coastal Engineering 43 (2), 105–129. Carrier, G., Greenspan, H., 1958. Water waves of finite amplitude on a sloping beach. J. Fluid Mech. 4, 97–109. Dodd, N., Stoker, A.M., Calvete, D., Sriariyawat, A., 2008. On the formation of beach cusps. Journal of Fluid Mechanics 597, 145–169. Harten, A., Lax, P., van Leer, B., 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25 (1), 35–61. Hibberd, S., Peregrine, D.H., 1979. Surf and run-up on a beach. Journal of Fluid Mechanics 95, 323–345. Hubbard, M.E., Dodd, N., 2002. A 2D numerical model of wave run-up and overtopping. Coastal Engineering 47 (1), 1–26. Hu, K., Mingham, C.G., Causon, D.M., 2001. Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations. Coastal Engineering 41 (4), 433–465.
505
Kennedy, A.B., Chen, Q., Kirby, J.T., Dalrymple, R.A., 2000. Boussinesq modeling of wave transformation, breaking and runup. I: 1D. Journal of Waterways Port Coastal and Ocean Engineering-ASCE 126, 39–47. Kim, A., Cho, Y.S., Kim, W.G., 2004. Weighted averaged flux-type scheme for shallowwater equations with fractional step methods. Journal of Engineering Mechanics 130 (2), 152–160. Lynett, P., Wu, T.-R., Liu, P.L.-F., 2002. Modeling wave runup with depth-integrated equations. Coastal Engineering 46 (2), 89–107. Pan, C.-H., Lin, B.-Y., Mao, X.-Z., 2007. Case study: numerical modeling of the tidal bore on Qiantang River, China. Journal of Hydraulic Engineering 133 (2), 130–138. Peregrine, D.H., Williams, S.M., 2001. Swash overtopping a truncated plane beach. Journal of Fluid Mechanics 440, 391–399. Press, W., Flannery, B., Teukolsky, S., Vetterling, W., 1992. Numerical Recipes: The Art of Scientific Computing. Cambridge Univ. Press, New York. 992 pp. Shen, M.C., Meyer, R.E., 1963. Climb of a bore on a beach. Part 3. Runup. Journal of Fluid Mechanics 16, 113–125. Strang, G., 1968. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 124, 506–517. Thacker, W.C., 1981. Some exact solutions to the nonlinear shallow water wave equations. Journal of Fluid Mechanics 107, 499–508. Titov, V.V., Synolakis, C.E., 1995. Modeling of breaking and nonbreaking long wave evolution and runup using VTCS-2. Journal of Waterway, Port, Coastal and Ocean Engineering 121 (6), 308–316. Toro, E.F., 1999. Riemann solvers and numerical methods for fluid dynamics. Springer, New York. Toro, E.F., 2001. Shock-capturing methods for free-surface shallow flows. Springer, New York. Watson, G., Peregrine, D.H., Toro, E.F., 1992. Numerical solution of the shallow water equations on a beach using the weighted average flux method. In: Hirsch, P., Kordulla (Eds.), Proceedings of the First European Computational Fluid Dynamics Conference. Elsevier, Brussels, pp. 495–502. September 7th to 11th, 1992.