Numerical analogs to Fourier and dispersion analysis: development, verification, and application to the shallow water equations

Advances in Water Resources 26 (2003) 649–662 www.elsevier.com/locate/advwatres

Numerical analogs to Fourier and dispersion analysis: development, verification, and application to the shallow water equations C.M. Szpilka *, R.L. Kolar School of Civil Engineering and Environmental Science, University of Oklahoma, 202 West Boyd Street, Room 334, Norman, OK 73019, USA Received 11 January 2002; received in revised form 3 February 2003; accepted 5 February 2003

Abstract Herein, we present numerical analogs to traditional Fourier and dispersion analyses and validate them with well-characterized phase behavior for classic finite difference and finite element (FE) discretizations of the shallow water equations. Basically, the procedure is to introduce a single wave with known amplitude and phase into the domain, propagate the wave approximately one wavelength using some discretization scheme, and then note its final amplitude and phase. The final state of the wave is then compared with the expected wave form predicted by the continuum equations to determine the propagation behavior of the discretization. After validating the technique, we then examine two case studies: (1) slope limiting schemes within the finite volume framework and (2) lumping coefficients within the selective lumping FE framework. Of the three common slope limiters that we examined, the Superbee limiter has the most promising phase behavior, as it is the least dissipative while maintaining minimal phase error. Using our numerical technique, we were also able to verify the range of values that has been found to be most accurate in practice for the selective lumping coefficient. Ó 2003 Elsevier Science Ltd. All rights reserved. Keywords: Shallow water equations; Finite volume method; Fourier analysis; Dispersion analysis; Numerical phase analysis

1. Introduction The shallow water equations (SWE) are derived from the depth-averaged Navier–Stokes equations under Boussinesq and hydrostatic pressure assumptions; they describe the motion of water bodies wherein the depth is small relative to the scale of the waves propagating on that body. The system of SWE is used to model the hydrodynamics of lakes, estuaries, tidal flats and coastal regions, as well as deep ocean tides. The SWE continuum equations have been successfully discretized using the finite difference (FD) method [13], a selective lumping finite element (SLFE) method [10], and the finite volume (FV) method [1,5]; in addition, modified equations, such as the wave continuity approach [12], have been used. Primitive finite element (FE) discretizations of the SWE were not as successful and often introduced

*

Corresponding author. E-mail addresses: [email protected] (C.M. Szpilka), [email protected] (R.L. Kolar).

numerical noise into the solution. Accurate solutions demand that the algorithm not introduce spurious modes or overdamp the system. In an effort to characterize these features of a successful algorithm, much work has been done to study the propagation characteristics of the FE and FD discretizations in both 1D and 2D frameworks [2,6,7,9,11,13–15,18,19]; but to our knowledge the FV and SLFE methods have not been similarly studied. Analytical Fourier and dispersion analyses can be used to study the propagation characteristics of many discretizations of the SWE. Such analyses provide information about the phasing of the numerical solutions, the damping characteristics and the stability of the discretization method. However, some discretization methods are not suitable for closed form analyses. For example higher-order FV discretizations, which use piecewise linear or higher interpolates, require some type of slope limiting procedure. These algorithms are transient and cannot be written in a closed form for all time. Consequently, we have developed numerical analogs to the traditional analysis techniques.

0309-1708/03/$ - see front matter Ó 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0309-1708(03)00028-9

650

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

With the exception of Kinnmark and GrayÕs 2 Dx-test [11], which only focused on a single wavelength, most of the previous propagation work has focused on the use of analytical analysis methods. In this study we examine numerical analogs for generating Fourier and dispersion plots. Several benefits exist for undertaking a numerical approach: (1) it is possible to study boundary effects, which are often ignored in analytical methods; (2) if the algorithm evolves with time, it is not possible to write the closed form equation required for analytical analysis; and (3) numerical phase results are more closely tied to the algorithm output and thus capture more of the actual simulation behavior relative to the predicted mathematical behavior (including any unforeseen effects of truncation and roundoff error). In the following section we present the system of SWEs and a short description of the FV method applied to this system. Section 3 contains a brief description of the traditional analytical analysis tools and the resulting propagation relationships for several discretizations of the SWE. The procedure for our numerical analogs to the traditional tools is presented in Section 4. Finally, we present some validation results and two case studies in Section 5: (1) slope limiters within the FV framework and (2) lumping coefficients within the SLFE algorithm, and we offer some closing remarks in Section 6.

2. General equations In this study, we examine the linearized one-dimensional system of equations governing shallow water flow. The system of SWE in this form consists of the primitive continuity equation ft þ hux ¼ 0

ð1Þ

and the non-conservative momentum equation ut þ su þ gfx ¼ 0

ð2Þ

where u is the depth-averaged velocity, f is the surface elevation, s is the linear bottom friction factor, g is the acceleration of gravity, h is the bathymetric depth of water (assumed constant), and the subscripts indicate partial derivatives. (Note that the linearized form of the conservative momentum equation is the same as that of the non-conservative form.) For completeness, we also include the linearized form of the generalized wave continuity (GWC) equation ftt þ Gft þ ðG  sÞhux  ghfxx ¼ 0

consider a no-flow velocity condition at the land boundary. All of the solution algorithms that will be compared in subsequent sections use (1) to solve for elevations and (2) to solve for velocities, with the exception of the GWC approach which uses (3) to solve for elevations. The FV method has only recently been applied to shallow water applications. This method invokes local conservation of mass and momentum by integrating over a discrete volume (or cell) and then solves for the boundary fluxes, which may be discontinuous, on each volume using an approximate Riemann solver, such as RoeÕs linearization. Advantages include its ability to capture shocks without introducing spurious oscillations, local and global mass conservation, utilization of the primitive form of the equations, and the ability to handle irregular meshes. The FV method can utilize any order of interpolate within each discrete cell, however piecewise constants (low-order FV) and piecewise linear (high-order FV) interpolates are most common. Fig. 1 shows a schematic of a typical FV discretization with piecewise linear interpolates, wherein each cell contains information about the cell-average value and the functional form of the slope for each variable. The arrows indicate the outward normal flux across the interfaces for cell j; notice that the interpolates need not be continuous across the interface. In the interest of brevity, the reader is referred to [8] for a more complete derivation of the FV discretization applied to the linearized one-dimensional SWE. Herein, we present only the generalized form of the discrete equations (before the application of some limiting procedure). These general equations are valid for both lowand high-order FV approximations, since the assumed form of the interpolates will only determine the representation of the left and right cell values; they can be written as Dt ½hðukj;R þ ukjþ1;L  ukj1;R  ukj;L Þ 2 Dx  aðfkj1;R  fkj;L  fkj;R þ fkjþ1;L Þ

