Instability in Leapfrog and Forward-Backward Schemes by Wen-Yih Sun

1. Instability in Leapfrog and Forward-Backward Schemes by Wen-Yih Sun Department of Earth and Atmospheric Sciences Purdue University West Lafayette, IN. 47907-2051, USA and Department of Atmospheric Sciences National Central University Chung-Li, Tao-Yuan, 320, Taiwan E-mail: wysun@purdue.edu

2. 1. Introduction In the linearized shallow water equations, the forward-backward in time and the leapfrog schemes can be unstable for 2x waves because of the repeated eigenvalues when Courant number (Co) is 0.5 in the leapfrog scheme (LF) or Courant number is 1 in the forward-backward schemes (FB) in C-grid because of the existence of repeated eigenvalues. The weak instability of the LF for a simplified wave equation at 2x wave and Co =0.5 has also been discussed by Durran (1999). Here we provide a rigorous proof of this weak instability for FB scheme at Co=1 for 2-dx waves. Usually, diffusion terms are added to control the shortwave instability in numerical simulation of wave equations. However, It is found that the instability will be amplified and spread to the longer waves if the diffusion terms are added in both schemes. On the other hand, Shuman smoothing can be applied to control the instability for both schemes in the shallow water equations. Numerical Simulations of dam-break and vortex-merge from the 2D shallow water equations using a new semi-implicit scheme will also be presented.

3. 2. Numerical Schemes & Eigenvalues for 1-D linearized shallow water equations H is the mean depth, g is gravity, h and u are depth perturbation and velocity,  is viscosity,  =0 or 1 In staggered C-grids, the FB scheme become If a wave-type solution at the nth time step is assumed and

4. The eigenvalue of FB is , , , For simplicity, we set g=1, H=1, and x =1. Without viscosity (i.e., =0), the FB is neutrally stable, =1, as long as the Courant number However, the FB has repeated eigenvalues = -1 . The generated eigenvector x2 can be An eigenvector corresponding to 1= -1 is and Because:

5. If we define a matrix , then, we can have which is in Jordan block form and since Hence, the FB becomes weakly unstable, the magnitude linearly increases with time-step and with 2t oscillation

6. . the LF scheme are The eigenvalue of LF is When Co=0.5, (i.e., R=1), and S2=0, (i.e. =0), The LF scheme has the repeated eigenvalues too, and Hence, it also becomes weakly unstable

7. 3. Results with viscosity Diffusion has been frequently applied to control the noise of the short waves in equations, because the amplification factor is when the forward-in-time and centered-in-space scheme is applied to the diffusion equation (Sun 1982) Growth rate for FB at Co=1.0 Growth rate for FB at Co=0.8

8. Growth rate for LF at Co=0.4 Growth rate for LF at Co=0.5

9. Magnitude of (dashed line) and (solid line) as function of nth time step, where =0.91 for =0.15

10. 4. Numerical Simulations of 2D Dam-Break and Vertex-Merge where =hu and =hv.

13. Shallow Water simulation of dam break (mass conserved) (a) Vertical cross sections of water depth at t=1, 3, 10, 15 and 25s from the 4th order scheme with α=.2. dt=0.04s, dx=1m, (color) and (b) Oh’s the planar regular hexagonal MPDATA Shallow Water Equations model. Domain size is 300m, dx = 0.5m, dt = 0.02s (black).

14. Height t=0.4s Velocity t=0.4s Comparison between current (in color) & Todd’s Riemann Solver, in black) simulations Velocity t=4.7s Height t=4.7s

15. Vortex-Merge in Shallow Water: PV(Potential Vorticity)=(VH t=40 t=0 t=20 Height t=40 t=20 t=60

16. Shallow Water Potential Vorticity=(VH t=100 t=60 t=80 Height t=80 t=60 t=100

18. Surface streamlines after 10h NTU-Purdue WRF

19. Surface streamlines after 20h NTU-Purdue WRF x

20. WSMR 10 m wind field 0800 January 25, 2004 Hydraulic Jump/ Trapped Lee Waves Wind Vectors: WRF NH Model Red Arrows: Observed Winds 199x199 grid 1 km grid spacing. This plot shows the central part. Color contours show the model terrain in m asl

21. WSMR 10 m wind field 0800 January 25, 2004 Hydraulic Jump/ Trapped Lee Waves 201x201 grid 1km Grid spacing. This plot shows the central part of the model domain Wind Barbs: Purdue/NTU NH Model Dark Arrows: Observed Winds Color contours show the model terrain in m asl. Grid number west to east direction

22. (a) wind inWhite Sands after 4-hr integration, (dx=dy=2km, and dz=300m). initial wind U= 5 m/s; (b) Streamline (white line) and virtual potential temperature (background shaded colors) at z=1.8km, warm color (red) indicates subsidence warming on the lee-side, and cold (blue) color adiabatic cooling on the windward side of the mountain.

23. Summary • 1. The 2nd-order schemes are weakly unstable for 2x wave when the Courant • number is 0.5 for LF and 1.0 for FB in shallow water equations. • Adding diffusion terms make both scheme more unstable, instability extend to • 3 x and/or 4 x waves, and Co=0.4 for LF and 0.8 for FB. • Shumann smoothing can be applied to control weakly instability in both schemes. • Numerical simulations of dam break and vortex-merge using a very weak • smoothing in finite-volume difference schemes are also presented. • Reference: • Sun, W. Y., 2010: Instability in leapfrog and forward-backward Schemes. Mon. Wea. Rev. , 138, 1497–1501. • Sun, W. Y., 2011: A Semi-Implicit Scheme Applied to Shallow Water Equations and Dam Break. (Computers & Fluids: journal homepage: www.elsevier.com/locate/compfluid)