ð4Þ

Dt ½gðfkj;R þ fkjþ1;L  fkj1;R  fkj;L Þ 2 Dx  aðukj1;R  ukj;L  ukj;R þ ukjþ1;L Þ  sukj Dt

ð5Þ

fkþ1 ¼ fkj  j

and ukþ1 ¼ ukj  j

ð3Þ

where G is a numerical coefficient that determines the balance between primitive (large values of G) and pure wave (small values of G) forms. Previous analyses have shown that (3) is preferable to (1) for equal order FE discretizations of the SWE [12]. For the boundary conditions we specify elevation at the open boundary and

j-1

j

j+ 1

Fig. 1. Typical finite volume discretization.

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

where a is the frictionless wave celerity (or the speed at which the wave would travel analytically in the absence 1=2 of friction) and is equal to ðghÞ , and the subscripts j, RðLÞ indicate the rightmost (leftmost) value of u or f within the jth cell. These values are calculated from the known cell-centered average and the assumed functional form of the interpolates within the cell. In this paper, we consider only piecewise constant and piecewise linear approximations where the raw slopes are calculated using values from neighboring cells (explicit slope reconstruction after the computation of the cell average quantities). For stability, some limiting procedure must be applied to the raw slopes so that two adjacent cells do not have slopes of opposite sign. (Refer to Fig. 1 for a discretization where the linear interpolates have been limited. Notice that the jth cell has a slope of zero.) These limiting procedures are also known as total-variation diminishing and essentially non-oscillatory schemes: they prohibit expressing the high-order FV algorithm in a completely closed-form equation since the slopes evolve with the solution. For this study, we have examined three common limiters (as presented in various applications of higher-order methods) and their behavior within the context of the high-order FV discretization of the SWE. These are the minmod limiter [16] expressed as minmod ðL; RÞ ¼ sL Min ðjLj; jRjÞ if LR > 0; else minmod ¼ 0

ð6Þ

the Superbee limiter [4] expressed as Superbee ðL; RÞ ¼ sL Max ½0; Minð2jLj; sL jRjÞ; MinðjLj; 2sL jRjÞ

ð7Þ

and the vanLeer limiter [3] expressed as vanLeerðL; R; CÞ ¼ 0:5sC Minð0:5C; minmodðL; RÞÞ

ð8Þ

where sx is the sign of x and the raw slopes are calculated using standard FDs: L ¼ ðbj  bj1 Þ=Dx R ¼ ðbjþ1  bj Þ=Dx C ¼ ðbjþ1  bj1 Þ=ð2 DxÞ

ð9Þ

where b can be either u or f, and the subscripts indicate the cell average values for the jth cell. Although mathematically each of these limiters takes on a slightly different form, they are essentially performing the same tasks: (1) to ensure a monotonic increase or decrease in the function values by prohibiting slopes of opposite sign in adjacent cells (i.e., LR > 0 tests that slopes have the same sign) and (2) to limit the magnitude of the boundary cellÕs slope to prevent overshoot and undershoot of the boundary conditions. It is instructive to compare the discrete FV equations with their continuum counterparts, Eqs. (1) and (2).

651

Returning to (4) and (5), and noting that for piecewise constants uj;L ¼ uj;R and fj;L ¼ fj;R , we can write the loworder discrete FV representations as Dt ½hðukjþ1  ukj1 Þ 2 Dx  aðfkj1  2fkj þ fkjþ1 Þ

ð10Þ

Dt ½gðfkjþ1  fkj1 Þ 2 Dx  aðukj1  2ukj þ ukjþ1 Þ  sukj Dt

ð11Þ

¼ fkj  fkþ1 j

and ukþ1 ¼ ukj  j

We note that the low-order FV is basically a central difference approximation (in space) that introduces a centered approximation of a second derivative into the discrete equations, viz. the terms of (10) and (11) that are preceded by the parameter a. Kinnmark [12] and Atkinson [2] have noted that a characteristic feature of the GWC and quasi-bubble algorithms, which produce non-oscillatory solutions, is the presence of a secondorder space derivative. This derivative is correlated with the ability of an algorithm to successfully propagate or damp out undesirable short wavelength oscillations. Chippada et al. [5] note that a centered spatial approximation to the first derivative, which is derived by using an arithmetic average for the flux approximation, results in an unstable algorithm, i.e., neglecting the terms preceded by a in Eqs. (10) and (11). Thus the source of the centered spatial approximations determines whether they improve or hinder the solution. Unlike the GWC and quasi-bubble algorithms, which manipulate the governing equations in continuum (GWC) or discrete (quasi-bubble) forms, a Taylor series analysis on Eqs. (10) and (11) reveals that the Godunov type flux approximation, using RoeÕs linearization, introduces the second-order spatial derivative in the leading term of the truncation error. In other words, the discrete equations presented in (10) and (11) are still consistent with the original continuum equations, since the second-order spatial derivative is introduced to the leading term of the error instead of in the equations themselves. This additional term has the effect of reducing the overall accuracy of the approximations to first-order and it adds diffusion to both the continuity and momentum equations. It also stabilizes the algorithm despite the fact that the first-order spatial derivatives come from the arithmetic average of the left and right states. The overly dissipative nature of this algorithm will be shown in Section 5.1. For a more detailed derivation of the FV applied to the non-linear SWE (see [1,5]). For discrete equations from the other discretization techniques (FE, SLFE, GWC, FD), please refer to the appropriate references [9,10,12,13].

652

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

3. Analytical propagation analysis Two frequently used tools for analytical propagation analysis exist: Fourier analysis, which uses the fully discretized equations; and dispersion analysis, which discretizes only the spatial variables, leaving time continuous via the harmonic form of the equations. As will be seen in Section 4.2, these two can be related through the phase speed. The desired propagation characteristics for any algorithm can be found by looking at the continuum equations, which, in this case, exhibit a monotonic dispersion relationship. It is also desirable for the algorithm to damp only those waves that propagate with significant phase error (typically high-frequency components on the order of 2–5 Dx).

where the subscript n has been left off to indicate that only one component of the series solution is being examined at a time due to the linearity of the governing equations. Here the analytical wave speed, or celerity, is expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð14Þ c ¼ gh  ðs=2rÞ From Eq. (13), we can see that in a single time step the wave will decrease in amplitude by a factor of expðs Dt=2Þ and travel a distance equal to c Dt. Given the analytical solution, we now turn to the discrete system and present the Fourier analysis for a general discretization of the SWE. For the fully discretized system, a single component of the Fourier series, represented as

3.1. Fourier analysis

bkj ¼ b0 eixk Dt eirj Dx

Fourier analysis is a useful tool to analyze the propagation behavior of both analytical and discrete solutions of homogeneous, linear difference or differential equations. It utilizes the assumption that separation of variables can be used to express the solution as periodic in time and space, and thus any solution that is periodic on an infinite interval or piecewise continuous on a finite interval can be represented by a Fourier series expansion. A solution (discrete or analytical) that satisfies these minimal constraints can be represented as a complex Fourier series

where j, k are the spatial and temporal discretization indices, respectively, and b0 is the amplitude, can be substituted into the discrete representation of Eqs. (1)– (3) for each independent variable. The discrete per-timestep propagation factor is defined as:

bðx; tÞ ¼

1 X

An eixn t eirn x

ð12Þ

k0

bkþ1 j ¼ eix Dt bkj

ð15Þ

ð16Þ

and definitions (15) and (16) are substituted into the set of equations for the algorithm under study. From this substitution a new system results,      M11 M12 f0 0 ¼ ð17Þ M21 M22 u0 0

n¼1

where, for the SWE, b ¼ u or f, i ¼ ð1Þ1=2 , xn is the temporal frequency of the solution, rn ¼ 2p=Ln is the wavenumber and Ln is the nth wavelength, and An is the Fourier coefficient for component n. The key to Fourier analysis is that the linearity of the difference or differential equation allows one to examine a single Fourier component of the solution at a time. Such an analysis, also known as von Neumann analysis, results in two related plots: phase error versus wavelength and damping ratio versus wavelength. Similarly, one can introduce a complex propagation factor, T , which is the ratio of the discrete wave to the analytical wave after the time it takes for the analytical wave to propagate one wavelength, as in [9,13]. The latter is the approach that will be taken in this study. The starting point for such an analysis is to determine the analytical wave properties from the continuum equations. Rather than rederive this relationship, we refer the interested reader to [9] and present their results herein. The per-time-step propagation factor, k, for the analytical wave is given as k ¼ esDt=2 e irc Dt

ð13Þ

where the form of the matrix coefficients will vary depending upon the spatial discretization. A non-trivial solution will exist only when the determinant of the coefficient matrix is zero, which results in a polynomial in the discrete per-time-step propagation factor k0 , whose complex roots can then be written as a function of r. The stability of the numerical algorithm can be assessed by computing jk0 j and verifying that its magnitude is less than or equal to 1 for all desired values of r. Meanwhile, damping and phase properties of the numerical solution, relative to the analytical solution, are determined by computing the complex (per-wavelength) propagation factor T , which is defined as !Nn k0 T ¼ ð18Þ k where Nn is the number of time steps required for the analytical wave to propagate exactly one wavelength (Nn ¼ Ln =c Dt). The magnitude of T is the ratio of the amplitude of the discrete wave to the amplitude of the continuum wave after propagating Nn time steps. If jT j is greater

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

than 1, the discrete wave has greater amplitude than the analytical wave, but it is not necessarily unstable, since T is not simply a measure of the discrete wave itself but is a comparison to the continuum wave. While if jT j is less than 1, the discrete wave has smaller amplitude than the analytical wave and the solution is damped. The argument of T is the phase of the discrete wave compared to the analytical wave, where a positive phase indicates lead and a negative phase indicates lag. For a numerical algorithm to perfectly match the analytical solution of the model problem, the magnitude of T should equal 1.0 and the phase of T should equal 0.0 for all components of the solution.

3.2. Dispersion analysis Again, this analysis tool utilizes the assumption that separation of variables can be used to express the solution of the differential equations as periodic in time and space. However, it differs from Fourier analysis in that the time variable remains continuous. Additionally, the analysis results in a relationship between the temporal frequency and the wavenumber, often referred to as a dispersion relation, as opposed to the propagation factor in Fourier analysis. In practice, this dispersion relation can be thought of as an intermediate step in the analytical solution obtained in Fourier analysis. It is this relation that allows one to write the per-time-step propagation factor, Eq. (13), without explicit dependency on x. Platzman [15] was the first to apply dispersion analysis to the SWE; he demonstrated that a monotonic curve indicates a numerical solution free from spurious 2 Dx oscillations. A folded curve indicates aliasing of wave components, with one wave corresponding to the long physical wave while the other corresponds to shortwavelength noise in the solution. It has since been demonstrated that for linear elements, but not for quadratic or other elements, a folded curve is a necessary and sufficient condition for the appearance of spurious modes in the numerical solution of the SWE [12]. However, a non-folded curve does not necessarily imply that the solution will be free from these spurious modes. Atkinson [2] has shown that a non-folded dispersion curve in 1D is not always indicative of 2D simulation behavior in that a scheme which exhibits a non-folded curve in 1D analysis can exhibit a folded curve in 2D analysis. However, his analysis did verify that a folded dispersion curve in 1D is a necessary and sufficient condition for a folded dispersion curve in 2D. Therefore, as a first step, it is instructive to examine the 1D dispersion relationships for new algorithms before proceeding with the analysis in higher dimensions (i.e., a folded curve in 1D indicates that further study is not warranted).

653

The starting point of dispersion analysis is to write the harmonic form of each equation by substituting the temporal harmonic components expressed as f ¼ feixt

and

u ¼ ueixt

ð19Þ

into the differential equations. Here f is the spatial harmonic of the elevation and u is the spatial harmonic of the velocity. The analytical dispersion relation can be derived by substituting the related continuum spatial harmonics (Fourier components) f ¼ f0 eirx

and

u ¼ u0 eirx

ð20Þ

into the resulting ordinary differential equations. The solution of the resulting system of equations for the magnitudes u0 and f0 determines the relationship between x and r. The magnitude of x is then plotted against r to develop the dispersion curve. (These plots can be represented either in physical or dimensionless variables; we adopt the former herein.) For the continuum equations, the above analysis steps result in the following system of equations: ixf0 þ irhu0 ¼ 0 and

ixu0 þ su0 þ irgf0 ¼ 0

which has the dispersion equation pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii 1h x ¼ is 4r2 gh  s2 2

ð21Þ

ð22Þ

For the discrete equations, the dispersion characteristics can be determined by substituting the discrete harmonic spatial solutions, expressed as f ¼ f0 eijr Dx

and

u ¼ u0 eijr Dx

ð23Þ

into the spatially discretized ordinary differential equations. Once Eq. (23) are substituted into the discrete equations, the analysis proceeds as for the continuum equations. The result is a polynomial expression for the discrete temporal frequency, x0 , as a function of r, which is then solved for its complex roots. The magnitude of the roots versus wavenumber is plotted as a dispersion curve. For more background information on dispersion analysis as applied to the SWE see [6,7,15,18]. 3.3. Analytical propagation relationships for the SWE In this section we examine the staggered finite difference (SFD) scheme of Leendertse, a leap-frog primitive finite element (PFE) method, the SLFE method of Kawahara, the GWC approach of Kinnmark, and the low-order FV method. (The high-order FV is the subject of Section 5.2.) A comparison matrix of the discretization features for these algorithms is provided below in Table 1; see the listed references for more details. Unless otherwise noted, the discretization features apply to both the continuity and momentum equations.

654

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

Table 1 Discretization features in time and space for study algorithms SFD [13] PFE [9] SLFE [10] GWC [12] Low-order FV [5]

Time-stepping algorithm

Spatial interpolates

Continuity of spatial interpolates

Crank–Nicolson Leap-frog Two-step explicit Three-level centered for Eq. (3), two-level centered for Eq. (2) One-step explicit

N/A Piecewise Piecewise Piecewise Piecewise

Continuous Continuous Continuous Continuous Discontinuous

linear linear linear constant

Table 2 Propagation factors from analytical Fourier analysis k ¼ expðs Dt=2Þ expð irc DtÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Shs DtT ðs DtT Þ2 þ4ST k0 ¼ 1 þ 2ðrT Sh2 Þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Continuum SFD

k0 ¼

PFE

1  s Dt  8K 2i 4K  ð4K þ s Dt=2Þ2

SLFE

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k0 ¼ 12 2E  s DtJ ðs DtJ Þ2 þ l2 ghF 2

GWC

P3 k0 þ P2 k0 þ P1 k0 þ P0 ¼ 0

Low-order FV

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k0 ¼ 12 2 þ laðA  2Þ  s Dt ðs DtÞ2 þ ghðlBÞ2

3

2

A ¼ 2 cosðr DxÞ, B ¼ 2i sinðr DxÞ, C ¼ 2 cosð2r DxÞ, D ¼ 2i sinð2r DxÞ

2 E ¼ 2þe þ 1e A þ l 8gh ðC  2Þ, F ¼ 2þe B þ 1e D  s 2Dt 2þf B þ 1f D 3 6 3 6 3 6



2þf 1e 2þe s Dt 2þf J ¼ 1f ðC þ 2Þ 1e  s 2Dt 1f  s Ds 2þf þ A 2þe þ 3 A 6 þ 3  2 3 6 6 3 6 3 K ¼ ½3al sinðr DxÞ=½4 þ 2 cosðr DxÞ2 , S ¼ 2i sinðr Dx=2Þgl2 hðeir Dx  1Þ, T ¼ eðir DxÞ=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ gh  s2 =ð4r2 Þ, r ¼ 1 þ s Dth, l ¼ Dt=Dx 2

2

2

K1 ¼ 6DxDt, K2 ¼ 3DxDt2 , K3 ¼ 3DxDt3







P0 ¼ K3 2 þ A2 þ ðK2 ðG þ sÞ  K1 GsÞ 1 þ A4 þ a3 gh 1  A2 s  Dt2







P1 ¼ 3K3 2 þ A2 þ ðK2 ðG þ sÞ þ K1 GsÞ 1 þ A4 þ gh 1  A2 Dt2 ða3  a2 Þ þ sða3 þ a2 Þ  gh8 B2 ðG  sÞ







P2 ¼ 3K3 2 þ A2 þ ðK2 ðG þ sÞ þ K1 GsÞ 1 þ A4 þ gh 1  A2 Dt2 ða2  a1 Þ þ sða2 þ a1 Þ  gh8 B2 ðG  sÞ







P3 ¼ K3 2 þ A2 þ ðK2 ðG þ sÞ  K1 GsÞ 1 þ A4 þ a1 gh 1  A2 s þ Dt2

The per-time-step propagation factor k0 and the dispersion relation that result from the above analytical treatment of the discretization schemes are summarized in Tables 2 and 3. The analytic propagation factor is as in Eq. (13) and the analytic dispersion relation is as in Eq. (22). Recall that the discrete phase characteristics are distinguished from the continuum characteristics with the prime notation, e.g., k0 versus k. In Table 2, note that the GWC formulation results in three roots due to the second order time derivative. In the one-dimensional context, two of these are physical roots (backward and forward traveling waves), while the third is a numerical artifact of the time differentiation. In the interest of brevity, we report only the cubic polynomial in k0 instead of the actual roots. For the SFD algorithm, h is the time weighting parameter from the Crank–Nicolson scheme. Also, for the SLFE algorithm,

f is the selective lumping coefficient for the bottom friction term in Eq. (2) and e is the selective lumping coefficient for all other terms that result in a mass matrix during the FE process. (We separated the lumping on the bottom friction term in order to allow independent study of this term.) These two coefficients take on values between 0, unlumped, and 1, fully lumped. For the dispersion relationships in Table 3, the GWC algorithm again results in three roots, where only two correspond to physical waves. As in the previous table, we have only reported the cubic polynomial in x0 . Note also that a dispersion relation cannot be derived for the SLFE algorithm due to the two-step temporal discretization (i.e., it is impossible to write the algorithm in a time-independent form). The coefficients A and B, which appear in this table, are the same as reported in Table 2.

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

655

Table 3 Relations for temporal frequency from analytical dispersion analysis h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii x ¼ 12 is 4ghr2  s2

Continuum

PFE

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x0 ¼ 2i s s2 þ 4ghðA  2Þ=Dx2 " # rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i2 B x0 ¼ 2i s s2 þ 36gh Dx2 3Aþ2B

GWC

i

SFD

3

2 ð4 þ AÞ x0 þ i6 ðG þ sÞð4 þ AÞ x0 þ 6i Gsð4 þ AÞ hgs

ihg hg 0 2 þ Dx ¼0 2 ð2  AÞx þ Dx2 ð2  AÞ þ 4 Dx2 ðs  GÞB  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x0 ¼ 2i s  aðA  2Þ=Dx s2 þ ðaB=DxÞ2 6

Low-order FV

4. Numerical propagation analysis tools Some discretization schemes evolve with time and cannot be written in the closed form required for analytical analysis techniques. One such algorithm is the high-order FV, which requires slope limiting procedures in order to prevent oscillations in the solution. Many standard limiters and hybrid variations on these exist, but currently, there is no clear method for choosing the ‘‘best’’ slope limiter. Often, this choice is based on a trial-and-error process, where the first successful limiter is used. In an effort to systematically study the properties of some of the more popular limiters as applied to the FV method, we have developed numerical analogs to traditional Fourier and dispersion analyses. Kinnmark and Gray [11] developed a ‘‘2 Dx-test’’ for the specific examination of spurious 2 Dx waves, but we desire to examine the propagation characteristics of the full spectrum and therefore developed the following procedure. We also acknowledge the work of Vichnevetsky and Shieh [17], who applied a related procedure to the scalar transport equations. The general theory behind the following numerical technique is to use the

1

numerical discretization scheme under study to propagate a single wave with known amplitude and wavelength exactly one wavelength forward in time and then examine the final location and amplitude of the wave. This information is then translated into phase and damping characteristics (for the equivalent Fourier analysis) and dispersion characteristics (for the equivalent dispersion analysis). The procedure is then repeated for each wavelength with the final product being numerically generated phase error, damping, and dispersion curves. In this way, the propagation characteristics of the linear algorithm can be obtained one wavelength at a time. As mentioned in the introduction, advantages of such an approach are the ability to study full discrete schemes, including boundary effects, the generality of working without a closed-form equation, and a direct relationship between the propagation behavior and actual simulation output. For simplification we only consider domains of constant bathymetric depth in this study, but with slight modifications the general procedure can be transferred to more complicated bathymetric domains.

Dnumeric

bI

0.75

bF

bI - initial peak amplitude bF - final peak amplitude xI - initial peak location xF - final peak location

Amplitude (u or ζ)

0.5 0.25 0

xI

- 0.25

xF

Dnumeric = xF - xI, is the distance that the wave has propagated

- 0.5 - 0.75 50000

100000

150000

200000

250000

300000

Horizontal position, x (meters) Fig. 2. Schematic of wave propagation properties ðLn =Dx ¼ 75Þ. Initial wave is represented by discrete points (  ) and final wave is represented with solid line (––).

656

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

A schematic representation of the general procedure is shown below in Fig. 2, where the illustrated curves were generated for the low-order FV algorithm using our method. The first wave (discrete points) is the initial wave introduced into the center of the domain at time, t ¼ 0, and the second wave (solid line) is the propagated wave after Nn time steps of length Dt. Note that the initial wave has been shifted such that the peak lies directly on a node point (to facilitate wave identification). Correspondingly, note that the final peak does not, in general, fall directly on a grid point (a further consequence of discrete numerical schemes representing continuous waves). We have chosen to use a local maximum value from the discrete grid points as an approximation for the peak instead of using interpolation to estimate the peak location. Note also that the trailing portion of the second wave is not identically zero and that there is some wiggle in the leading portion as well; this is an indication of the inherent problems in numerical discretizations. For reference, the notation used in this figure will be carried throughout the following derivations. As mentioned above, there is some error introduced by not knowing the exact location of the final peak, since it most likely lies between grid points. In order to characterize these errors, we examined two other methods of tracking the wave position; namely, a center of mass (COM) approach, by which we numerically integrated the initial and final waves and used first moment principles to determine the ‘‘center’’ of the wave and a least squares harmonic (LSH) analysis. For the LSH analysis we ‘‘fit’’ the final waveform with a wave of the appropriate wavelength and determined the amplitude and phase shift in a least squares sense. In general, we found that the errors of the numerical scheme relative to the analytic propagation behavior were more pronounced when a COM approach was used and that the LSH approach occasionally produces smoother results in the lower wavelengths, but in other cases it is just as noisy as peak tracking. Ultimately we chose to use peak tracking to determine the phase characteristics for this paper.

4.1. Numerical Fourier analysis Step 1: Initialize with a single wave of known amplitude, wavelength, and position. The domain size will depend upon the wavelength being simulated and must be large enough that the wave will not run into the ‘‘land’’ boundary (and reflect back into the domain); a domain length of three wavelengths is sufficient when the wave is initialized at the open boundary and four wavelengths is sufficient if the wave is initialized near the center of the domain. The wave may be shifted, such that the initial peak falls on a grid

point, for ease of comparison. In this study we use a simple sine wave represented as sin

2pðx  /Þ Ln

ð24Þ

where / is an appropriate shift. Only Ln þ 1 grid points are initialized, such that a single wave is within the domain, and the remaining grid points are set to zero. In general, there are two initialization options available for introducing the wave into the domain: (1) use a single wave as the initial condition at t ¼ 0, or (2) use the open ocean boundary condition to force the domain with the appropriate wavelength for one full period. In developing our procedure, it was not clear at the outset which of these two methods to introduce the wave into the domain was best. Thus, we undertook a systematic study in an effort to determine the importance of the initialization procedure and any possible boundary effects on the final phase results. We examined the following methods for initializing the starting wave: Method 1. Initialize wave near the open ocean boundary. Method 2. Initialize wave near the center of domain. Method 3. Use elevation forcing of appropriate wavelength (period) at the open ocean boundary to bring the wave into the domain. For Methods 1 and 2, the elevation and velocity profiles were identically initialized with one wavelength of the wave at time t ¼ 0 and then this wave was allowed to propagate for one full period. For Method 3, we used the open ocean boundary condition to force the system with a wave of the desired wavelength for one full period, then turned the forcing off and allowed the wave to propagate for another full period. In order to compare these initialization schemes, we used both qualitative and quantitative measures. Qualitatively, we compared the numerically generated propagation characteristics (damping, phase, celerity, and dispersion) with plots of the analytical results from Section 3.3 to determine if any of the initialization methods were grossly in error (i.e., did not capture the general shape and trend of the analytical plots). Quantitatively, we used an L2 norm for the relative error in the phase and damping characteristics of each algorithm. In general, Methods 1 and 2 are not significantly different from one another and Method 3 more closely captures the dispersion behavior than it does Fourier behavior. After careful examination of the qualitative and quantitative measures, we determined that there was not one initialization technique that would characterize all discretization schemes equally well. We thus recommend using Method 2 for all discretization schemes except for the GWC, which must use Method 3, when one is in-

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

terested in Fourier behavior. For dispersion behavior, Method 3 is appropriate for all discretization schemes. All results in future sections will be generated following these recommendations. Step 2: Execute SWE algorithm and propagate wave for Nn timesteps. Once the wave has been introduced to the domain, the next step is to propagate the wave forward approximately one wavelength by running the SWE solution algorithm for enough time steps so that the wave would travel one wavelength analytically. We know from the analytical Fourier analysis in Section 3.1 that the wave should travel exactly one wavelength in Nn ¼ Ln =ðc DtÞ time steps where c is defined in (14) and depends upon the wavelength. Of course, within the framework of any computational procedure the machine cannot compute the results for fractional time steps and the mathematical ceiling (or floor) of Nn must be used for the number of time steps in the actual code. Although this does not introduce a large error into the numerical phase scheme, it is easily corrected in the analytical wave solution and should be taken into consideration before computing the complex propagation factor T . In order to account for this additional propagation time in the analytical results, we define the corrected propagation time as Nn0 ¼ ceiling Nn and note that the analytical solution predicts a final amplitude of 0

jkj ¼ esNn Dt=2

ð25Þ

and a distance travelled of Danalytic ¼ Nn0 c Dt

ð26Þ

for an initial wave of amplitude 1. Using Nn0 instead of Nn in Eq. (18) allows an equal comparison between the numerical and analytical results. Step 3: Compute the final position and amplitude and compare to the initial wave from Step 1. From Step 1, the exact location and magnitude of the initial peak are known. Determining the location of the peak after propagation is somewhat more difficult since it will most generally not lie on a grid point. However, a local maximum can be determined from the discrete data and, for the larger physical wavelengths, little error arises from not knowing the exact location of the peak. Since we have moved the wave forward approximately one wavelength, the propagation properties are most naturally reported on a per wavelength basis, which also allows comparison with the complex propagation factors computed in Section 3.3. Referring to Fig. 2, the propagation properties for each wavelength are then computed as magnitude of T ¼ jT j ¼ ðbF =bI Þ=jkj

ð27Þ

657

and phase of T ¼ ðDnumeric  Danalytic Þ2p=Ln

ð28Þ

where bI is the initial peak amplitude, bF is the final peak amplitude for u or f (approximated by the local maximum as discussed above) after the wave travels Nn0 time steps and the numeric distance traveled is quantified as Dnumeric ¼ xF  xI

ð29Þ

where xI and xF are the initial and approximate final peak locations. From the expression for the phase, it is apparent that if the change in peak position is exactly equal to the analytic distance, then the phase error is zero; and if the peak travels more (less) than the analytic distance, then there is a phase lead (lag). Similarly, if the final peak amplitude is greater (less) than that of the analytic wave, then the numeric wave is amplified (damped) relative to the analytic solution. Similarly, one can examine the relative wave speed or celerity of the numerically propagated wave by computing the ratio of distances traveled in Nn0 time steps, since the time is the same for both the numerical and analytical waves. This ratio, represented as c0 Dnumeric ¼ c Danalytic

ð30Þ

where c0 is the numerical wave speed, is also an indicator of the phase properties. When the numerical wave speed is equal to the analytical wave speed, indicated by a ratio of 1.0, the wave has propagated at the same speed as predicted analytically and there will be no phase error. In the same way, a ratio greater (less) than 1.0 would indicate a phase lead (lag) in the solution. Step 4: Repeat Steps 1–3 over range of wavelengths to generate plots. Propagation data is output for each wavelength and discrete plots of phase and amplitude for the complex propagation factor versus wavelength are generated. 4.2. Numerical dispersion analysis Assuming constant bathymetry and bottom friction, and hence constant phase speed, a connection between the propagation phase from Fourier analysis and the dispersion relationship can be established through the definition for phase speed, given as c ¼ x=r. By rearranging this definition and estimating the phase speed from the waveÕs final position after Nn0 time steps, an expression for x0 can be established as    2p xF  xI x0 ¼ rc0 ¼ ð31Þ Ln Nn0 Dt Thus, a discrete dispersion plot can be generated one wavelength at a time.

658

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662 Table 4 Parameters used in verification study

5. Results 5.1. Verification of numerical tools

0. 8 0. 6 0. 4 PFE

0

| T | (dimensionless)

phase of T (degrees)

1

0. 2 20

40

60

Ln /∆x

80

0. 8 0. 6 0. 4 GWC

0

20

40

60

80

300 200 100 0 - 100 - 200 - 300

100

1

0. 2

three discretization schemes captures the full range of possible propagation behaviors. The Fourier characteristics for all three discretization schemes are shown below in Fig. 3 and the dispersion curves for only the PFE and GWC algorithms are shown in Fig. 4. For the PFE discretization, we notice that the phase of T is well-characterized by our numerical technique, as evidenced by the close agreement between the discrete data points and the analytical curve. We also notice that the damping behavior for this discretization scheme approaches the analytical damping curve asymptotically but contains some wiggles. This is due, in part, to the difficulty expressing a wave with so few data points in the lower wavelength range. It is also a manifestation of the PFE discretization schemeÕs tendency to produce wiggles in the solution. From the phase behavior, we

phase of T (degrees)

| T | (dimensionless)

A comparison of the analytical propagation characteristics, as computed in Section 3.3, to the numerically generated results from Section 4 is essentially a comparison of the theoretical propagation behavior predicted by the analytical analyses to the actual propagation behavior one might expect in practice. Hence, one should not expect an exact match between the numerical scheme and the analytical results (e.g., the former does not consider boundary effects). Nevertheless, using the methods as described herein, our numerical techniques are able to predict the general shape of the propagation curves for all of the discretization schemes under study. The parameter values used in this study are summarized in Table 4. As discussed in Step 1 of Section 4.1, we rigorously compared the numerical propagation results to the analytical characteristics for the five discretization algorithms in Table 1. For the purpose of illustration, we will present the results for only three of these: the PFE with leap-frog time-stepping (PFE), the GWC, and the low-order FV approaches. The general behavior of these

g ¼ 9:81 m/s2 h ¼ 10 m Dt ¼ 1:0 s s ¼ 0:0001 s1 G ¼ 0:001 s1 Dx ¼ 1000 m

PFE

0

20

40

60

80

0

20

40

60

80

300 200 100 0 - 100 - 200 - 300

100

0. 6 0. 4 0. 2

FV

0

20

40

60

Ln /∆x

80

100

100

Ln /∆x phase of T (degrees)

| T | (dimensionless)

0. 8

100

GWC

Ln /∆x 1

Ln /∆x

300 200 100 0 - 100 - 200 - 300

FV

0

20

40

60

80

100

Ln /∆x

Fig. 3. Fourier damping and phase error characteristics for the PFE, GWC, and FV discretization methods. Discrete points (  ) generated using our numerical techniques and solid lines (––) generated from analytical results in Table 2.

0.03 0.025 0.02 0.015 0.01 0.005 0 0

temporal frequency, ω

temporal frequency, ω

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

PFE

0.001 0.002 wave number, σ

0.003

0.03 0.025 0.02 0.015 0.01 0.005 0 0

659

GWC

0.001 0.002 wave number, σ

0.003

Fig. 4. Dispersion characteristics for the PFE and GWC. Discrete points (  ) generated using our numerical techniques, solid lines (––) generated from analytical results in Table 3, and dashed lines (- - -) generated from continuum results in Table 3.

would expect that the PFE discretization would produce wiggles, on the order of 2–3 Dx, as evidenced by the phase error of 2p with a corresponding magnitude of 1.0, which indicates that these spurious waves will not be damped. For the GWC discretization, the numerical technique captures the magnitude and phase behavior in the longer wavelengths nearly exactly but has some difficulty in the 2–5 Dx range. While there is a slight phase lead evidenced in the low wavelength range of the GWC data, there is also sufficient damping for these same wavelengths. Again for the FV discretization, the numerical technique is able to capture the general trend of the phase and damping characteristics. We notice that the FV exhibits a phase lag of 2p at the 2 Dx wavelength, but the damping characteristics indicate that this wave will be completely damped out (i.e., jT j is zero at this wavelength). However, we observe that there is considerable damping in the long-wavelength physical waves as well and it is seen that the method is overly dissipative as jT j never reaches unity. This result confirms our observations in Section 2 where we noted that the introduction of the second-order spatial derivative adds excessive diffusion into the solution. However, this diffusion can be controlled by using higher-order interpolates and introducing appropriate limiting functions, as we will show in Section 5.2. Comparison of the analytical and numerical dispersion plots in Fig. 4 shows similar behavior to the Fourier results. Here also the lower wavelengths (i.e., higher wavenumbers––rightmost points on the x-axis) deviate most from the analytical curves while the higher wavelengths are nearly exact. However, despite the inexactness of the data points in the 2–5 Dx range, the analytical behavior is still captured and the discrete data points are able to predict a folded curve (for the PFE), indicating aliasing of two wavenumbers for a single temporal frequency. The monotonic curve for the GWC scheme is also successfully predicted. It appears then, that although our scheme does have certain limitations, it is able to correctly capture the overall propagation behavior and dispersion characteristics for a range of possible behaviors (i.e., monotonic versus folded dispersion).

In response to the observed errors in the low wavelength field of the spectrum, we considered alternative methods for capturing the phase behavior in this range. In particular, we examined regression curves through the discrete data points in the 2 Dx–20 Dx range, in order to smooth out the oscillations, with the hope of obtaining a more accurate approximation. This did remove the oscillations in the discrete propagation data and proved to generate more accurate dispersion results in the low wavelength data range; however, we did not see significant enough improvements to warrant such an approach for the Fourier results (phase and damping). Consequently, we decided to use a regression approach for the dispersion curves within this data range and then use the raw data for the higher wavelengths. (Raw data was used for the entire phase and damping curves.) The graphs shown in Fig. 4 were generated in this manner. Additionally, the analytical Fourier and dispersion analyses generate a single graph from the determinant of the equation system, while two graphs, one each for elevation and velocity, are generated with our numerical techniques. Upon examination of the results, however, we noticed that the elevation and velocity results either enveloped the analytical curves or were nearly identical. Thus we chose to report the average of the individual elevation and velocity curves in order to have a single curve for each SWE discretization scheme, as in the analytical analysis technique. Again, the graphs shown in Figs. 3 and 4 were generated using this procedure.

5.2. Case study––slope limiters within FV method We now apply our numerical propagation analysis techniques to higher-order (piecewise linear) FV discretizations, wherein the slope limiting procedures prohibit a traditional analytical analysis. In particular, we have examined the three slope limiters presented earlier in Section 2: minmod, vanLeer and Superbee. The propagation curves, as generated by our numerical techniques, are presented in Fig. 5; notice that all of the phase and magnitude data is discrete since there are no

1 0. 8 0. 6 0. 4 0. 2 A

0

temporal frequency, ω

0 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0

20

40

60

Ln /∆x

80

100

phase of T (degrees)

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

| T | (dimensionless)

660

100 50 0 -50 -100 -150

B

0

20

40

60

Ln /∆x

80

100

Minmod vanLeer Superbee C

0

0.001

0.002

wave number, σ

Continuum

0.003

Fig. 5. Propagation characteristics of higher-order FVM methods using slope limiters.

analytical results to present, but that the continuum dispersion behavior has been plotted for reference. The minmod slope limiter scheme exhibits a slight phase lag in the low wavelengths, but these would be damped out, as shown in the magnitude plot. However, this scheme appears to be overly dissipative as the magnitude only reaches a maximum value of approximately 0.9, implying that the long physical waves will also be damped, despite their correct phase. The vanLeer limiting scheme exhibits an asymptotic phase lag where the long wavelength, physical waves are approximately 25° out of phase and reach a relative magnitude of only 0.8 of the continuum waves. This scheme would not appear to be a good choice in practice. In contrast to these two limiting methods, the Superbee exhibits nearly perfect phase behavior and much less damping, as it approaches the desired magnitude of 1.0 within the first 20 Dx wavelengths. This is within the common range used in application (i.e., in practice we do not usually simulate waves smaller than 20 Dx but rather adjust grid spacing accordingly). Note that even with the higher-order FV method, Eqs. (4) and (5), the additional second-order spatial derivative introduced via the Riemann problem still exists (see Section 2 for discussion). However, if we are judicious about which slope limiting algorithm is used to reconstruct the slopes, we only add enough diffusion to keep the algorithm stable (‘‘controlled’’ diffusion). From the results in Fig. 5, it is apparent that the Superbee class of limiters is more successful than the other two at controlling the amount of dissipation that is introduced by the second-order spatial derivative. Thus, with our numerical technique, we now have a tool to judge slope limiters for their ability to damp noise without damping the physical waves as well.

Meanwhile, the dispersion results (Fig. 5c) indicate that the minmod and Superbee limiters produce monotonic dispersion relationships, while the vanLeer limiter begins to fold in the low wavelengths, which is consistent with the significant lag in this region of the phase plot. Considering both Fourier and dispersion results, we find that the Superbee class of limiters is the most promising and thus warrants further study in higher dimensions. A comparison of 1D simulation results using the Superbee and vanLeer limiters with the commonly used GWC method is shown below in Fig. 6 and verifies the phase characteristics shown in Fig. 5. The simulation parameters are the same as in Table 4, and we used the main lunar semidiurnal tide, period ¼ 12.42 h, with a 1.0 m amplitude as the open-ocean boundary condition. The curves shown are an expanded view of the final peak after three full periods of simulation time. Notice that the peak of the vanLeer wave occurs at a later time and is thus lagged, as predicted. The predicted values given in the table in Fig. 6 were taken from the physical wave range of the phase plots in Fig. 5a, and the measured values were taken from the observed simulation behavior in Fig. 6. Notice that the actual phase errors seen in the simulation are predicted quite closely by the numerical phase technique. 5.3. Case study––lumping parameters in the SLFE method As a second example to show the utility of our numerical techniques, we also studied the effect of the lumping parameter in the SLFE method. After observing the analytical Fourier behavior of the SLFE dis-

Elevation amplitude (m)

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

661

1.22 1.215

Phase error in degrees (+ lead, - lag) relative to the GWC method

1.21

Measured

Predicted

1.205 1.2 1.195

vanLeer Superbee GWC 101500

102000

102500

Superbee

7.9

vanLeer

-25.8

5.3 -27.0

103000

Time (sec) Fig. 6. Comparison of simulation results for three discretization schemes after three full periods of the M2 tide. Measured errors are from this figure and predicted errors are from analysis results shown in Fig. 5.

cretization scheme, as characterized by the equations in Table 2, we discovered that the second lumping parameter, f , does not significantly alter the solution. (Recall that f is the lumping parameter for the bottom friction term in Eq. (2).) We examined values of bottom friction ranging from 105 to 103 s1 and held e constant, while allowing f to vary from 0 (unlumped) to 1 (fully lumped). There was no change in phase or amplitude behavior. Given that this second lumping parameter has no effect, we chose to set e and f equal for all further studies. (Recall that the dispersion behavior of the SLFE discretization scheme cannot be studied analytically due to the split-level time marching algorithm.) Kawahara et al. [10] note that in applications, the lumping parameter, e, should be between 0.8 and 0.95 and that the maximum stable time step should be used in order to reduce the numerical damping; we tested a parameter range between 0.75 and 1.0. Other than the time step, Dt, which was set to the maximum stable value for each value of the lumping parameter, all other parameters are as in Table 4. In Fig. 7 we show the continuum dispersion relation and the numerically generated dispersion curves for this range of lumping parameters. Observe that for values of the lumping parameter outside of the suggested range (i.e. for e ¼ 0:75

temporal frequency, ω

0.03 0.025 0.02

continuum ef 0.75 ef 0.8 ef 0.85 ef 0.95 ef 1.0

and 1.0) the dispersion curves begin to fold away from the continuum dispersion relation. Additionally, we note that using a time step significantly less than the maximum stable step for any value of e results in a folded dispersion curve where even the long wavelength physical waves are not propagated correctly. Thus our numerical techniques are able to verify not only the suggested range of parameter values (obtained via trial and error), but also that this algorithm is extremely sensitive to the time step.

6. Conclusions We have presented numerical analogs to traditional Fourier and dispersion analyses and shown that they are able to capture the overall propagation behavior for a variety of shallow water discretization schemes. They also provide a structured framework in which to study higher-order algorithms and other schemes that do not lend themselves to traditional analysis techniques. For the higher-order FVM discretization schemes, we note that the Superbee class of slope limiters has the most promising behavior and warrants further study in higher spatial dimensions. We also note that the practical range for the lumping coefficient in the SLFE discretization that was suggested by Kawahara et al. [10], 0.8–0.95, is verified with these numerical tools, provided that the maximum stable time step is used.

Acknowledgements

0.015 0.01 0.005 0.0005

0.001 0.0015 0.002 0.0025 wave number, σ

0.003

Fig. 7. Dispersion curves for SLFE lumping parameter study. The continuum line was generated from Eq. (22).

Financial support for this research was provided, in part, by the National Science Foundation under contract ACI-9623592 and an NSF Graduate Fellowship, and the Graduate College of the University of Oklahoma. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the funding agencies.

662

C.M. Szpilka, R.L. Kolar / Advances in Water Resources 26 (2003) 649–662

References [1] Alcrudo F, Garcia-Navarro P. A high-resolution Godunov-type scheme in finite volumes for the 2D shallow water equations. Numer Meth Fluids 1993;16:489–505. [2] Atkinson JH, Westerink JJ, Luettich RA. Two-dimensional dispersion analysis of finite element approximations to the shallow water equations. Int J Numer Meth Fluids, in press. [3] Bell JB, Shubin GR. Higher order godunov methods for reservoir simulation. Exxon Production Research Company Technical Report EPR.12PR.84, Houston, 1984. [4] Causon DM, Ingram DM, Mingham CG, Yang G, Pearson RV. Calculation of shallow water flows using a cartesian cut cell approach. Adv Water Resour 2000;23(5):545–62. [5] Chippada S, Dawson CN, Martinez ML, Wheeler MF. A Godunov-type finite volume method for the system of shallow water equations. Comput Meth Appl Mech Engng 1998;151(1– 2):105–29. [6] Foreman MGG. An analysis of the ‘‘wave equation’’ model for finite element tidal computations. J Comput Phys 1983;52:290–312. [7] Foreman MGG. A two-dimensional dispersion analysis of selected methods for solving the linearized shallow water equations. J Comput Phys 1984;56:287–323. [8] Gossard CM, Kolar RL. Phase behavior of a finite volume shallow water algorithm. In: Bentley LR et al., editors, Proceedings of the XIII International Conference on Computational Methods in Water Resources, Calgary, 2000. p. 921–8. [9] Gray WG, Lynch DR. Time-stepping schemes for finite element tidal model computations. Adv Water Resour 1977;1(2):83–95.

[10] Kawahara M, Hirano H, Tsubota K, Inagaki K. Selective lumping finite element method for shallow water flow. Int J Numer Meth Fluids 1982;2:89–112. [11] Kinnmark I, Gray WG. The 2 Dx-test: a tool for analyzing spurious oscillations. Adv Water Resour 1985;8:129–35. [12] Kinnmark I. The shallow water wave equations formulation analysis and application. Lecture notes in engineering, vol. 15. Berlin: Springer-Verlag; 1986. [13] Leendertse JJ. Aspects of a computational model for long period water-wave propagation. Rand Corporation Technical Report RM-5294-PR, Santa Monica, 1967. [14] Lynch DR, Gray WG. A wave equation model for finite element tidal computations. Comput Fluids 1979;7(3):207–28. [15] Platzman GW. Some response characteristics of finite-element tidal models. J Comput Phys 1981;40:36–63. [16] Sweby PK. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J Numer Anal 1984;21(5): 995–1011. [17] Vichnevetsky R, Shieh YS. On the numerical method of lines for one dimensional water quality equations. Department of Computer Science Technical Report #20, Rutgers University, 1972. [18] Walters RA. Numerically induced oscillations in finite element approximations to the shallow water equations. Int J Numer Meth Fluids 1983;3:591–604. [19] Westerink JJ, Luettich RA, Wu JK, Kolar RL. The influence of normal flow boundary conditions on spurious modes in finite element solutions to the shallow water equations. Int J Numer Meth Fluids 1994;18:1021–60.