i
 
 
CRWR Online Report 08-04 
 
 
A Time-Centered Split for Implicit Discretization of Unsteady 
Advection Problems 
 
 
by 
Shipeng Fu, B.S.; M.S. 
 
and 
 
Ben R. Hodges, Ph.D. 
 
 
 
The University of Texas at Austin 
May 13, 2008 
This document is available online via World Wide Web at 
http://www.crwr.utexas.edu/online.html  
 ii
 
 
 
 
 
 
Copyright 
by 
Shipeng Fu 
2008 
 iii
 
 
 
 
Acknowledgements 
 
I would like to begin by expressing my gratitude to my advisor Dr. Ben Hodges.  I 
am thankful for his guidance and support, and I appreciate the freedom he has given me 
to explore and pursue my own research interests.  I would also like to thank Dr. Graham 
Carey, Dr. Randall Charbeneau, Dr. Spyros Kinnas and Dr. Daene Mckinney for serving 
on my committee. 
The faculty members, facilities, and the world class program of EWRE have made 
my research possible.  I am grateful to Dr. Kinnas and Dr. Carey for their input in the 
development of the numerical methods in this dissertation. I have gone to Dr. Liljestrand, 
Dr. Lawler and Dr. Kinnas many times for help and advice, and they have always been 
happy to help. I want to thank them for their encouragement and advice. I would like to 
thank Dr. Charbeneau for giving me the opportunity to work with lab experiments (so I 
could finally play with the real water!). I would like to thank Sharon Bernard, Marcy 
Betak, and Susan Swanson-Cartwright at the CRWR Office for helping with all my 
administrative paperwork and questions. 
During my Ph.D study at UT, my friends in Austin have been extremely helpful 
and supportive. I could not finish this without their friendship. I want to particularly 
thank Paula Kulis, Li-Jung Chen and their wonderful husbands for feeding me meals and 
listening to my whines. I also want to thank Paula for many discussions in fluid 
mechanics and beyond. I want to thank Yihsiang Yu and Vimal Vinayan for their 
insightful discussions on numerical methods and wave mechanics. I want to thank Becky 
Teasley and John Allen, Rebekah and Nate Johnson, and Jeremy Seibert for graciously 
reading and providing valuable feedbacks on my dissertation. I want to thank all my 
friends in EWRE and CRWR. Our lunch topic at CRWR is always intellectually 
stimulating, and sometimes wild. I also want to thank all the members in the BUSTED 
(Best Underpaid Students Turning into Environmental Doctors) club. I really benefited 
from the support from everyone. I will miss you girls. I want to thank my running buddy, 
Susan De Long, for many interesting discussions about Yoga and philosophy.  
I am thankful to my parents. I thank them for teaching me good values, and 
encouraging me to think independently since I was very young. My diligent father and 
curious mother have always been my role models and inspiration.   
My deepest thanks go to my best friend and husband, Yunzhi Yang. I am really 
fortunate to have such a wonderful life partner accompanying me for this journey. I want 
to thank him for truly understanding me, encouraging me to strive for excellence, and 
bringing life into perspective. 
 iv
Abstract 
 
Environmental flows (e.g. river and atmospheric flows) governed by the shallow 
water equations (SWE) are usually dominated by the advective mechanism over multiple 
time-scales. The combination of time dependency and nonlinear advection creates 
difficulties in the numerical solution of the SWE. A fully-implicit scheme is desirable 
because a relatively large time step may be used in a simulation. However, nonlinearity in 
a fully implicit method results in a system of nonlinear equations to be solved at each 
time step. To address this difficulty, a new method for implicit solution of unsteady 
nonlinear advection equations is developed in this research.  This Time-Centered Split 
(TCS) method uses a nested application of the midpoint rule to computationally decouple 
advection terms in a temporally second-order accurate time-marching discretization. The 
method requires solution of only two sets of linear equations without an outer iteration, 
and is theoretically applicable to quadratically-nonlinear coupled equations for any 
number of variables.  
To explore its characteristics, the TCS algorithm is first applied to one-
dimensional problems and compared to the conventional nonlinear solution methods.  
The temporal accuracy and practical stability of the method is confirmed using these 1D 
examples.  It is shown that TCS can computationally linearize unsteady nonlinear 
advection problems without either 1) outer iteration or 2) calculation of the Jacobian. A 
family of the TCS method is created in one general form by introducing weighting factors 
to different terms. We prove both analytically and by examples that the value of the 
weighting factors does not affect the order of accuracy of the scheme. In addition, the 
TCS method can not only computationally linearize but also decouple an equation system 
of coupled variables using special combinations of weighting factors. Hence, the TCS 
method provides flexibilities and efficiency in applications. 
 v
Table of Contents 
List of Tables ........................................................................................................ vii 
List of Figures...................................................................................................... viii 
Chapter 1 Introduction .............................................................................................1 
1.1 Background...............................................................................................1 
1.2 Challenge and Motivation of the Numerical model..................................3 
1.3 Objective...................................................................................................4 
1.4 Approach...................................................................................................4 
1.5 Broad Application.....................................................................................5 
Chapter 2 Review of Finite Difference Schemes for Solving the Unsteady Nonlinear 
Advection in the Shallow Water Equations....................................................6 
2.1 Explicit Method ........................................................................................6 
2.2 Semi-Implicit Method...............................................................................7 
2.3 Implicit Method ........................................................................................7 
2.4 Summary.................................................................................................11 
Chapter 3 Theoretical Development of the Time-Centered Split (TCS) Method..12 
3.1 The Theory of Time-Centered Splitting .................................................12 
3.2 Derivation of the TCS method in a 1D Advection-Diffusion Equation .16 
3.3 TCS for Coupled Momentum and Scalar Transport...............................19 
3.4 Summary.................................................................................................20 
Chapter 4 Implementation of the TCS Method in 1D Problems ...........................22 
4.1 Application of the TCS Method to the 1D Conservative Burgers? Equation
..............................................................................................................22 
4.2 Applications of the TCS Method to the 1D Non-conservative Burgers? 
Equation ...............................................................................................42 
4.3 Application of the TCS Method to a 1D Nonlinear Ordinary Differential 
Equation ...............................................................................................59 
4.4 Summary.................................................................................................63 
 vi
Chapter 5 The TCS Family Method.......................................................................65 
5.1 Derivation of the TCS Family Method ...................................................65 
5.2 Application of the TCS Family Method to the 1D Non-conservative 
Burgers? Equation ................................................................................68 
5.3 Results and Discussion ...........................................................................70 
5.4 Computational Decoupling .....................................................................78 
5.5 Summary.................................................................................................81 
Chapter 6 Application of the TCS Method to a 2D Depth Averaged Shallow Water 
Equations (SWE) ..........................................................................................83 
6.1. The 2D Depth Averaged SWE...............................................................83 
6.2 The TCS Discretized SWE .....................................................................84 
6.3 Decoupling the SWE...............................................................................87 
6.4 Characteristics of the TCS Decoupled Equation System........................89 
6.5 Numerical Tests ......................................................................................91 
6.6 Summary...............................................................................................109 
Chapter 7 Conclusions and Recommendations....................................................110 
7.1 Summary of Discussion ........................................................................110 
7.2 Conclusions...........................................................................................113 
7.3 Recommendations for Future Work......................................................113 
Appendix A Test of a Progressive Wave in an Open Boundary System.............115 
A.1 Initial and Inlet Boundary Condition ...................................................115 
A.2 Outlet Boundary Condition..................................................................116 
A.3 Results and Discussion.........................................................................119 
References............................................................................................................124 
 vii
List of Tables 
Table 5.1 The mathematical meaning of the value of 
1
? and 
2
? ...................................... 70 
Table 5.2 Weighting factors for computational decoupling the 2D Burgers? equation.... 80 
Table 6.1 Weighting factors and the solution orders for decoupled Equations (6.11) 
through (6.16). .......................................................................................................... 89 
Table 6.2 Solution procedures of the three TCS discretizations....................................... 91 
Table 6.3 Parameters of simulations of a 1D standing wave............................................ 93 
Table 6.4 Parameters of simulations of the 2D standing wave......................................... 97 
  
 viii
List of Figures 
Figure 1.1 Sulphur River (Texas) large woody debris at low-flow conditions .................. 2 
Figure 4.1 Solution of 1D Burgers? equation evolving in time for t?{0, 0.1, 0.2, 0.3, 0.6, 
1, 2, 3} using 0.05?=  along with the initial and boundary conditions of Equations   
(4.43) through (4.45)................................................................................................. 34 
Figure 4.2 Absolute error time evolution for numerical solutions of 1D conservative 
Burgers? equation for TCSF and TCSF-D using { t? =0.01, x? =1/50, ?=0.05}, 
where 
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
 and 
analytical
u (x, t)  is numerically 
calculated from Equation  (4.46) with K=30. ........................................................... 35 
Figure 4.3 
RMS
? vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and t0.3/? =?with 
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 37 
Figure 4.4 
2
L  norm vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and t0.3/? =?with 
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 38 
Figure 4.5 L
?
norm vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and t0.3/? =?with 
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 39 
Figure 4.6 
RMS
?  vs. CFL number of TCSF, RK2 and RK4 methods for solution of 1D 
conservative Burgers? equation, where 0.05?= x? =1/50, and t0.3/? =?with 
??{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 5}...................................................... 40 
Figure 4.7 Ideal operations per grid point for one time step using various nonlinear 
solution methods for 1D conservative Burgers? equation. R is the number of outer 
iterations taken by Newton method.  It is assumed that the Picard method 
convergences in R
2
 outer iterations........................................................................... 42 
Figure 4.8 Absolute error time evolution for numerical solutions of 1D non-conservative 
Burgers? equation for the TCSF and TCSG methods using 
{ }t 0.01, x 1/ 50, 0.05?= ?= ?=  where 
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
 
and 
analytical
u(x,t) is numerically calculated from Equation  (4.46) with K=30........ 54 
 ix
Figure 4.9 
RMS
?  vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50 
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 55 
Figure 4.10  
2
L  norm vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50 
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 56 
Figure 4.11 L
?
 norm vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50 
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 57 
Figure 4.12 Temporal accuracy of the TCS and RK methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50 
and t? =0.6/? with ? ?{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 5}.................... 58 
Figure 4.13 Ideal operations per grid point for one time step using various nonlinear 
solution methods for 1D Non-conservative Burgers equation. R is the number of 
outer iterations taken by Newton method. It is assumed that the Picard method 
convergences in 
2
R outer iterations.......................................................................... 59 
Figure 4.14 Time evolution of absolute errors for TCS methods applied to the ODE, 
Equation  (4.101), for t? =0.1, where 
abs model analytical
(t) y (t) y (t)?= ? . .................... 62 
Figure 4.15 RMS error from Equation  (4.117) for discrete solutions of the ODE, 
Equation (4.101), computed using a range of ? time steps. .................................... 63 
Figure  5.1 Temporal accuracy of various combinations of 
1
?  and 
2
? for solution of the 
Burgers? equation at t=0.6, where x150/? = , 005.?=  and t06./? =? 
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 72 
Figure  5.2 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=0.6, where x150/? = , 005.?=  and t06./? =? 
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 73 
Figure  5.3 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=0.3, where x150/? = , 005.?=  and t03? =?./  
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 75 
Figure  5.4 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=1, where x150? = / , 005?= .  and t1? =?/  
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 76 
 x
Figure 5.5 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=2, where x150/? = , 005.?=  and t2? =?/  
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 77 
Figure 6.1 Illustration of Arakawa C grid......................................................................... 85 
Figure 6.2 A standing wave in a rectangular basin........................................................... 92 
Figure 6.3 Simulations in the inviscid flow with different wave steepness...................... 94 
Figure 6.4 Viscous damping effects in the shallow water model ..................................... 95 
Figure 6.5 The initial 2D standing wave and its decomposed x and y direction 1D wave96 
Figure 6.6 Computational domain (not to scale) .............................................................. 97 
Figure 6.7 Normalized surface elevation contours in the computational domain at time 
4T. ........................................................................................................................... 100 
Figure 6.8 Normalized surface elevation contours in the computational domain at time 
15T. ......................................................................................................................... 101 
Figure 6.9 Normalized U-velocity contours in the computational domain at time 4T... 102 
Figure 6.10 Normalized U-velocity contours in the computational domain at time 15T.
................................................................................................................................. 103 
Figure 6.11 Normalized V-velocity contours in the computational domain at time 4T. 104 
Figure 6.12 Normalized V-velocity contours in the computational domain at time 15T.
................................................................................................................................. 105 
Figure 6.13 Surface displacement at points (0.5, 0) and (0, 0.5) simulated using three 
TCS solutions.......................................................................................................... 106 
Figure 6.14 Snapshot of the standing wave at 0.5 T. The surface displacements at the two 
monitored points (0.5, 0) and (0, 0.5) are circled. .................................................. 107 
Figure 6.15 Snapshot of the standing wave at 15 T. The surface displacements at the two 
monitored points (0.5, 0) and (0, 0.5) are circled. .................................................. 108 
Figure A.1 Initial water level, H (x, 0), in the rectangular open channel ....................... 115 
Figure A.2 Inlet boundary condition H (0, t).................................................................. 116 
Figure A.3 Shape function of viscosity........................................................................... 117 
Figure A.4 An open boundary rectangular channel with a sponge layer........................ 117 
Figure A.5 Schematic illustrations of grids .................................................................... 118 
Figure A.6 Water level at point A with/without the sponge layer.................................. 118 
 xi
Figure A.7 Water level evolutions inside the test section and sponge layer................... 119 
Figure A.8 Comparison of wave shape at time 7T and 14T for TCS solution1 and TCS 
solution 2 respectively ............................................................................................ 120 
Figure A.9 Comparison of wave shape simulated from TCS solution 1 and 2 at time 7T 
and 14T respectively............................................................................................... 121 
Figure A.10 Time evolution of a progressive wave at different wave periods from TCS 
solution 1................................................................................................................. 122 
Figure A.11 Time evolution of a progressive wave at different wave periods from TCS 
solution 2................................................................................................................. 123 
 
 1
Chapter 1 Introduction 
 
Models of river flow using mathematic tools have been developed since last half 
century (Cao and Carling 2002). Various numerical river models have been created to 
study flood control (Dutta et al. 2007; Liao et al. 2007), water allocation (Fleckenstein et 
al. 2006; Luo 2007), sediment and morphological evolution (Giri and Shimizu 2006; Le 
et al. 2006) and so on. The river flow is usually modeled at a coarse-grid scale due to the 
limitations of the computational power (Dedong et al. 2007; Shaw et al. 2005; Thouvenin 
et al. 2007). However, there exist some particular river systems such as a river with many 
large woody debris (LWD), as shown in Figure 1.1. The flow information around LWD 
is at subgrid-scale, which can?t satisfactorily be modeled at the conventional coarse-grid. 
In addition, the river flow is often simulated using two-dimensional (2D) depth 
averaged shallow water equations (SWE) as the governing equations (Chen and Peng 
2006; Nguyen et al. 2006; Weerakoon et al. 2003). However, the unsteadiness and the 
dominant advective mechanism of a river flow create difficulties in SWE based 
numerical simulation because of the combination of the ?stiffness? (i.e. the unsteady 
term) and ?nonlinearity? (i.e. the nonlinear advection term) in the momentum equations. 
Therefore, in this research, we develop a new numerical algorithm to address the 
difficulties in simulating the unsteady nonlinear advection. This numerical scheme can 
also serve as a platform to perform the new conceptual model which can integrate the 
subgrid-scale physics into a coarse-grid scale model. The detailed development of the 
conceptual model can be found in Fu and Hodges (2005).  First, let us start with the 
background information of the LWD.  
1.1 BACKGROUND  
Large woody debris (LWD) refers to woody material such as fallen tree trunks or 
root balls that become lodged in stream channels. Figure 1.1 gives an example of a river 
laden with large woody debris (LWD).  The conventional size of LWD is diameters 
larger than 0.1 m and lengths greater than 1.0 m (Keller and Swanson 1979; Andrus et al. 
1988). LWD accumulations have been historically called ?snags?, a pejorative term 
reflecting a perception of LWD as a nuisance to river navigation and efficient water use. 
River ?improvement? schemes typically involve removing LWD (Shields & Nunnaly 
1984; Gippel 1995), to improve water conveyance, rejuvenate channels, lessen the risk of 
damage to bridges, improve recreational amenity, and remove barriers to fish migration 
(Harmon et al. 1986). However, from an ecological perspective, LWD provides a stable 
substratum for microorganisms, algae and invertebrates (Minshall 1984; Brown & May 
2000; Statzner & Higler 1986). This, in turn, creates grater local productivity for higher 
trophic levels such as invertebrates and fish (Price & Lovett 2002). Furthermore, LWD 
enhances hydraulic diversity (i.e. a wide range of flow conditions), which in turn 
enhances diversity in fish habitat, e.g. providing low-velocity refuges sought by many 
fish species (Benke et al. 1985; Matthews & Hill 1980), and also providing cover from 
predators (Angermeier & Karr 1984, Everett & Ruiz, 1993), foraging habitat, and 
 2
spawning substratum - all of which vary dynamically with flow rate (Bao & Mathews 
1991; Mathews & Tallent 1997; Marzolf 1978). In short, LWD creates greater habitat 
complexity (O?Conner 1991), which should produce greater biodiversity (Gerhard & 
Reich 2000).  
 
 
 
 
 
 
Figure 1.1 Sulphur River (Texas) large woody debris at low-flow conditions 
(Photo courtesy of Texas Water Development Board) 
 
Over the past two decades, a wide variety of studies demonstrate how woody 
debris may be a dominant player in aquatic habitat. Angermeier & Karr (1984) 
selectively modified a stream by inducing and removing woody debris, subsequently 
showing a greater abundance of fish and benthic invertebrates correlated with the 
introduced debris. Benke et al. (1985) examined a low-gradient stream, showing that only 
4% of habitat surface was in LWD and yet this supported 60% of the invertebrate 
biomass and 16% of overall production. Similar results are noted in Benke et al. 1984 and 
Jacobi & Benke1991. Effects of flooding on fish movements and distribution (Hill & 
Grossman 1987; Harvey et al. 1999) have confirmed the importance of habitat 
complexity (including LWD) for population persistence. More subtle and indirect effects 
on habitat have been attributed to LWD?s modification of the flow field and water level 
(Triska & Cromack 1980; Sedell et al. 1982; Wallace & Benke 1984). Near LWD, 
reduced velocities allow retention of organic matter, building up ecologically-important 
debris dams, while water lever effects at high flow rates can influence the seed dispersal 
from riparian vegetation (Merritt & Wohl 2002). Studies drawing similar conclusions 
have been conducted over a variety of locales and river types, ranging from salmon 
spawning streams in Alaska (Dolloff 1986), low-gradient rivers in West Virginia (Lobb 
& Orth 1991) to an ephemeral river in Africa (Jacobson et al. 1999). The common theme 
is that LWD plays a varied and critical role in aquatic habitat for rivers and streams ? a 
role determined by the flow field near LWD. Moreover, water management agencies are 
 3
interested in including the effects of LWD in aquatic habitat assessments used for water 
resource allocations.   
The flow structure around LWD is very complicated and varied, which has been 
revealed by several previous field studies (e.g. Mutz 2000; Beebe 2001; Daniels & 
Rhoads 2003). Moreover, there is a large difference in scales between hydraulic 
processes and the critical ecological processes. Because of this scale difference, the 
small-scale flow heterogeneity (important to biota) is poorly represented in the channel-
scale numerical models currently used for aquatic habitat analysis (e.g. Waddle 2001). To 
address this problem, a new conceptual model is proposed to integrate the subgrid-scale 
physics into a coarse-grid scale model (Fu and Hodges, 2005).  
 
1.2 CHALLENGE AND MOTIVATION OF THE NUMERICAL MODEL 
In addition to the scale difference, another challenge for solving the river flow 
with LWD is to build a numerical model to simulate the unsteady river flow. River flow 
is usually characterized as shallow-water flow. As shown in Figure 1.1, the flow around 
LWD is indeed shallow as most of the LWD are exposed out of the water surface. This 
kind of river system is usually modeled by 2D depth averaged shallow water equations 
(SWE) (Arega et al. 2008; Burguete et al. 2006; Gejadze and Monnier 2007). 2D depth 
averaged SWE can be obtained by integrating the Navier-Stokes equation over the water 
depth with the underlying hydrostatic assumption. One of the difficulties in solving SWE 
numerically is the combination of the ?stiffness? (i.e. the unsteady term) and 
?nonlinearity? (i.e. the nonlinear advection term) in the momentum equations. Numerous 
finite-difference numerical models for solving SWE have been developed over the past 
several decades. Among these developed numerical models, explicit and semi-implicit 
schemes are the most widely used. For example, in the atmospheric or weather research 
community, explicit methods such as Leapfrog and Lax-Wendroff-type are the most 
popular ones (Mendez-Nunez and Carroll 1993). Casulli (Casulli 1990) proposed a semi-
implicit scheme, on which a number of  variations are based. A common property for 
both the explicit and semi-implicit methods is that they explicitly discretize the advection 
term. Hence, the nonlinearity does not arise in an explicit or a semi-implicit sachem. 
However, their stability is restricted by the Courant-Friedrichs-Lewy (CFL) condition 
because of the explicitness. A fully implicit scheme is desirable because we can use a 
relative large time step in the simulation. A fully implicit solution for an unsteady 
nonlinear advection problem involves solving a system of nonlinear equation at each 
step. The conventional techniques for solving a system of nonlinear equations require 
either 1) iterations of root-finding procedure or/and 2) calculations of the Jacobian. These 
two processes make a fully implicit solution more complex and computationally 
expensive. These motivate this research to develop a new implicit scheme to solve 
unsteady nonlinear advection problems.  
 
 
 
 
 4
 
1.3 OBJECTIVE 
The primary objective of this research is:  
Develop a new implicit solution for unsteady nonlinear advection problems. The 
new numerical algorithm should have the advantage in both stability and efficiency while 
keeping the 2
nd
-order temporal accuracy as in most existing SWE models. In addition, 
this new method would provide a novel approach in decoupling a coupled equation 
system. 
 
1.4 APPROACH 
To achieve the objective of this research, several steps need to be taken to guide 
this study. The major steps are as follows: 
Literature review 
? Review current finite-difference models for the 2D depth averaged SWE. 
? Examine how the nonlinear advection term is treated in different temporal 
discretizations.   
? Review the merits and limitations of the existing models.  
This part of the dissertation work is presented in Chapter 2.  
 
Analytical development of the new numerical algorithm 
Illustrate the theoretical principle and concept of the new algorithm. This part of 
the dissertation work is presented in Chapter 3. 
 
Verify the new method using 1D test case 
? Test the new algorithm using a 1D PDE with time dependent nonlinear 
advection term.  
? Test the new algorithm using 1D ordinary differential equation (ODE).  
? Explore the characteristics of the method in accuracy, stability and 
efficiency.   
This part of the dissertation work can be found in Chapter 4. 
 
Verify the new method using a multi-variable and multi-dimensional problem 
? Theoretically develop the new algorithm in multi-variable and multi-
dimensional problems.  
? Demonstrate and analyze the decoupling process of the new method.  
? Conduct initial test cases of the numerical algorithm in 2D depth averaged 
SWE.  
This work can be found in Chapter 5 and Chapter 6.  
 
 5
1.5 BROAD APPLICATION 
The initial motivation and application for this research is simulating a river flow 
with LWD. The newly developed numerical algorithm for solving unsteady nonlinear 
advection problems is not limited only in simulating river flows. In most of the 
environmental flows such as atmospheric flow, ocean circulation, storm surge and etc., 
advection is the dominant mechanism and we are often interested in simulating these 
flows with a range of different time scales. Therefore, the difficulties arising from the 
combination of the unsteadiness and nonlinear advection exist in modeling many types of 
environmental flows. This new numerical method provides a new approach to address 
this difficulty in these and many other research areas. 
 6
Chapter 2 Review of Finite Difference Schemes for Solving the 
Unsteady Nonlinear Advection in the Shallow Water Equations 
 
Environmental flows (e.g. river and atmospheric flows) governed by the SWE are 
usually dominated by the advective mechanism over multiple time-scales (Vreugdenhil 
1994). The combination of  time dependency and nonlinear advection creates difficulties 
in the numerical solution of the shallow water equations (SWE) (Bourchtein and 
Bourchtein 2006). To address this difficulty, numerous numerical schemes have been 
developed over the past several decades. Iskandarani et al (2005) reviewed the finite 
element/finite volume methods in solving advection equations. More reviews about finite 
element methods can be found in the article by Thomee (2001).  The research herein 
focuses on developing a finite difference algorithm for solving unsteady nonlinear 
advection problems. To review the existing finite difference techniques, we first classify 
time-marching algorithms into three categories: explicit, semi-implicit and implicit 
methods. In the explicit and semi-implicit methods, the nonlinear advection terms are 
treated explicitly. Hence, the issue of solving a nonlinear matrix inversion does not arise 
in these two temporal discretizations. However, to solve a fully-implicit system, different 
computational linearization methods are required treat the nonlinearity, so that the 
inversion of a nonlinear matrix can be obtained. Therefore, we briefly review the explicit 
and semi-implicit models for solving the SWE and examine different linearization 
approaches in the implicit system.  Comparisons among different existing methods are 
summarized in the last section of this chapter.  
 
2.1 EXPLICIT METHOD 
Fully explicit methods are popular in atmospheric or weather research 
community. The most widely used explicit methods in atmospheric research are Leapfrog 
and Lax-wendroff-type schemes (Mendez-Nunez and Carroll 1993). The Leapfrog 
method is a one-step method, requiring two time levels of known values to compute an 
unknown level. Modifications of the Leapfrog methods are used in simulating ocean 
circulation (Cho and Yoon 1998; Fujima and Shigemura 2000). The Lax-Wendroff-type 
method is a two-step Predictor-Corrector method. One of the popular variations of the 
Lax-wendroff scheme is developed by MacCormack (1969). The MacCormack scheme is 
extended to many research areas such as aerospace engineering (Hirose et al. 1991; 
Kurbatskii and Mankbadi 2004) and hydrological science in simulating free surface flow 
(Fennema and Chaudhry 1990; Garcia and Kahawita 1986; Kazezy?lmaz-Alhana et al. 
2005; Keshari and Koo 2007; Li and Jackson 2007), groundwater flow (Keshari and Koo 
2007) and dam-break shock wave (Li and Jackson 2007). Another commonly used 
explicit Predictor-Corrector method is the Runge-Kutta method (Delis and Katsaounis 
2005; Zhou et al. 2007). The Leapfrog and the two-step Predictor-Corrector methods can 
maintain 2
nd
-order temporal accuracy. Although simpler explicit schemes such as 
 7
forward-time scheme is also used in solving SWE (Murillo et al. 2008), it can only 
maintain 1
st
-order temporal accuracy.  
Explicit schemes are straightforward and easy to implement. Nonlinear advective 
terms in partial differential equations (PDE) are readily-computed from known time-
levels using fully-explicit time-marching methods. However, the implementations of 
explicit methods are restricted by stability requirements. Runge-Kutta and MacCormack 
methods typically require an excessively small time step to satisfy the Courant-
Friedrichs-Lewy (CFL) condition(Mendez-Nunez and Carroll 1993); the Leapfrog 
approach must be applied with appropriate numerical strategies such as alternate time 
levels of the solution (Agoshkov et al. 1994; Peyret and Taylor 1983; Zhou 2002) or  
sophisticated mode-splitting (e.g.Blumberg and Mellor 1987). 
 
2.2 SEMI-IMPLICIT METHOD 
Semi-implicit methods have been developed and widely used in temporal 
discretized SWE. In most semi-implicit methods, the nonlinear advection terms are 
discretized explicitly (Bonaventura and Rosatti 2002; Bourchtein and Bourchtein 2007; 
Casulli 1990; Casulli and Cattani 1994; Casulli and Cheng 1992; Kar 2006; Spitaleri and 
Corinaldesi 1997) so that the issue of solving a nonlinear matrix inversion does not arise. 
For instance, Environmental Fluid Dynamics Computer code (EFDC) uses a three-level 
semi-implicit method, in which advection terms are explicitly discretized using upwind 
scheme.  As a result, the stability is constrained by the explicit advection terms (Hamrick 
1992). Some other models (Kar 2006) use the explicit Runge-Kutta method to calculate 
the advection term. It is obvious that the stability of this treatment is constrained by the 
CFL number.  
To overcome the stability constraint, different numerical treatments are 
introduced to calculate the linearized advection term. One of the most widely used 
techniques is the semi-Lagrangian method (Cullen 2001; Rosatti et al. 2005). Casulli 
(Casulli 1990; Casulli and Cattani 1994; Casulli and Cheng 1992) used this method in 
UnTRIM and the earlier version, TRIM. Some latest models including ELCIRC (Zhang 
et al. 2004) and ELCOM(Hodges 2000) followed this idea. Although Casulli?s approach 
has been successful and stable, its accuracy is relatively poor (Hodges, 2004).  
The main advantage of the semi-Lagrangian method is that the stability is not 
constrained by the CFL condition (Barros and Garcia 2007).  However, a set of trajectory 
equations need to be solved at each time step and the maximum allowable time step is 
restricted by the convergence criteria of the iteration of the trajectory equations 
(Bourchtein and Bourchtein 2007).  
 
2.3 IMPLICIT METHOD 
An obvious advantage of a fully-implicit method is that the time step is not 
restricted by the CFL number. Various fully-implicit finite difference schemes have been 
developed. Steppeler (2006) solved a fully implicit, SWE based meteorological model 
using the Fourier transformation method at each grid point and each time step. Some 
 8
researchers derived a fully-implicit scheme based on Alternating Direction Implicit (ADI) 
method (Szymkiewicz 1992; Wilders et al. 1988). Yuan and Wu solved implicit Navier-
Stokes equations using a staggered finite difference Crank-Nicolson scheme (Yuan and 
Wu 2004).  Burguete and Garcia-Navarro (Burguete and Garcia-Navarro 2004) 
implemented a first order upwind implicit scheme to simulate river hydraulics problem.   
Although the fully-implicit methods are more accurate and robust, they are less 
common, possibly due to computational complexity and expense (Turek 1996). 
Nonlinearity in a fully-implicit method results in a system of nonlinear algebraic 
equations to be solved at each time step. Existing strategies for solving implicit nonlinear 
equations include iterative and non-iterative methods (Moin 2001). In the following, we 
will review the conventional techniques for implicit nonlinear solutions.  
The Newton method and Picard method (Ferziger and Peric 1999; Lehmann and 
Ackerer 1998; Paniconi et al. 1991) are the most widely used two-level iterative 
techniques for implicit time-marching of nonlinear equations. For steady-state nonlinear 
problems, Newton and Picard iterative algorithms have proven quite successful (Paniconi 
and Putti 1994).  These methods may be thought of as successive linear solutions of 
=Ax b
%
, where A
%
 is an approximation of A that is iteratively refined to solve the 
nonlinear problem.  However, for time-evolving CFD problems, each time-marching step 
using the Newton or Picard method requires an outer iteration applied over an inner 
solution of a linear equation set.  When the inner problem also requires an iterative 
solution (as is common in CFD), the time march for this doubly-iterative approach is 
computationally expensive.  The principal differences between the Picard and Newton 
methods are that the former is easier to implement and requires fewer computations per 
outer iteration but has only 1st-order convergence, whereas the latter is more difficult to 
implement but provides 2nd-order convergence.  Thus, which method is appropriate 
depends upon the difficulty in implementing the Newton method compared to the slower 
convergence of the Picard method. A further implicit technique, local linearization, has 
not been widely used in CFD but does provide time-marching with 2nd-order accuracy 
(Lomax et al. 1999) without an outer iteration.  By providing a single linear 
approximation of the nonlinear problem, local linearization side-steps the convergence 
issue of the outer iteration for Newton/Picard techniques.  However, local linearization 
requires discretization of the Jacobian, which is often difficult to derive and implement 
for typical CFD applications.   
To better illustrate these existing implicit linearization methods, we use a simple 
scalar ODE with a quadratic nonlinearity as an example. The nonlinear ODE is written as 
 
d
f( ,t)
dt
?
= ?  (2.1) 
Using a finite-difference Crank-Nicolson discretization (the simplest 2
nd
-order implicit 
method), the above can be approximated as 
 
()(){}()
n1 n n1 n 3
t
ff t
2
++
?
?=?+ ? +?+?O   (2.2) 
 9
where superscripts indicate the discrete time step and 
( )
n1
f
+
?  implies a quadratic 
nonlinear relationship in 
n1+
?  (e.g. 
n1+
?
n1+
? or 
n1
x
+
?
n1+
? ). Equation  (2.2) is our example 
of an implicit nonlinear equation. In the following, we will demonstrate the principles of 
existing techniques to linearize Equation (2.2) computationally. 
 
2.3.1 Newton methods 
The Newton method is one kind of root-finding algorithm. It starts with an initial 
guess and iteratively estimates the root of the function. The Newton method can be 
derived from a Taylor series expansion or a geometric proof (Amat et al. 2003). One of 
the advantages of the Newton method is that it can be applied to equation systems with 
complex nonlinearities and keeps the quadratic convergence rate. In petroleum 
engineering and other research areas concerning flow through porous media, the Newton 
method is a standard approach because of the complex nonlinearities (Cao and Sun 2005; 
Dettmer and Peric 2007; Kwok and Tchelepi 2007). Using the Newton method, the 
linearized equation system of Equation (2.2) can be written as 
 
()()
( )
k
n1
k1 k
n1 n1
k
n1
g
g
+
+
++
+
? ?
?
? ?
?=??
???
??
??
??
 (2.3) 
where the outer superscript indicates the iteration number and  
 
() ()( ){}
n1 n1 n n n1
t
gff
2
++ +
?
?=???? ?+?  (2.4) 
In general, we can use the value at the previous time step as the value at the first iteration, 
which means 
 
( )
1
n1 n+
? =?  (2.5) 
Equation (2.3) is solved repeatedly until the difference between two successive iterations 
satisfies the pre-defined convergence criteria. In Equation (2.3), a first order derivative 
()
n1
g/
+
??? must be calculated and updated at each iteration. If ?is a vector and 
discretized in space, the resulting 
( )
n1
g/
+
???  is a matrix and called the Jacobian of the 
linearized equation system. However,  the Jacobian matrix may be computationally 
expensive and difficult to calculate analytically (Ferziger and Peric 1999; Niet et al. 
2007).   
A number of approaches have been proposed to modify Newton method. For 
example, it can be modified with a Chebyshev approximation to accelerate the 
convergence (Bagatur 2007). To simplify the Newton method, much research has been 
carried out to simplify the calculation of the Jacobian matrix. Niet et al. (2007) evaluated 
the Jacobian matrix using a partial-analytical and partial-numerical technique in an 
ocean-climate model. Li (1993) attempted to reduce the effort to compute the derivatives. 
 10
A Jacobian-free Newton-Krylov (JFNK) method has been developed recently to solve 
nonlinear equation systems (Brown and Saad 1990; Chan and Jackson 1984; Knoll and 
Keyes 2004; Mousseau et al. 2002; Reisner et al. 2005; Wubs et al. 2006). The key to 
JFNK solver is approximating the Jacobian-vector product iteratively instead of 
evaluating each element of the Jacobian matrix. Although all these modifications reduce 
the computation effort of Jacobian, inconvenient iteration still remains in the Newton 
method.  
 
2.3.2 Picard Iteration 
The Picard iteration is more straightforward than the Newton method. In the 
Picard iteration, the previous outer iteration value, 
( )
k
n1+
? , is substituted for 
n1+
?  on the 
RHS of Equation (2.2), resulting in  
 
() () ( )
{ }
k1 k
n1 n n n1
t
ff
2
+
++
?
? ?
?=?+?+?
? ?
? ?
 (2.6) 
As in the Newton iteration, Equation (2.6) is solved repeatedly until the difference 
between two successive iterations satisfies pre-defined convergence criteria. The same 
initial condition, Equation (2.5) is used to start the solution. Due to the simplicity, Picard 
iteration has been widely used in fully implicit nonlinear solver in Computational Fluid 
Dynamics (CFD)(Clement et al. 1994; Kwag 2000; Webster 2007). Different 
modifications of the Picard method have been reported in the literature. For example, 
Celia et al.(1987) proposed a modified Picard iteration that is a combination of Picard 
iteration and Newton iteration. However, this and other modified Picard methods remain 
linearly convergent (Lehmann and Ackerer 1998). A detailed comparison of the Picard 
and Newton methods is provided by Paniconi and Putti (1994).  In general, Newton 
methods require more computations per iteration, but converge more rapidly (quadratic) 
than Picard methods (linear).  However, for advection equations, the overall 
computational cost of Newton methods is typically greater than Picard methods because 
of the computational cost of the Jacobian (Ferziger and Peric 1999).   
 
2.3.2 Local linearization 
In contrast to the above-mentioned Newton and Picard methods, local 
linearization is a non-iterative method.  Let?s first review the derivation of the local 
linearization techniques, which are based on a Taylor-series expansion for 
n1
f
+
about 
n
t in 
Equation (2.2). The following derivation is based on Lomax et al.(1999), 
 
()
n
n
n1 n 2
ff
ff t t
t
+
????
??
=+?? +? +?
????
?? ?
??
??
O  (2.7) 
 11
where 
n1 n+
??=? ?? . For advection equations, f is generally not an explicit function 
of t, so ()
n
tf/t???  in Equation (2.7) is zero.  Substituting Equation  (2.7) into Equation 
(2.2)  provides 
 
()
n
n1 n n n 3
tf
fft
2
O
+
??
??????
?=?+ +?? + +?
??
??
??
??
??
 (2.8) 
which is equivalent to   
 
n
n
tf
1tf
2
??
????
? ??=???
??
??
??
??
 (2.9) 
Equation  (2.9) is a linear second-order discrete form of the original nonlinear ODE and 
can be solved by any number of standard techniques once 
()
n
f/? ??  is explicitly 
computed. However, as in the Newton method, the evaluation of the Jacobian ()
n
f/??? 
complicates the overall computation.  
 
2.4 SUMMARY 
SWE is widely used in simulating environmental flows. Mathematically, the 
stiffness and nonlinear advection increases the difficulty in the numerical solutions of 
SWE. Explicit methods can easily solve the nonlinearity but are restricted by the CFL 
stability criteria. Semi-implicit methods treat the advection term explicitly but need 
additional numerical treatment such as the semi-Lagrangian method to calculate 
advection terms for a higher CFL number. A fully implicit method is desirable but 
requires further computational linearization to solve a nonlinear system of equations. 
Existing implicit linearization techniques (including Newton, Picard and local 
linearization methods) require either 1) additional explicit derivative evaluations to 
provide an approximate linear problem, or 2) an outer iteration that converges an inner 
approximate linear problem.  These two characteristics make implicit systems complex 
and expensive to solve. To address the nonlinearity in a fully implicit system, in the next 
chapter we propose a new computational linearization method that can linearize the 
nonlinear advection term without either 1) the outer iteration or 2) computation of the 
Jacobian. 
 12
Chapter 3 Theoretical Development of the Time-Centered Split (TCS) 
Method 
 
In Chapter 2 we reviewed the current implicit linearization methods that require 
either an outer iteration or a computation of the Jacobian. In this chapter, we develop a 
new method that is similar to local linearization in that it allows non-iterative 
discretization of a temporally 2
nd
-order approximation of the nonlinear equation set.  
Instead of requiring the Jacobian, the new method splits the time-marching nonlinear 
problem into two sets of linear problems that are solved in succession.  The new method 
has both the computational efficiency of a non-iterative local linearization method and 
the implementation simplicity of a Picard iterative method. We call this the Time-
Centered-Split (TCS) method. 
The theory of the TCS method is first derived using a generic single variable 
nonlinear equation in this chapter. The TCS method can generate different discrete forms 
by introducing the time-centered split to different terms. This advantage is demonstrated 
by applying the TCS method to a 1D advection-diffusion equation (Burgers? equation). 
Four different TCS formats such as TCSF (split the flux term), TCSG (split the gradient 
term), TCSF-D (split the flux and diffusion terms) and TCSG-D (split the gradient and 
the diffusion terms) are derived and presented for the 1D Burgers? equation. One of the 
key advantages of the TCS method is that it provides non-iterative coupling between 
multiple variables. A 1D coupled advection diffusion equation and scalar transport 
equation are used as an example to illustrate this advantage. A summary of the principles 
of the TCS method and its advantages are provided in the last section. 
 
 
3.1 THE THEORY OF TIME-CENTERED SPLITTING  
3.1.1 Computational Splitting of the Nonlinear Term 
The TCS method is based on a nested application of the midpoint rule (i.e. a 
centered-time approximation).  Midpoint rule discretizations are often used for time-
marching to obtain second-order temporal accuracy (Ferziger and Peric 1999).  A 
common approach for nonlinear time-marching in meteorology and oceanography is 
application of the midpoint rule in an explicit formulation known as the 3-level Leapfrog 
method (Dubois et al. 2005; Fujima and Shigemura 2000).  The explicit Leapfrog method 
can be written as 
 
( )
n1 n1 n n n 3
f2tt
+?
?=?+????+?,O() (3.1) 
where superscripts represent time levels and t?  is the model time step.  Equation (3.1) 
can be seen as a generic form of a nonlinear equation. The function f is a linear operator 
of ?? and ?. Inside the function f, the product term ?? represents a generic quadratic 
nonlinear term. This quadratic nonlinear term in a flow and transport problem is of the 
 13
interest in this research. To develop the TCS method, we note that in the same vein as 
Equation (3.1), the midpoint rule can be written across only two time levels as  
 
( )
n1 n n12 n12 n12 3
ftt
++++
?=?+? ? ? ?+?
// /
,O() (3.2) 
Equation (3.2) has the desirable property that the time n+1/2 information is retained only 
within the computation of the n to n+1 time step so that the advance is not leapfrogging 
over alternating data.  By introducing a time-centered linear approximation of  
 
( ) ( )
n1/2 n n1 2
/2 tO
++
?=?+? +? (3.3) 
 in one part of the nonlinear product 
n 1/2 n 1/2++
??in Equation (3.2), the discrete equation 
becomes computationally linear in time (i.e. no products with the same time level): 
 
()
nn1
n1 n 2 n12 n12 3
fOt tt
2
+
+++
??
???+?
?=?+ +? ? ? ?+?
??
??
????
//
,O() (3.4) 
Reorganizing Equation (3.4) provides  
 
()()
n1 n n n12 n12 n1 n12 n12 3
tt
ff t
22
+++++
? ?
?=?+?? ? +?? ? +?
// //
,,O() (3.5) 
Equation (3.5) is computationally split into two steps and an intermediate variable 
*
?  is 
defined as,   
 
()
* n n n 1/2 n 1/2
t
f,
2
++
?
?=?+ ?? ?  (3.6) 
Subtracting Equation (3.6) from Equation (3.5) provides 
 
()
n1 n1 n12 n12 3
t
ft
2
++++
?
?=?+?? ? +?
*//
,O() (3.7) 
After introducing the above time-centered split, Equations (3.6) and (3.7) are both 
computationally linear. Furthermore, Equation (3.6) is similar to the implicit Euler 
approximation of 
n1/2+
?  written as: 
 
()
2
n12 n n12 n12 n12
tt
f
22
////
,O
++++
? ?
??
?=?+??? +
??
??
 (3.8) 
Using the Taylor expansion, 
n1/2+
? can be expanded as: 
 
2
n12 n
t
tt
22
/
O
+
??
??
?=?+?+
??
??
 (3.9) 
Substituting Equation (3.9) into the nonlinear term of Equation (3.8), provides 
 
2
n12 n n n12 n12
ttt
f
22
///
O,O
+++
?????
??
?? ??
?=?+?+ ?? +
??
?? ??
??
?? ??
????
 (3.10) 
Grouping the higher order terms,  
 14
 
()
2
n12 n n n12 n12
tt
f
22
///
,O
+++
? ?
??
?=?+??? +
??
??
 (3.11) 
Substituting Equation (3.6) into Equation (3.11), results in 
 
2
n12
t
2
+
?
??
?=?+
??
??
/*
O  (3.12) 
Substitution of Equation (3.12)  into Equations (3.6) and (3.7) provides a two-step 
method  
 
()
nn 3
t
ft
2
?
?=?+ ??? + ?
***
,O() (3.13) 
 
()
n1 n1 3
t
ft
2
++
?
?=?+??? +?
***
,O() (3.14) 
Equations (3.13) and (3.14) are a computationally linearized implicit equation system. 
The summation of Equations (3.13) and (3.14) is a 2
nd
-order computational equivalent of 
Equation (3.2). 
Remark: 
Although the derivation above provides a two-step direct linearization method, an 
outer iteration might also be combined for a big t?  value. This process can be illustrated 
as the following:  
1. obtain 
n1+
?  from Equations (3.13) and (3.14). 
2. calculate  
n1/2+
? using 
( )
n1 n
/2
+
?+? , where 
n1+
?  is calculated from step 1. 
3. calculate a new 
n1+
?  using Equations (3.13) and (3.14) with the updated 
n1/2+
?  in step 2. 
Repeat the procedure until the difference between the old and new 
n1+
? is 
acceptable.  
 
3.1.2 Computationally Splitting the Linear Term  
In the previous section, we introduced time-centered split only into the nonlinear 
term of ??. In addition to splitting the nonlinear product term, one can also use the same 
splitting idea in the linear term?, following 
() ()
nn1 nn1
n1 n 2 n12 2 3
ft tt(t)
22
/
O,O O
++
++
??
?????+? ?+?
?=?+ +? ? +? ?+?
??
????
??????
 (3.15) 
Equation (3.15) is guaranteed 2
nd
-order accurate in time because of Equation (3.3).   
Reorganizing Equation (3.15) and grouping the higher order terms, 
 15
 
{}
n1 n n n12 n1 n12 n n1 3
t
f()
2
//
,,, O
++++
?
?=?+?? ?? ?? +? (3.16) 
In Equation (3.16), not only the nonlinear term but also the linear term is split into two 
parts.  Therefore, a different intermediate variable 
*
? is defined as 
 
()
*n nn1/2n
t
f,
2
+
?
?=?+ ?? ?  (3.17) 
The second step of the splitting system then becomes, 
 
()
n1 * n1 n1/2 n1
t
f,
2
++++
?
?=?+?? ?  (3.18) 
Equations (3.17) and (3.18) are different from the split system in the previous section. As 
a result, the correspondence between 
*
? and 
n12/+
?  is proved differently. Instead of the 
implicit Euler equation, the explicit Euler approximation for 
n12/+
?  is introduced,    
 
()
2
n12 n n n n
tt
f
22
/
,O
+
??
??
?=?+??? +
??
??
 (3.19) 
Substituting the Taylor expansion of 
n12/+
?  into Equation (3.17) and grouping the higher 
order terms,   
 
()
2
*n nnn
tt
f, O
22
??
??
?=?+ ??? +
??
??
 (3.20) 
Substituting Equation (3.20) into Equation (3.19), we obtain 
 
2
n12
t
2
+
?
??
?=?+
??
??
/*
O  (3.21) 
Thus, we can use 
*
?  to replace 
n12/+
? in Equations (3.17) and (3.18). A different set of 
computationally linearized two-step equations can be written as: 
 
()
*n n*n
t
f,
2
?
?=?+ ???  (3.22) 
 
()
n1 * n1 * n1
t
f,
2
+++
?
?=?+???  (3.23) 
The key to the TCS method is the second-order time-splitting of quadratic 
nonlinear terms to two different time levels, i.e. 
n*
? ? and 
n1 *+
? ? . The result is discretely 
linear in any time-level of information.  A further application of the same splitting in the 
linear terms will give another set of discrete linearized equations. The above derivation 
for the single variable ? can be readily extended to a vector of variables, 
12 N
[ , ,... ]?? ? , or 
variables and linear operators, e.g. 
112 2N 2
[ ,L( ), ,L( )... ,L( )]? ?? ? ? ? .   However, as 
additional variables (or operators) are introduced, the discretization method has multiple 
implementations.  For example, even with a 1D advection diffusion equation, we can 
 16
obtain at least four different TCS discrete formats. In the next section, we will present 
these four different TCS formats by applying the TCS method to a 1D Burgers? equation. 
 
3.2 DERIVATION OF THE TCS METHOD IN A 1D ADVECTION-DIFFUSION EQUATION 
The Burgers? equation is the simplest advection-diffusion test case and can be 
viewed as a prototype of the Navier-Stokes equation or the SWE. The non-conservative 
Burgers? equation can be written as: 
 
2
2
uu u
u
txx
? ??
+=?
? ??
 (3.24) 
To develop the TCS method, we note that in the same vein as Equation (3.7), the 
midpoint rule can be written across only two time levels as  
 
()( )
n1 n n12 n12 2 n12 3
xx
uu u u u t t
// /
O
++++
=+? ? +?? ?+? (3.25) 
where 
x
? is the shorthand notation of the generic discretized spatial derivative of ?u?.  
 
3.2.1 The TCSF Discrete Format 
The nonlinear advection term in Equation (3.24) is constructed by the flux part ?u? 
and the gradient part, ? ux/???. The time-centered splitting maybe applied to either part 
of the advection term. We first introduce a time-centered linear approximation of 
( ) ( )
n1/2 n n1 2
uuu/2tO
++
=+ +? to the flux term in the nonlinear product in Equation 
(3.25). The discrete equation becomes computationally linear in time (i.e. no products 
with the same time level): 
 
()
n1 n
n1 n n12 2 n12 3
xx
uu
uu u u t t
2
//
O
+
+++
??
??+
=+? ? +?? ?+?
??
??
????
 (3.26) 
 
Multiplying the products in Equation (3.26) provides 
 
()()()()
n1 n n1 n1/2 n n1/2 2 n1/2 3
xxx
tt
uu uu uu u tOt
22
+++ + +
? ?
=+?? +?? +?? ?+? (3.27) 
A computational-splitting technique is used to obtain a numerically-solvable set of 
equations from the discrete form of Equation (3.27).  Defining a generic intermediate 
variable 
*
u  as 
 
()()
n n n12 2 n12
xx
tt
uu uu u
22
*//++
? ?
=+?? +??  (3.28) 
and subtracting Equation (3.28) from Equation (3.27) provides 
 
()()()
n1 n1 n12 2 n12 3
xx
tt
uu uu u t
22
*/ /
O
+++ +
? ?
=+?? +?? +? (3.29) 
 17
Using an implicit Euler approximation of 
n12
u
/+
and the same mathematical derivation as 
in Equation (3.12), we can prove 
 
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
 (3.30) 
Substitution of Equation (3.30) into Equations  (3.28) and (3.29) provides a two-step 
method, which is called the TCSF method: 
 
{}
nn2
xx
t
uu uu u
2
***
?
=+ ??+??  (3.31) 
 
{}()
n1 n1 2 3
xx
t
uu uu u t
2
***
O
++
?
=+ ??+?? +? (3.32) 
The summation of Equations  (3.31) and (3.32) is equivalent to the original Equation 
(3.25). The first step is an implicit solution for the variable u
*
 involving only 
n
u and with 
an implicit discretization of the diffusion term. 
n1
u
+
 in the second step can be explicitly 
calculated. 
 
3.2.2 The TCSG Discrete Format 
Splitting the gradient term instead of the flux term in Equation (3.25) results in a 
second basic form of TCS.  Thus, instead of Equation (3.26), we obtain  
 
n1 n
n1 n n1/2 2 n1/2 3
xx
uu
uu u u t(t)
2
O
+
++ +
??
??+
=+? ? +?? ?+?
??
??
????
 (3.33) 
For this second form, we define a slightly different intermediate variable as 
 
()()
* n n1/2 n 2 n1/2
xx
tt
uu u u u
22
++
? ?
=+? ? +??  (3.34) 
Similarly in Equation (3.30), we have 
 
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
 (3.35) 
Subtracting Equation (3.34) from (3.33) and substituting Equation (3.35) in the same 
manner as the transition from Equation (3.27) through (3.32), we obtain a second discrete 
split form that we will call TCSG 
 
{}
*n *n 2*
xx
t
uu uu u
2
?
=+ ??+??  (3.36) 
 
{}
n1 * * n1 2 *
xx
t
uu uu u
2
++
?
=+ ?? +??  (3.37) 
The summation of Equations (3.36) and (3.37) is equivalent to the original Equation 
(3.25). The first step is an implicit solution for the variable u
*
 involving
n
u  and with an 
 18
implicit discretization of the diffusion term. The second step is an implicit solution for 
the variable 
n1
u
+
 involving only u
*
and with an explicit discretization of the diffusion 
term. 
 
3.2.3 The TCSF-D Discrete Format 
The TCSF and TCSG forms are obtained by introducing the time-centered split 
into the advection term. Further approximation of the diffusion term using the same time 
splitting techniques will create different discretizations. For instance, in Equation (3.26), 
we can substitute the time splitting into both the flux and diffusion terms,  
 
()
n1 n n1 n
n1 n n12 2 3
xx
uu uu
uu u t t
22
/
O
++
++
??
?? ??
=+? ? +?? ?+?
??
?? ??
?? ????
 (3.38) 
Defining  u
*
 as, 
 
()()
nnn12 2n
xx
tt
uu uu u
22
*/+
? ?
=+?? +??  (3.39) 
the second step follows as, 
 
()()()
n1 n1 n12 2 n1 3
xx
tt
uu uu u t
22
*/
O
+++ +
? ?
=+?? +?? +? (3.40) 
Using an explicit Euler approximation of 
n12
u
/+
and the same mathematical derivation as 
in Equation (3.21), we can prove that 
 
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
 (3.41) 
Substituting Equation (3.30) into Equations (3.39) and (3.40), we obtained the TCSF-D 
format as 
 
{}
nn2n
xx
t
uu uu u
2
**
?
=+ ??+??  (3.42) 
 
{}()
n1 n1 2 n1 3
xx
t
uu uu u t
2
**
O
+++
?
=+ ??+?? +? (3.43) 
We call Equations  (3.42) and (3.43) the TCSF-D method because we apply the time-
centered splitting to both the flux term and the diffusion term. The summation of 
Equations  (3.42) and (3.43) is equivalent to the original Equation (3.25). The first step is 
an implicit solution for the variable u
*
 involving
n
u  and with an explicit discretization of 
the diffusion term. The second step is an implicit solution for the variable 
n1
u
+
 involving 
only u
*
and with an implicit discretization of the diffusion term. 
  
 19
3.2.4 The TCSG-D Discrete Format 
Similar to the development of the TCSF-D method, substitution of time-centered 
splitting into the diffusion term in the TCSG provides another set of two-step equations 
as 
 
{}
*n *n 2n
xx
t
uu uu u
2
?
=+ ??+??  (3.44) 
 
{}
n1 * * n1 2 n1
xx
t
uu uu u
2
+++
?
=+ ?? +??  (3.45) 
Equations (3.44) and (3.45) are called the TCSG-D discrete form. The summation of 
Equations  (3.44) and (3.45) is equivalent to the original Equation (3.25). u
*
in the first 
step can be explicitly calculated and the diffusion term is explicitly discretized. The 
second step is an implicit solution for the variable 
n1
u
+
 involving only u
*
and with an 
implicit discretization of the diffusion term. 
We have derived four discrete forms of the TCS method applied to a 1D Burgers? 
equation in the previous sections. More possible discretizations will be generated when 
we apply the TCS method to an equation system with multiple variables. We will discuss 
this characteristic in Chapter 5. 
 
3.3 TCS FOR COUPLED MOMENTUM AND SCALAR TRANSPORT 
A key advantage of the TCS method is that it provides non-iterative coupling 
between multiple variables.   The method is best understood through application to a 
simple model problem.  First, we will consider the 1D advection-diffusion equation for a 
scalar ? 
 
2
2
u
txx
??????
+=?
? ??
 (3.46) 
where u is the velocity and ? is a diffusion coefficient.  Equation (3.46) can be coupled to 
the 1D Burgers? equation for momentum 
 
2
2
uu u
u
txx
? ??
+=?
? ??
 (3.47) 
Applying TCSF to Equation (3.46) provides 
 
2*
nn
2
t
u
2xx
?
?
? ?? ?? ? ?
?=?? ??
? ?
??
??
 (3.48) 
 
2*
n1 n1
2
t
u
2xx
?
+? +
? ?? ?? ? ?
?=?? ??
? ?
??
??
 (3.49) 
The first step, Equation (3.48), is an implicit equation for 
*
?  involving only time ?n? 
values of u, whereas the second step, Equation (3.49), linearly couples solution of 
n1+
?  to 
 20
n1
u
+
.  To complete the algorithm, Equation (3.47) can be similarly discretized so that the 
coupled TCSF method for both scalar transport and momentum can be written as linear 
operators, 
 
() ()
2
nn
2
t
1u
2x x
?
??
????
? ??
+???=?
????
??
??
??
 (3.50) 
 
() ()
2
nn
2
t
1u uu
2x x
?
??
????
???
+??=
????
??
??
??
 (3.51) 
 
2*
n1
2
2*
n1
2
tt
1
2x 2 x
tu t u
u01 u
2x 2 x
?
+?
?
+ ?
????????? ? ??
??+?
??????
=
??
?? ? ?
++?
??
??????
 (3.52) 
where the parentheses indicate a spatial derivative operator on the term to the right of the 
brackets.  The first step of the TSCF method for N variables over Q grid cells results in N 
independent linear problems of order Q.  The second step of TCSF results in a single 
linear problem of order NQ.  This series of linear problems is a 2
nd
-order temporal 
equivalent of the original time-marching, coupled, nonlinear momentum-advection 
problem.  As the TCS uses an intermediate solution ( ,u )
? ?
?  followed by final solution 
n1 n1
(,u)
++
? , it resembles predictor-corrector schemes.  However, classic predictor-
corrector methods are formulated with an explicit predictor of the time n+1 values 
followed by an implicit corrector to the time n+1 values, which makes the classic 
predictor-corrector methods restricted by the CFL condition; furthermore, predictor-
corrector methods do not provide an avenue for a simple nonlinear solution (see 
discussions in Lomax et al, 1999 and Tannehill et al, 1997).  In contrast, the new TCS 
method provides an implicit linearized predictor of the time n+1/2 values that are used in 
a coupled implicit solution of the time n+1 values.   This new approach provides a simple 
method of linearly-coupling equations that are nonlinearly-coupled in the original 
problem.  The method requires only linear matrix solutions and does not require the outer 
iteration of Picard or Newton methods.  As compared to the functional Jacobians required 
for local linearization and Newton iteration, the TCS coefficient matrices are relatively 
easy to derive and have forms very similar to the original model problem.   
 
3.4 SUMMARY 
A new computational linearization method, the TCS, is derived and analyzed in 
this chapter. Without iteration and calculation of the Jacobian, the TCS method splits the 
quadratic nonlinear term into two steps so that each step is computationally linear. That 
is, for time marching from known state 
n
x  to unknown state 
n1+
x  we use a linear solution 
of 
n
()=Ax x b
% %
%  followed by 
n1
()
+
=Axx b
% %
% %
%  that is a second-order equivalent to 
 21
n1 n1
()
++
=Ax x b.  In addition to the linearization, the TCS method can generate different 
TCS discrete forms.  All different TCS formats are mathematically equivalent to the 
original implicit midpoint rule discretization. As a result, they share the same 2nd-order 
temporal accuracy. Furthermore, we also demonstrated that the TCS method can couple 
multiple variables without iterations. Characteristics of the TCS method such as 
accuracy, stability and efficiency will be explored in the next chapter using examples 
with analytical solutions. 
 22
Chapter 4 Implementation of the TCS Method in 1D Problems 
 
The theory of the TCS method is illustrated in Chapter 3. To investigate the 
characteristics of the new method, we apply the TCS method to three different test cases: 
a 1D conservative Burgers? equation, a 1D non-conservative Burgers? equation and a 1D 
nonlinear ordinary differential equation (ODE).  For each test case, the TCS algorithm is 
compared to the conventional implicit nonlinear solution methods (local linearization, 
Picard iteration and Newton iteration) applied to Crank-Nicolson discretization. The 
temporal accuracy of different TCS discretizations is verified by all three test cases. The 
practical stability of the TCS method is confirmed using the unsteady flow test case with 
an analytical solution in both conservative and non-conservative forms.  The method is 
shown to require computational effort similar to local linearization, but does not require 
discrete computation of a functional Jacobian for solution. 
 
4.1 APPLICATION OF THE TCS METHOD TO THE 1D CONSERVATIVE BURGERS? 
EQUATION 
4.1.1 Discrete formats of the 1D Conservative Burger?s equation using various 
methods 
Burgers? equation provides a useful model problem for comparing the TCS to 
other nonlinear solution methods. We will begin from the conservative form of Burgers? 
equation. The application of the TCS method to the non-conservative Burgers? equation 
will be discussed in the next section. The conservative form of the 1D Burgers? equation 
is: 
 
 
()
2
2
2
u1 u
u
t2x x
? ??
+=?
?
 (4.1) 
 
Time-Centered Split Method 
For an equation with one variable and a simple quadratic nonlinearity, TCSF and 
TCSG collapse into a single form. For Equation (4.1), the TCSF and TCSG methods are 
identical, and can be presented as linear operators 
 
() ()
n2
n
2
u
t1
1uu
22 x x
?
??
????
? ??
+??=
????
??
??
??
 (4.2) 
 
( )
2*
n1
2
u
ttu
1uu
4x 2x
?
+?
???
????
+=+
??
??
??
 (4.3) 
 23
For simplicity in the following discussion, we use the generic symbols A, x and b to 
represent a coefficient matrix, the left-hand-side (LHS) variable vector and the right-
hand-side (RHS) known vector in a matrix equation of the form of Ax = b for various 
solution methods.  The number of grid points in space is Q.  A central difference spatial 
discretization is applied to derivatives in Equation (4.2) and(4.3). It is useful to define a 
viscous scale,?, 
 
2
t
x
??
??
?
 (4.4) 
and the ?i? grid cell CFL number for the 
{ }
Ln,,n1? ?+time level for 
{ }
i 1,2,...Q=  as 
 
L
L i
i
ut
C
x
?
?
?
 (4.5) 
Similarly, it will be useful to also define a diffusion operator at any time level as 
 
( )
LLLL
ii1ii1
Du2uu
+?
?? ? +  (4.6) 
and a nonlinear adjective gradient operator 
 
LLL LL
ii1i1i1i1
GCu Cu
+ +??
?? (4.7) 
It follows that Eq, (4.2), the first step of TCSF can be written in the form Ax = b over Q 
grid points where the A matrix is tridiagonal such that 
 
n
2
nn
13
nn
Q2 Q
n
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
?? +? ?+
??
?
?? +?
??
A OOO  (4.8) 
 
 
T
**
12 Q
u ,u ,...u
? ?
=
? ?
x  (4.9) 
 
 24
 
n
n
0
10
n
2
n
3
n
Q1
n
Q1n
QQ1
C
uu
28
u
u
u
C
uu
28
?
?
+ ?
+
? ?
???
++
? ?
??
??? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
??
?? ?
+?
??
? ?
??
? ?
b
M
 (4.10) 
with 
nn
00Q1
C,u,C
?
+
and 
Q1
u
?
+
 implemented as Dirichlet boundary conditions. Neumann 
boundary conditions can also be readily invoked, but are not presented here for brevity.   
The second step of TCS for the 1D Burgers? equation can be as evaluated from 
another Ax = b problem using 
 
2
13
Q2 Q
Q1
C
10
8
CC
1
88
CC
1
88
C
01
8
?
??
??
?
?
?
??
+
??
?+
=
??
?+
?
??
A OOO  (4.11) 
 
T
n1 n1 n1
12 Q
u,u,..u
++ +
? ?
=
? ?
x  (4.12) 
 25
 
0
11
22
33
Q1 Q1
Q1
QQ
1C
uD
28
1
uD
2
1
uD
2
1
uD
2
C
1
uD
28
?
??
??
??
??
??
?
+??
? ?
++
? ?
? ?
? ?
+
? ?
? ?
? ?
+
? ?
=
? ?
? ?
? ?
+
? ?
? ?
? ?
+?
? ?
? ?
b
M
 (4.13) 
Thus, the TCSF method for the 1D conservative Burgers? equation requires solution of 
two tridiagonal linear problems at each time step. 
If we further split the diffusion term, we obtain the TCSF-D form for Equation 
(4.1) as, 
 
( )
n
2n
n
2
uttu
1uu
4x 2x
?
???
????
+=+
??
? ?
??
 (4.14) 
 
() ()
2
n1
2
u
t1
u
22 x x
?
+ ?
????
?
+??=
??
??
????
??
 (4.15) 
The first step is a matrix equation Ax=b with  
 
n
2
nn
13
nn
Q2 Q
n
Q1
C
10
8
CC
1
88
CC
1
88
C
01
8
?
?
??
+
??
?+
=
??
?+
?
??
A OOO  (4.16) 
 
 
 26
 
T
** *
12 Q
u,u,...u
? ?
=
? ?
x  (4.17) 
 
 
 
n
nn
0
11
nn
22
nn
33
nn
Q1 Q1
n
Q1nn
QQ
1C
uD
28
1
uD
2
1
uD
2
1
uD
2
C
1
uD
28
??
+
? ?
++
? ?
? ?
? ?
+
? ?
? ?
? ?
+
=? ?
? ?
? ?
? ?
+
? ?
? ?
? ?
+?
? ?
? ?
b
M
 (4.18) 
In the second step, we have 
 
*
2
**
13
**
Q2 Q
*
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
??
?? +? ?+
?
?? +?
??
A OOO  (4.19) 
 
 
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?
=
? ?
x  (4.20) 
 
 
 27
 
n
*n1
0
10
*
2
*
3
*
Q1
n
Q1*n1
QQ
C
uu
28
u
u
u
C
uu
28
+
?
+ +
+
? ?
???
++
? ?
??
??
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
??
?
? ?
+?
??
? ?
??
? ?
b
M
 (4.21) 
 
Similar to TCSF, the TCSF-D method for the 1D conservative Burgers? equation also 
requires solution of two tridiagonal linear problems at each time step. However, TCSF-D 
reverses the solution process in TCSF. The first step in TCSF-D is in the same structure 
as the second step in TCSF; the second step in TCSF-D is in the same structure as the 
first step in TCSF. 
To compare the above TCS methods to other temporally 2
nd
-order accurate 
implicit nonlinear solution methods, we apply Crank-Nicolson 2
nd
-order temporal and 
central difference spatial 2
nd
-order discretization to Equation (4.1), resulting in 
 
nn nn n n n
n1 n
i1 i1 i1 i1 i1 i i1
ii 2
n1 n1 n1 n1 n1 n1 n1
i1 i1 i1 i1 i1 i i1
2
tuu uu u 2u u
uu
24x x
uu uu u 2u u
4x x
+
++ ?? + ?
++ ++ + + +
++ ?? + ?
??? ?+
=? ??
?
??
?
?
+??
?
??
?
 (4.22) 
Conventional linearization methods such as Picard iteration, Newton iteration and local 
linearization methods are used to solve Equation (4.22). 
Picard Iteration 
The simplest approach to implement for a nonlinear equation such as Equation  
(4.22) is a lagged-coefficient iteration (Tannehill et al. 1997,pg 450) , which is a form of 
Picard iteration. A linear inner equation is formed by estimating the time n+1 flux as 
()
k
n1
i
u
+
, where the additional superscript ?k? is introduced as the outer iteration counter. 
Equation (4.22) can then be represented as the Picard method 
 28
 
()
()()()()
() ()()
nn nn n n n
k1
n1 n
i1 i1 i1 i1 i1 i i1
ii 2
kk1 kk1
n1 n1 n1 n1
i1 i1 i1 i1
k1 k1 k1
n1 n1 n1
i1 i i1
2
tuu uu u 2u u
uu
24x x
uu uu
4x
u2u u
x
+
+
++ ?? + ?
++
++ ++
??
+++
+?
=?
??? ?+
??
?
??
?
?
+
?
?
?+
?
??
?
?
?
?
 (4.23) 
which is solved for k={1,2,3...} until an appropriate convergence criterion is reached.  .  
Thus, each outer iteration requires an inner linear solution of the form Ax = b.  For 
comparison with other methods, it is useful to define an 
kk1 k+
= +Ax b c  form in which 
the b vector requires only a single computation in each time step whereas the c vector is 
recomputed in each outer iteration. For Equation (4.23), the result is 
 
 
k
3
kk
24
k
kk
Q2 Q
k
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
??
?? +? ?+
?
?? +?
??
A OOO  (4.24) 
  
 
()() ()
T
k1 k1 k1
k1 n1 n1 n1
12 Q
u , u ,... u
++ +
++ + +??
=
??
??
x  (4.25) 
 29
 
nn n
11 1
nn n
22 2
nn n
33 3
nn n
Q1 Q1 Q1
nn n
QQ Q
11
uGD
82
11
uGD
82
11
uGD
82
11
uGD
82
11
uGD
82
???
? ?
?+
? ?
? ?
? ?
?+
? ?
? ?
? ??+
=
? ?
? ?
? ?
? ?
?+
? ?
? ?
? ?
?+
? ?
b
M
 (4.26) 
 
 
k
0
k
k
Q1
1
C
28
0
0
0
1
C
28
+
?
? ?
+
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ??
?
? ?
? ?
c
M
 (4.27) 
For the Picard iteration, the first outer iterative solution is started with 
()
1
n1 n
ii
uu
+
=  .  The 
outer iteration is stopped when 
()()
{ }
k1 k
n1 n1
+
++
? <?uuL  where L is an appropriate 
linear norm acting on the vector u and ? is the desired convergence.  The above method 
can be expected to have no better than first order convergence for the outer iteration.  
Solution for each time step requires multiple tridiagonal solutions of the =+Ax b c  
problem where the A matrix and boundary conditions (possibly) in c are re-calculated at 
each outer iteration. 
 
Newton Iteration 
The Newton method is commonly used for accelerated convergence in iterative solution 
of nonlinear problems. Using classic Newton iteration method, the time-march from ?n? 
to ?n+1? of Equation (4.22) can be written as the root-finding problem for the function g 
as 
 30
 
()
n1 n1 n1 n1 n
ii1i i1 i i
n1 n1 n1 n1 n1 n1 n1
i1 i1 i1 i1 i1 i i1
2
nn nn n n n
i1 i1 i1 i1 i1 i i1
2
g u ,u ,u u u
tuu uu u 2u u
24x x
uu uu u 2u u
0
2x x
+++ +
?+
++ ++ + + +
++ ?? + ?
++ ?? + ?
=?
??? ?
+??
?
??
?
???+
+ ?? =
?
??
?
 (4.28) 
Applying Newton iteration(Ferziger and Peric 1999), Equation (4.28) can be written as 
 
()
k
k1
n1 ki
ji
n1
j
g
ug
u
+
+
+
??
?
? =?
??
?
??
 (4.29) 
where  
 
()( ) ( )
k1 k1 k
n1 n1 n1
jjj
uuu
++
+
?= ? (4.30) 
and
k
i
n1
j
g
u
+
??
?
??
?
??
 is the Jacobian matrix, evaluated using term by term discretizations of 
Equation  (4.28).  The result can be presented as a linearized equation system 
kk1 k+
=+Ax b c  with 
 
k
2
kk
13
k
kk
Q2 Q
k
Q1
C
10
42
CC
1
42 42
CC
1
42 24
C
01
42
?
?
? ??
+? ?
? ?
? ?
??? ?
?? +? ?
? ?
? ?
=
? ?
? ?
??
? ?
?? +? ?
? ?
? ?
?
? ?
? ?+?
? ?
? ?
A
OO O
 (4.31) 
 
 
()()()
T
k1 k1 k1
k1 n1 n1 n1
12 Q
u , u ,... u
++ +
++ + +??
=? ? ?
??
??
x  (4.32) 
 
 31
 
nn
n 11
1
nn
n 22
2
nn
Q1 Q1n
Q1
nn
QQn
Q
GD
u
82
GD
u
82
GD
u
82
GD
u
82
??
?
? ?
+?
? ?
? ?
? ?
+?
? ?
? ?
=?
? ?
? ?
? ?
+?
? ?
? ?
? ?
+?
? ?
? ?
b
M
 (4.33) 
 
kkk
k 110
1
kk
k 22
2
k
kk
Q1 Q1k
Q1
kkk
QQQ1k
Q
GDc
u
8242
GD
u
82
GD
u
82
GDc
u
82 42
??
?
+
???
+?++
??
+?
??
=?
+?
??
?
+?? +
??
c
M
 (4.34) 
In the present 1D single-variable example, the Jacobian is relatively easy to compute and 
thus it is easy to form the A matrix.  However, in multi-dimensional, multi-variable 
equation systems, the element-by-element calculation of the A matrix from g / u?? is 
generally difficult to derive and code.  As the A matrix must be recomputed for every 
outer iteration, the resulting outer iterations can be computationally expensive.   
 
Local Linearization 
Local linearization is not widely used for time-marching CFD problems, but is 
presented here because it has some similarities to the TCS method.  Lomax et al. (1999) 
provides an example of solving an ODE using Local linearization  which we extend to 
solving Equation (4.22). Local linearization approximates the time march to the ?n+1? 
time level using ?n? time level and a Taylor expansion.  The result is a linear equation set 
without an outer iteration. To apply Local Linearization, Equation  (4.22) can be written 
as, 
 
()()
n1 n n1 n1 n1 n1 n n1 n1 n1
i i i i1i i1 i i1i i1
t
u u F u ,u ,u F u ,u ,u
2
++++++
?? ??
?
??
=+ +
??
 (4.35) 
where   
 32
 
()
LL LL L L L
LL LL i1 i1 i1 i1 i1 i i1
ii1ii1
2
uu uu u 2u u
Fu,u,u
4x x
+ +?? + ?
?+
??+
=? +?
? ?
 (4.36) 
for L {n,n 1}?+.  Following Lomax et al. (1999), apply a Taylor expansion onto 
n1
i
F
+
 , 
 
(){
()
()
()} ()
n1 n n n1 n1 n1
ii ii1ii1
n n1n1n1i1
ii1i i1
n1 n
jj n
j
ji1
n n1n1n1 3
ii1i i1
t
u u F u ,u ,u
2
F u ,u ,u
uu
u
F u ,u ,u O t
+++
??
++++
??
+
=?
+++
??
?
=+
?
+?
?
++?
?
 (4.37) 
Equation (4.37) can be reorganized as 
 
n
n1 ni
ji
n
j
tF
1utF
2u
+
??
??
? ?=?
??
?
??
 (4.38) 
where 
nn
ij
F/u?? is the Jacobian and 
 
n1 n1 n
jjj
uuu
++
? =? (4.39) 
Therefore, the Local Linearization system can be written as 
n1+
=Ax b : 
 
 
n
2
nn
13
nn
Q2 Q
n
Q1
C
10
42
CC
1
42 42
CC
1
42 24
C
01
42
?
?
? ??
+? ?
? ?
? ?
??
? ?
?? +? ?
? ?
? ?
=
? ?
? ?
??
? ??? +? ?
? ?
? ?
?
? ?? ?+?
? ?
A
OO O
 (4.40) 
 
T
n1 n1 n1 n1
12 Q
u,u,..u
++++
? ?
=? ? ?
? ?
x  (4.41) 
 33
 
nn
n
10
1
n
n
2
2
n
Q1 n
Q1
nn
QQ1n
Q
GC
D
442
G
D
4
G
D
4
GC
D
442
?
?
+
? ??
?+++
? ?
? ?
? ?
?+
? ?
? ?
=
? ?
? ?
? ??+
? ?
? ?
?
? ??+? +
? ?
b
M
 (4.42) 
Much like the TCS method, Local Linearization provides a linear equation system that is 
a second-order approximation of the original nonlinear system and does not require an 
outer iteration.  In the above 1D single variable example, the Jacobian is relatively easy 
to compute and derive the A matrix, however, in multi-dimensional, multi-variable 
systems the Jacobian may be difficult to derive and compute.  Indeed, the difficulty of 
deriving the Jacobian for general multi-dimensional systems is arguably the principal 
reason that local linearization has not been widely used in CFD.  
 
4.1.2 Results and Discussion  
Accuracy 
The accuracy of the TCS method is demonstrated and compared to the other 
nonlinear methods for the unsteady 1D Burgers? equation.  Over the domain 0x1?? 
with initial and boundary conditions of 
 ( )u0,t 0=  (4.43) 
 
( )u1,t 0=  (4.44) 
 u(x,0) sin( x) : 0 x 1= ?<< (4.45) 
the 1D Burgers? equation has a solution constructed from truncated Fourier series:  
 
()
K
22
m
k1
analytical K
22
0k
k1
a exp( k t)ksin(k x)
u (x,t) 2 k 1,2,3,...K
a a exp( k t)cos(k x)
=
=
??? ?
=?? =
+????
?
?
 (4.46) 
where 
 () ()
{}
1
1
0
0
a exp 2 1 cos x dx
?
=???????
??
?
 (4.47) 
 34
 
() ()
{}
()
1
1
k
0
a 2 exp 2 1 cos x cos k x dx
?
=????? ???
??
?
 (4.48) 
The solution becomes exact as K ??.  This specific case of Burgers? equation has been 
frequently used to verify numerical methods (Kadalbajoo and Awasthi 2006; Kutluay et 
al. 1999).  A discrete approximation of Equation  (4.46) for a truncated series of K=30 
terms is plotted in Figure 4.1. The integrals in Equations (4.47) and (4.48) are 
numerically evaluated using ?quad? function in MATLAB with x? =0.01 and 0.05?= .   
 
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
u
t = 0
t = 0.6
t = 3
t = 0.3
t = 1
 
Figure 4.1 Solution of 1D Burgers? equation evolving in time for t?{0, 0.1, 0.2, 0.3, 0.6, 
1, 2, 3} using 0.05?=  along with the initial and boundary conditions of 
Equations   (4.43) through (4.45). 
 
We have solved this 1D Burgers? equation for the various nonlinear methods 
discretized in the previous session.  All methods use 2
nd
-order central difference spatial 
discretization and a sufficiently fine uniform mesh ( x0.02? = ) such that error is 
controlled by temporal accuracy.  As a basis for error comparison, we use the time of 
maximum absolute error found in these two TCS methods, which is determined from 
Figure 4.2 as t = 0.3 for both the TCSF and TCSF-D. Note that the maximum error 
 35
occurs when the solution reaches its maximum steepness in Figure 4.1. With a smaller ?, 
the similar asymptotic behavior is obtained when t? is sufficiently small.  
  
 
0 2 4 6
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
?4
t
?
ab
s
 
 
TCSF
TCSF?D
 
Figure 4.2 Absolute error time evolution for numerical solutions of 1D conservative 
Burgers? equation for TCSF and TCSF-D using { t? =0.01, x? =1/50, 
?=0.05}, where 
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
 and 
analytical
u (x, t)  is numerically calculated from Equation  (4.46) with K=30. 
 
 
To evaluate the model error, we examine performance with nine different time 
steps over the range 0.006 t 0.15???  and with grid cell spacing of 
{ }
23 3
x 2 10 5 10 1 25 10
?? ?
?? ? ? ?,,. . Based on the Burgers? equation solution for 
analytical
u, 
the test conditions cover the range 0.25 CFL 50? ? , where  
 36
 
analytical
t
CFL max u (t)
x
?
=
?
 (4.49) 
is the maximum CFL number. The error for the different methods has been estimated 
using the RMS error, 
2
L  and L
?
over Q grid points.  These approaches are similar to 
Rueda and Schladow (2002), where 
 
()
Q
2
RMS j j
j1
1
uu
Q
=
?= ?
?
%  (4.50) 
 
()
()
1/2
M
2
kk
k1
2 1/2
M
2
k
k1
uu
L
u
=
=
? ?
?
? ?
? ?
? ?
=
??
??
??
?
?
%
%
 (4.51) 
 
kk
k
max u u
L:1kM
max u
?
?
=??
%
%
 (4.52) 
and u%  is the model solution using an extremely fine time step such that 
4
0.3/10 , which is 
0.5% of the smallest  time step used in the other simulations.  Using a fine-time step 
solution is a standard approach (e.g. Ferziger and Peric, 1999) and is preferred over the 
approximate numerical solution of Equation (4.46), so that the accuracy measure is not 
distorted by the differences in the numerical approximations. Similarly, for consistency in 
the inter-model comparison, each model is compared to the small time-step solution for 
that model (e.g., the Picard method is compared to the small time-step solution of the 
Picard method, not to a small time-step TCSF solution).   Figure 4.3 shows the RMS 
error associated with different methods over the tested ranges of CFL.  It can be seen that 
all the methods provide 2nd-order accuracy for fixed x? and decreasing t? . Error 
magnitudes are substantially similar using the same x? and t?  except for the TCSF-D 
method. The TCSF-D method has a noticeably smaller error than all the other methods.  
2
L and L
?
 norms behave similar to RMS error and are plotted in Figure 4.4 and Figure 
4.5 respectively.  
 37
10
0
10
1
10
2
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
?
RMS
 
 
TCSF
TCSF?D
Picard
Newton
Local linearization
2nd order slope
?x=1/50
?x=1/200
?x=1/800
  
Figure 4.3 
RMS
? vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and 
t0.3/?= ?with ??{50, 30, 25, 20, 10, 8, 5, 4, 2}.  
 
 38
10
0
10
1
10
2
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
L
2
 
 
TCSF
TCSF?D
Picard
Newton
Local linearization
2nd order slope
?x=1/200
?x=1/800
?x=1/50
 
Figure 4.4 
2
L  norm vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and 
t0.3/?= ?with ??{50, 30, 25, 20, 10, 8, 5, 4, 2}. 
 
 
 39
10
0
10
1
10
2
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
L
?
 
 
TCSF
TCSF?D
Picard
Newton
Local linearization
2nd order slope
?x=1/800
?x=1/50
?x=1/200
 
Figure 4.5 L
?
norm vs. CFL number of various methods for solution of 1D conservative 
Burgers? equation with three different x? , where 0.05?=  and 
t0.3/?= ?with ??{50, 30, 25, 20, 10, 8, 5, 4, 2}. 
 
Stability 
Nonlinear stability cannot be generally proven, but is instead demonstrated by 
examples.  As shown in Figure 4.3, all five numerical methods remain stable over the full 
range of tested conditions. The stability of Picard iteration, Newton iteration and Local 
linearization may be expected because all the underlying discretizations begin with the 
Crank-Nicolson scheme. The TCS method is based on a midpoint rule discretization that 
is split, which cannot be guranteed stable even if the split matrices are linearly stable. Our 
tests thus far have shown very promising stability characteristics; however, such result 
cannot be presumed definitive. The new TCS method may theoretically be applied for 
any class of quadratically-nonlinear coupled equations, but its stability characteristics for 
different equations remains a subject for future investigation.   
 40
The TCS method is in an implicit format, thus it is expected to have stability 
advantage over the explicit methods. To compare the TCS method with explicit methods, 
we solve the same 1D conservative Burgers? equation using Runge-Kutta 2
nd
-Order 
(RK2) and Runge-Kutta 4
th
-Order (RK4) methods. The RMS error over a range of CFL 
numbers is plotted in Figure 4.3 for the TCSF, TCSF-D, RK2 and RK4 methods. In this 
comparison, we study the performance with 12 different time steps over the range 
0 0006 t 0 06..?? ?  and with grid cell spacing of x002.? = . As a result, the test 
conditions cover the range, 0 0013 CFL 2 57..? ? . The points of th RK2 and RK4 
methods are only shown in the left part of Figure 4.6 because these two explicit methods 
become unstable when CFL?0.3. Contrarily, the TCSF and TCSF-D methods remain 
stable in the whole range of CFL numbers. 
 
 
10
?2
10
?1
10
0
10
1
10
?15
10
?14
10
?13
10
?12
10
?11
10
?10
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
?
RMS
 
 
TCSF
TCSF?D
RK2
RK4
2nd order slope
unstable
unstable
 
Figure 4.6 
RMS
?  vs. CFL number of TCSF, RK2 and RK4 methods for solution of 1D 
conservative Burgers? equation, where 0.05?= x? =1/50, and 
t0.3/?= ?with ??{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 5}. 
 
 41
Computational Requirements 
To gain a better understanding of the computational requirements of the different 
methods, it is useful to compare the ideal operation counts using tridiagonal solutions of 
the Ax = b problems.  We focus on idealized operation counts rather than CPU time as 
we solve the above discrete problems using a MATLAB script, for which the software 
overhead tends to dominate the model runtime.   Based on algorithms in Press et al 
(1992) and Pozrikidis et al (2005, pg50), a tridiagonal solution ideally requires 8Q 
operations and a pentadiagonal solution requires 24Q operations, where Q is the number 
of grid points.  Here we do not distinguish between addition, multiplication or division 
operations.  If the Newton method requires R outer iterations for convergence, the Picard 
method should require R
2
 outer iterations for the same convergence level.  A tridiagonal 
Matrix algorithm (TDMA) is used for each method. A comparison of operations per grid 
point for the different methods is shown in Figure 4.7.   Although these operation counts 
are idealized and only for the 1D Burgers? equation with direct inner solutions, they 
represent the general trends that should be expected in any comparison of nonlinear 
solution methods.  The operation counts for both the TCSF and TCSF-D methods are 
identical because the solution processes for both methods are the same although they are 
in a reversed manner. For this 1D test case with an analytical Jacobian, the TCS method 
does not have any computational efficiency advantage over Local Linearization.  
However, the TCS derivation does not require a functional Jacobian computation, so 
extending the method to multiple dimensions and multiple variables is arguably easier, an 
idea that will be explored in next chapter.  
 
 
 
1 2 3 4 5
0
50
100
150
200
250
300
350
R
operations per grid point
 
 
TCSF
TCSF?D
local linearization
Newton iteration
Picard iteration
 
 42
Figure 4.7 Ideal operations per grid point for one time step using various nonlinear 
solution methods for 1D conservative Burgers? equation. R is the number of 
outer iterations taken by Newton method.  It is assumed that the Picard 
method convergences in R
2
 outer iterations. 
 
4.2 APPLICATIONS OF THE TCS METHOD TO THE 1D NON-CONSERVATIVE 
BURGERS? EQUATION 
Unlike in the conservative Burgers? equation, the TCSF and TCSG formats 
remain distinct for the non-conservative Burgers? equation. Applying the TCS method to 
the non-conservative Burgers? equation can create four different discretizations: the 
TCSF, TCSG, TCSF-D and the TCSG-D. In this section, we analyze all four different 
discretizations in the same vein as in the conservative case.   
 
4.2.1 Discrete formats of the 1D Non-Conservative Burger?s equation using various 
methods 
The TCSF Discrete Form 
The non-conservative Burgers? equation can be written as:  
 
2
2
uu u
u
txx
? ??
+=?
? ??
 (4.53) 
Introducing the approximation in the flux term results in the TCSF form as below: 
 
*2*
*n n
2
tuu
uu u
2xx
? ????
=+ ? +?
? ?
??
??
 (4.54) 
 
*2*
n1 n1
2
uut
uu u
2xx
+? +
? ????
=+ ? +?
? ?
??
??
 (4.55) 
Similar to the analysis in the conservative Burger?s equation, in the following discussion, 
we examine each method in a generic matrix equation, Ax=b. A central difference 
discretization is used for all the spatial derivatives. The same parameters such as the 
number of spatial grid points Q, the viscous scale?, the grid cell CFL number
L
i
C, the 
diffusion operator 
L
i
D , and the nonlinear adjective gradient operator 
L
i
G are also used in 
the following discussion. In the first step of the TCSF,  
 43
 
n
2
nn
13
nn
Q2 Q
Q1
C
1000
24
CC
10
24 24
00
CC
00 1
24 24
c
00 0 1
24
?
?
???
+? ? +
??
??
?? +? ?+
=
??
?? +? ?+
??
?
?? +?
??
A OO O  (4.56) 
 
 
T
** *
12 Q
u ,u ,...u? ?=
? ?
x  (4.57) 
 
 
n
n*
0
10
n
2
n
3
n
Q1
n
Q1n*
QQ1
C
uu
24
u
u
u
C
uu
24
?
+
+
? ?
???
++
? ???
??
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
??
?
+?
? ???
? ?
??
? ?
b
M
 (4.58) 
The second step of the TCSF is in an explicit format, which can be represented as  
 
()
()
***
ii1ii1
n1
i
**
i1 i1
uu2uu
2
u ; i {1,2,...Q}
1
1cc
4
?
+?
+
+?
?
+?+
==
+?
 (4.59) 
Thus, the TCSF solution for the 1D non-conservative Burgers? equation only needs one 
tridagonal linear problem at the first step and an explicit solution for 
n1
u
+
 at the second 
step.  
  
The TCSG Discrete Form 
If the linear approximation is introduced in the gradient rather than the flux part, the 
resulting discrete equations are the TCSG form:  
 
n2*
*n *
2
uut
uu u
2xx
? ????
=+ ? +?
? ?
??
??
 (4.60) 
 44
 
n1 2 *
n1 * *
2
uut
uu u
2xx
+
+
? ????
=+ ? +?
? ?
??
??
 (4.61) 
The first step of the TCSG has  
 
()
()
()
()
nn
31
nn
42
nn
QQ2
nn
Q1 Q1
CC
10
42
CC
1
242
CC
1
24 2
CC
01
24
?
+?
? ?
?
?
++??? ?
? ?
? ?
?
??
? ?
?++??
? ?
? ?
=
? ?
?? ?
??
?+ +? ?
? ?
? ?
? ?
?
?
? ?? ++?
? ?
A OO O (4.62) 
 
 
T
** *
12 Q
u ,u ,...u? ?=
? ?
x  (4.63) 
and 
 
n*
11
n
2
n
3
n
Q1
n*
QQ1
uu
2
u
u
u
uu
2
?
+
?
? ?
+
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
?
? ?+
? ?
b
M
 (4.64) 
 
 
 
 
 
 
 
 
 
 
 
 
In the second step of the TCSG,   
 45
 
*
1
**
22
**
Q1 Q1
**
QQ
C
10
4
CC
1
44
CC
1
44
CC
01
44
??
??
??
?
=
??
?
?
??
A
OO O
 (4.65) 
 
 
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?
=
? ?
x  (4.66) 
and 
 
**
*n1
11
11
*
*
2
2
*
Q1*
Q1
**
QQ*n1
QQ
DC
uu
24
D
u
2
D
u
2
DC
uu
24
+
?
?
+
? ?
++
? ?
? ?
? ?
+
? ?
? ?
=
? ?
? ?
? ?+
? ?
? ?
? ?+?
? ?
b
M
 (4.67) 
Thus the TCSG method for the 1D non-conservative Burgers? equation requires solution 
of two tridiagonal linear problems at each time step. 
 
 
 
The TCSF-D Discrete Form 
If we further split the diffusion term in the TCSF, we obtain the TCSF-D form. The 
TCSF-D form for Equation (4.53) can be written as, 
 
*2n
*n n
2
tuu
uu u
2xx
? ????
=+ ? +?
? ?
??
??
 (4.68) 
 
*2n1
n1 n1
2
tuu
uu u
2xx
+
+? +
? ????
=+ ? +?
? ?
??
??
 (4.69) 
In the first step, 
 46
 
n
1
nn
22
nn
Q1 Q1
nn
QQ
C
10
4
CC
1
44
CC
1
44
CC
01
44
??
??
??
?
=
??
?
?
??
A
OO O
 (4.70) 
 
T
** *
12 Q
u ,u ,...u? ?=
? ?
x  (4.71) 
and 
 
nn
n*
11
11
n
n
2
2
n
Q1n
Q1
nn
QQn*
QQ
DC
uu
24
D
u
2
D
u
2
DC
uu
24
?
?
? ?
++
? ?
? ?
? ?
+
? ?
? ?
=
? ?
? ?
? ?+
? ?
? ?
? ?+?
? ?
b
M
 (4.72) 
 
In the second step, 
()
()
()
()
**
31
**
42
**
QQ2
**
Q1 Q1
CC
10
42
CC
1
242
CC
1
24 2
CC
01
24
?
+?
? ?
?
?
++??? ?
? ?
? ?
?
??
? ?
?++??
? ?
? ?
=
? ?
?? ?
??
?+ +? ?
? ?
? ?
? ?
?
?
? ?? ++?
? ?
A OO O (4.73) 
 
 
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?=
? ?
x  (4.74) 
and 
 47
 
*n1
11
*
2
*
3
*
Q1
*n1
QQ
uu
2
u
u
u
uu
2
+
?
+
+
?
? ?
+
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
?
? ?+
? ?
b
M
 (4.75) 
It can be observed that the TCSF-D method reverses the solution process of the TCSG 
method. The first step in TCSF-D is in the same structure as the second step in TCSG; the 
second step in TCSF-D is in the same structure as the first step in TCSG. 
 
 
The TCSG-D Discrete Form 
The TCSG-D is obtained by further splitting the diffusion term in the TCSG 
method written as 
 
n2n
*n *
2
tuu
uu u
2xx
? ????
=+ ? +?
? ?
??
??
 (4.76) 
 
n1 2 n1
n1 * *
2
tu u
uu u
2xx
++
+
? ??? ?
=+ ? +?
? ?
??
??
 (4.77) 
u
*
 in the first step of Equation (4.76) can be explicitly calculated, 
 
()
()
nnnn
ii1ii1
*
i
nn
i1 i1
uu2uu
2
u ; i {1,2,...Q}
1
1cc
4
+?
+?
?
+?+
==
+?
 (4.78) 
The second step is a matrix equation with 
 48
 
*
2
**
13
**
Q2 Q
*
Q1
C
1000
24
CC
10
24 24
00
CC
00 1
24 24
C
00 0 1
24
?
?
???
+? ? +
??
??
?? +? ?+
=
??
?? +? ?+
??
?
? ?+?
??
A
OO O
 (4.79) 
 
 
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?
=
? ?
x  (4.80) 
 
 
*
*n1
0
10
*
2
*
3
*
Q1
*
Q1*n1
QQ
C
uu
24
u
u
u
C
uu
24
+
?
+ +
+
? ?
???
++
? ???
??
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
??
?
+?
? ???
? ?
??
? ?
b
M
 (4.81) 
The TCSG-D method reverses the solution process of the TCSF. The first step in TCSG-
D is in the same structure as in the second step in TCSF; the second step in TCSG-D is in 
the same structure as in the first step in TCSF. 
Similar to the last section, we first apply CN 2
nd
-order discretization to Equation 
(4.53), such that Equation  (4.53) is 
 
nn n nn
n1 n n
i1 i1 i1 i i1
ii i 2
n1 n1 n1 n1 n1
n1
i1 i1 i1 i i1
i 2
tuu u2uu
uu u
22x x
uu u2uu
u
2x x
+
+? + ?
++ + ++
+
+? ?
??? ?+
=? ??
?
??
?
???
+??
?
??
?
 (4.82) 
Equation (4.82) is then linearized using conventional Picard, Newton and local 
linearization methods.  
 
 
 49
 
 
 
Picard iteration 
Applying Picard iteration to Equation (4.82), the linearized equation can be written as, 
 
()
()
()()
() ()()
nn n nn
k1
n1 n n
i1 i1 i1 i i1
iii 2
k1 k1
n1 n1
k
i1 i1
n1
i
k1 k1 k1
n1 n1 n1
i1 i i1
2
tuu u2uu
uuu
22x x
uu
u
2x
u2u u
x
+
+
+? + ?
++
+?
+
+++
+?
=?
??? ?+
??
?
??
?
?
+
?
?
?+
?
??
?
?
?
?
 (4.83) 
where the superscripts n and k indicate the time step and iteration number respectively. 
with 
 
k
3
kk
24
k
kk
Q2 Q
k
M
C
10
24
CC
1
24 24
CC
1
24 24
C
01
24
?
???
+? ? +
??
??
?? +? ?+
=
??
?? +? ?+
??
?
?? +?
??
A OOO  (4.84) 
 
 
()() ()
T
k1 k1 k1
k1 n1 n1 n1
12 Q
u,u,.u
++ +
++ + +
??
=
??
??
x  (4.85) 
 
 
 50
 
nn k
n
11 1
1
nn
n
22
2
nn
Q1 Q1n
Q1
nn k
QQ Qn
Q
DG C
u
2424
DG
u
24
DG
u
24
DG C
u
2424
??
?
? ??
??++
? ?
? ?
? ?
??
? ?
? ?
=
? ?
? ?
? ???
? ?
? ?
?
? ???+?
? ?
b
M
 (4.86) 
 
 
k
1
k
k
Q1
C
24
0
0
C
24
+
? ??
+
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
?
? ?
?
? ?
? ?
c M  (4.87) 
 
 
Newton Iteration 
Applying the Newton iteration to Equation (4.82), the linearized equation can be written 
as: 
 
()
k
k1
n1 ki
ji
n1
j
g
ug
u
+
+
+
??
?
? =?
??
?
??
 (4.88) 
with 
 
()
n1 n1 n1 n1 n
ii1i i1 i i
n1 n1 n1 n1 n1
n1
i1 i1 i1 i i1
i 2
nn n nn
n
i1 i1 i1 i i1
i 2
g u ,u ,u u u
tuu u2uu
u
22x x
uu u2uu
u0
2x x
+++ +
?+
++ + ++
+
+? ?
+? + ?
=?
??? ?
+??
?
??
?
???+
+ ?? =
?
??
?
 (4.89) 
and 
 
()( ) ( )
kk1k
n1 n1 n1
uu u
+
++ +
?= ?  (4.90) 
Equation (4.88) is a matrix equation Ax=b with  
 51
 
kk
23
kk k
23 4
k
kk k
Q2 Q1 Q
kk
Q1 Q
GC
10
424
CG C
1
24 4 24
CG C
1
24 4 24
CG
01
24 4
??
?
???
++??+
??
??
?? + +? ?+
=
??
?? + +? ?+
??
?
? ?+ +?
??
A
OO O
 (4.91) 
 
 
()()()
T
k1 k1 k1
k1 n1 n1 n1
12 Q
u , u ,... u
++ +
++ + +
??
=? ? ?
??
??
x  (4.92) 
 
nn
n
11
1
nn
n
22
2
nn
Q1 Q1n
Q1
nn
QQn
Q
GD
u
42
GD
u
42
GD
u
42
GD
u
42
??
?
? ?
+?
? ?
? ?
? ?
+?
? ?
? ?
=?
? ?
? ?
? ?+?
? ?
? ?
? ?+?
? ?
b
M
 (4.93) 
 
kk k
k
11 0
1
kk
k
22
2
k
kk
Q-1 Q-1k
Q-1
kk k
QQ Q1k
Q
GD C
u+
4224
GD
u+
42
GD
u+
42
GD C
u+
4224
+
? ??
?++
? ?
? ?
? ?
?
? ?
? ?
=
? ?
? ?
? ??
? ?
? ?
?
? ??+?
? ?
C
M
 (4.94) 
 
Local linearization 
Applying Local linearization to Equation (4.82), the linearized equation system can be 
written as: 
 52
 
n
nn
i
ji
j
1F
1t utF
2u
??
??
?
??
? ??=?
??
?
??
??
 (4.95) 
where  
 
()
LL L LL
LL LL L
i1 i1 i1 i i1
ii1ii1 i 2
uu u2uu
F u ,u ,u u
2x x
+ ?+ ?
?+
??+
=? +?
??
 (4.96) 
for L{n,n1}?+ and 
 
nn1n
ii i
uu u
+
? =- (4.97) 
Therefore, Equation (4.95) is a matrix equation Ax=b with  
 
n
23
nn
23 4
n
Q2 Q1 Q
Q1 Q
GC
10
424
CG C
1
24 4 24
CG C
1
24 4 24
CG
01
24 4
??
?
???
++??+
??
??
?? + +? ?+
=
??
?? + +? ?+
??
?
? ?++?
??
A OOO  (4.98) 
 
 
T
n1 n1 n1 n1
12 Q
u,u,..u
++++
? ?=? ? ?
? ?
x  (4.99) 
 
n
nn
0
11
nn
22
nn
Q1 Q1
n
Q1nn
QQ
1C
GD
224
1
GD
2
1
GD
2
C
1
GD
224
??
+
? ??
?+++
? ?
? ?
? ?
?+
? ?
? ?
=
? ?
? ?
?+
? ?
? ?
? ??
?++?
? ?
? ?
b M
 (4.100) 
 
 53
4.2.2 Results and Discussion 
Accuracy 
As in the analysis of the conservative Burgers? equation, we first locate the 
maximum absolute error for all four discrete forms of the TCS method. As shown in 
Figure 4.8, the maximum absolute error occurs at t=0.6 for all four forms. The exact 
solution of the Burgers? equation at t=0.6 can be found in Figure 4.3.  
To evaluate the model error, we examine performance with seven different time 
steps over the range, 0 012 t 0 12..??? , and with grid cell spacing of 
{ }
23 3
x 2 10 5 10 1 25 10,,.
?? ?
?? ? ? ? . As a result, the test conditions cover the range, 
041 CFL 662..??. Similar to section 4.1, the fine-time step solution u% is used in the 
analysis of RMS error, 
2
L  norm and L
?
norm.  Here u%  is the model solution of an 
extremely fine time step such that 
4
t0.6/10?= , which is 0.5% of the smallest  time step 
used in the other simulations. The 
RMS 2
, L norm, and L
?
? norm associated with different 
methods over the tested ranges of CFL numbers are shown in Figure 4.9, 4.10 and 4.11 
respectively. All methods are proven to be 2
nd
-order temporal accurate for fixed x? and 
decreasing t? . For the same x? and t? , it can be seen that the difference in relative error 
among all methods is within one order of magnitude except for the TCSF method. The 
TCSF method has a noticeably smaller error than all the other methods, and about 1 
magnitude smaller error than the TCSG-D form.  
  
 54
0 2 4 6 8 10
0
1
2
3
4
5
6
7
8
x 10
?4
t
?
ab
s
 
 
TCSF
TCSG
TCSF?D
TCSG?D
 
Figure 4.8 Absolute error time evolution for numerical solutions of 1D non-conservative 
Burgers? equation for the TCSF and TCSG methods using 
{ }t 0.01, x 1/ 50, 0.05?= ?= ?=  where 
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
 and 
analytical
u (x, t)  is numerically 
calculated from Equation  (4.46) with K=30. 
 
 
 55
10
0
10
1
10
2
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
10
1
10
2
CFL
?
RMS
 
 
TCSF
TCSG
TCSF?D
TCSG?D
Picard
Newton
Local linearization
2nd order slope
?x=1/50
?x=1/200
?x=1/800
 
Figure 4.9 
RMS
?  vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, 
x? =1/50 and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}. 
 
 56
10
0
10
1
10
2
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
10
1
10
2
CFL
L
2
 
 
TCSF
TCSG
TCSF?D
TCSG?D
Picard
Newton
Local linearization
2nd order slope
?x=1/50
?x=1/200
?x=1/800
 
Figure 4.10  
2
L  norm vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, 
x? =1/50 and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}. 
 
 57
10
0
10
1
10
2
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
10
1
10
2
CFL
L
?
 
 
TCSF
TCSG
TCSF?D
TCSG?D
Picard
Newton
Local linearization
2nd order slope
?x=1/50 ?x=1/200
?x=1/800
 
Figure 4.11 L
?
 norm vs. CFL numbers of various methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, 
x? =1/50 and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}. 
 
 
Stability  
As shown in Figure 4.9, Figure 4.10 and Figure 4.11, all the methods remain 
stable over the full range of tested conditions.  The RK2 and RK4 methods are also used 
to solve this 1D non-conservative Burgers? equation.  The RMS error over a range of 
CFL for the TCSF, TCSG, TCSF-D, TCSG-D, RK2 and RK4 methods is plotted in 
 58
Figure 4.12. In this comparison, 0 0006 t 0 06..??? and x002.? =  provide test 
conditions 0 013 CFL 2 57..??. The RK methods become unstable when CFL>0.1. 
However, the TCS methods remain stable in the whole range of CFL number. 
 
10
?2
10
?1
10
0
10
1
10
?15
10
?14
10
?13
10
?12
10
?11
10
?10
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
?
RMS
 
 
TCSF
TCSF?D
RK2
RK4
2nd order slope
unstable
unstable
 
Figure 4.12 Temporal accuracy of the TCS and RK methods for solution of 1D non-
conservative Burgers? equation at different CFL numbers, where ?=0.05, 
x? =1/50 and t? =0.6/? with ? ?{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 
5}. 
Computational requirement 
The same assessment of operation counts is used in this non-conservative 
Burgers? equation. All the methods perform similar to the non-conservative Burgers? 
equation. The TCSF and TCSG-D have the identical operation counts because they have 
the same solution process only in a reverse mode. The same correspondence can be found 
between the TCSG and TCSF-D for the same reason. In this 1D non-conservative 
Burgers? equation test case, the TCSF and TCSG-D have the least operation counts.  
 
 59
1 2 3 4 5 6 7
0
100
200
300
400
500
600
R
operations per grid point
 
 
TCSF
TCSG
TCSF?D
TCSG?D
local linearization
Newton iteration
Picard iteration
 
Figure 4.13 Ideal operations per grid point for one time step using various nonlinear 
solution methods for 1D Non-conservative Burgers equation. R is the 
number of outer iterations taken by Newton method. It is assumed that the 
Picard method convergences in 
2
R outer iterations. 
 
4.3 APPLICATION OF THE TCS METHOD TO A 1D NONLINEAR ORDINARY 
DIFFERENTIAL EQUATION 
The TCS method can be extended to any equation system with quadratic 
nonlinearities. To test this general capability, we use a nonlinear ODE as our example:   
 
dy
y(1 y) 0
dt
+ ?= (4.101) 
For the initial condition y(0) 1/ 2= , Equation  (4.101) has the analytical solution (Moin 
2001) 
 
[ ]
1
y 1 exp(t)
?
=+  (4.102) 
 
4.3.1 The TCS discretizations  
The TCSF Discretization 
As in the 1D conservative Burgers? equation, the TCSF and TCSG are the same 
for Equation (4.101). Applying the TCSF method to Equation (4.101) provides 
 60
 
()
nn**
t
yy yyy
2
?
?
=+ ?  (4.103) 
 
()
n1 * n1 * *
t
yyyyy
2
++
?
=+ ?  (4.104) 
The discrete form can be written from Equation  (4.103) and  (4.104) as 
 
n
n
y
y
t
1y1
2
?
=
?
? ?? ?
? ?
 (4.105) 
 
**
n1
*
t
yy
2
y
t
1y
2
+
?
?
=
?
?
 (4.106) 
The TCSF-D Discretization 
Similarly, the TCSF-D discrete format for Equation (4.101) is written as 
 
nn
n
t
yy
2
y
t
1y
2
?
?
?
=
?
?
 (4.107) 
 
()
*
n1
*
y
y
t
1y1
2
+
=
?
? ?
 (4.108) 
It can be seen that Equation (4.108) has the same structure as Equation (4.105) and 
Equation (4.107) has the same structure as Equation (4.106). Therefore, the TCSF and 
TCSF-D discrete forms for Equation (4.101) are similar but reversed.  
 
4.3.2 Conventional linearization methods 
Crank-Nicolson (CN) discretized Equation (4.101) is,  
 
() ()
n1nnnn1n1
tt
yy y1y y1y
22
+++
? ?
=+ ? + ?  (4.109) 
To linearize the nonlinear term on right hand side of Equation (4.109), we use the 
Newton iteration, Picard iteration and local linearization methods.  
 
 61
 
Newton method 
Applying the Newton iteration approach, Equation  (4.109) can be written as: 
 
() ()
( )
k
n1
k1 k
n1 n1
n1
k
gy
yy
g
y
+
+
++
+
=?
???
??
?
??
 (4.110) 
where the outer superscript indicates the outer iteration number and 
 
() () ()
n1 n1 n1 n1 n n n
tt
gy y y 1y y y 1y
22
++ + +
? ?
=? ? ?? ? (4.111) 
Substituting Equation (4.111) into Equation (4.110),  the original ODE in the Newton 
linearized CN discretized form is 
 
() ()
() () () ()
()
kkk
n1 n1 n1 n n n
k1 k
n1 n1
k
n1
tt
yy1yyy1y
22
yy
t
12y 1
2
+++
+
++
+
? ?
??
?????
??
??
=?
?
??
??
??
??
 (4.112) 
n1 n
0
yy
+
=  can be used as the initial condition to start the solution. 
 
Picard method 
 A second conventional approach, Picard iteration, in Equation (4.109), results in 
 
() () () ()
k1 k k
n1 n n n n1 n1
tt
yyy1yy1y
22
+
+++
? ?
??
=+ ? + ?
??
??
 (4.113) 
The same initial condition, 
n1 n
0
yy
+
= can be used to start the solution.  
 
Local linearization 
Applying the third conventional method, local linearization method, Equation 
(4.109) can be written as    
 
n
n1 n
n
tf
yy
tf
1
2y
+
?
=+
? ?
????
?
? ?
??
?
??? ?
? ?
 (4.114) 
where  
 ( )fy1y=?  (4.115) 
 62
Substituting Equation  (4.115) into Equation (4.114) provides the locally-linearized 
Equation  (4.109) as 
 
( )
()
nn
n1 n
n
ty 1y
yy
t
12y1
2
+
??
=+
?
? ?
??
? ?
? ?
 (4.116) 
In this specific case, the local linearization, TCSF and TCSF-D are mathematically 
equivalent. If we substitute Equation (4.105) to Equation (4.106) and substitute Equation 
(4.107) to Equation (4.108), we will have two expressions which are exactly the same as 
Equation (4.116).  
    
4.3.3 Comparisons between the TCS and other methods  
The evolution of the absolute errors for two TCS methods applied to the ODE is 
shown in Figure 4.14.  The maximum absolute errors for both TCS discrete formats occur 
at time t=1.5.  
 
0 2 4 6 8 10
0
1
x 10
?4
t
?
abs
 
 
TCSF
TCSF?D
 
Figure 4.14 Time evolution of absolute errors for TCS methods applied to the ODE, 
Equation  (4.101), for t? =0.1, where 
abs model analytical
(t) y (t) y (t)?= ? .  
 
Figure 4.15 provides a comparison of root-mean-square (RMS) error over ? time 
steps prior to maximum error, defined as 
 63
 
1/2
2
RMS mod el i analytical i
i1
1
y (t) y (t)
?
=
??
??
???? ?
??
??
?
??
?
 (4.117) 
 where 
max
t1.5=  and { }10, 50, 100, 200, 500, 600, 800, 1000, 1200, 1500, 2000?? .  All 
tested methods provide second-order accuracy. In this particular case, the iterative 
methods have slightly smaller error than non-iterative methods. The errors from the three 
non-iterative methods are identical since TCSF, TCSF-D and local linearization are 
mathematically equivalent. 
 
10
1
10
2
10
3
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
?
?
RMS
 
 
TCSF
TCSF?D
Picard
Newton
Local linearization
2nd order
 
Figure 4.15 RMS error from Equation  (4.117) for discrete solutions of the ODE, 
Equation (4.101), computed using a range of ? time steps. 
 
4.4 SUMMARY 
The TCS method has been tested on three different test cases: a 1D conservative 
Burgers? equation, a 1D non-conservative Burgers? equation and an ODE.  The TCS 
method has multiple discrete forms for each test case. All the TCS discretizations have 
been compared to conventional implicit linearization methods, including Picard iteration, 
Newton iteration and local linearization methods.  All the TCS discretizations 
 64
demonstrate 2
nd
-order temporal accuracy and are shown stable up to a CFL O (10).   All 
the TCS discretizations are proved stable when the CFL value is in the order of 10.  The 
stability advantage of the TCS methods has been further demonstrated by the comparison 
with the RK methods. The TCS methods require fewer operation counts than Picard and 
Newton method, but are similar to local linearization.  The principal advantage of the 
TCS method over local linearization is the relative ease with which the TCS method can 
be derived and implemented as it does not require discrete evaluation of a function 
Jacobian. The TCS method is applicable to any quadratically-nonlinear problems, which 
is verified by the ODE case.   
For each test case, different TCS discretizations perform similarly in temporal 
accuracy, stability and efficiency. However, difference can be observed among different 
discretizations. In the 1D conservative Burgers? equation, the TCSF and TCSG collapse 
into one single form. The TCSF-D reverses the solution process of the TCSF. 
Nevertheless, the TCSF-D has a noticeably smaller relative error than the TCSF and other 
methods. In the 1D non-conservative Burgers? equation, the TCSF and TCSG hold 
distinct forms. The TCSF-D reverses the solution process of the TCSG, and the TCSG-D 
reverses the solution process of the TCSF. Among the four discretizations, the TCSF has 
the smallest relative error, which is about one order of magnitude smaller than the TCSG-
D. In the 1D nonlinear ODE case, the TCSF and TCSG remain the same. The TCSF-D 
reverses the solution process of the TCSF, but they have the same relative error.    
 
 65
Chapter 5 The TCS Family Method 
 
In Chapter 3, we revealed that the key to the TCS method is the 2nd-order 
accurate time-centered splitting. Applying this splitting to different terms, different TCS 
discretizations can be generated. In Chapter 4, we compared the TCS discretizations with 
conventional linearization methods using various examples. In this chapter, we derive the 
TCS family method in one general form by introducing weighting factors to different 
terms. A significance of this TCS family method is that the value of the weighting factors 
does not affect the order of accuracy of the scheme. Since weighting factors can be any 
value between zero and one while still providing essentially the same advantages, it may 
provide a flexible approach for a variety of problems. To examine the properties of 
weighting factors, the 1D non-conservative Burgers? equation is solved using different 
combinations of weighting factors. Furthermore, we demonstrate that the TCS method 
can computationally decouple an equation system of coupled variables using special 
combinations of weighting factors. This advantage makes the TCS method more efficient 
in solving coupled multi-variable equation systems.  
 
5.1 DERIVATION OF THE TCS FAMILY METHOD 
Let us start a time marching differential equation from a simplest form: 
 ()() ()
{}
fFM L N
t
??
= ?? ???
??
?
,  (5.1) 
where ? s a generic variable, F, M, L and N represent generic linear operators. The 
function f is a linear function of ( ) ( )FM L? ?? ?
? ?
 and ( )N ? . Equation (5.1) then 
represents a time marching differential equation with a quadratic nonlinearity of 
()()FM L????
??
 and a linear operator of the variable, ( )N ? . This quadratic nonlinearity 
can be: the product of the variable itself, 
2
?  if F, M and L are equal to 1; the product of 
the variable and the gradient of the variable, / x??? ? ,  if F=1, M=1 and L= ()/x??; the 
gradient of the product of the variable, ( )/x? ?? ?  if  M=1, L=1 and F= ()/x??. Thus, 
Equation (5.1) represents a generic time marching differential equation with any kind of 
quadratic nonlinearity. Equation (5.1) can be discretized using the midpoint rule, written 
as  
 
( ) ( ) ( )
{ }
( )
n1 n n12 n12 n12 3
fFM L N t t
++++
??
? =? + ?? ??+ ?
??
// /
,O (5.2) 
In the previous chapters, we demonstrated that the time-centered split can be introduced 
to the quadratic nonlinear term and/or the linear term. Introducing two weighting factors, 
 66
1
? and 
2
? , and combining with the time-centered linear approximation, 
()
nn1
n12 2
t
2
/
O
+
+
?+?
?= +? , Equation (5.2) can be written as a general discretization 
 
()()()
()() ()
nn1 nn1
n1 n n12 n12
11
nn1
n12 3
22
fFM L 1 FM L
N1N tt
2
++
+++
+
+
?
????
?? ???+? ?+?
?
?=?+? ? +?? ?
? ????
?? ??
?? ???
????
?
?
???+?
?? +?? ?+?
?
??
??
?
//
/
,
O
(5.3) 
 
where the weighting factor 
i
?  ( { }i12,? ) is required that 
i
01??? . Equation (5.3) is 
mathematically equivalent to Equation (5.2). Similar to the derivation in Chapter 3, an 
intermediate variable 
*
? is introduced as 
 
()( ) ()( )()
{
()()()}
nnn12 n12n
11
n12 n
22
fFM L 1 FM L
t
N1N
2
++
+
????
?=?+ ? ? ? + ?? ? ?
????
?
?? +?? ?
*//
/
,
 (5.4) 
Subtracting Equation (5.4) from Equation (5.3), the second step is written as: 
 
()( ) ()( )( )
{
()()()} ()
n1 n1 n12 n12 n1
11
n12 n1 3
22
fFM L 1 FM L
t
N1N t
2
+++ ++
++
????
?=?+? ? ? +?? ? ?
????
?
?? +?? ? +?
*//
/
,
O
 (5.5) 
Expanding 
n12/+
? in Equation (5.4) using the Taylor expansion,  
 
()
() ()
()()
2
nnn
1t
2
nn
1t
2
2t 2
tt
fFM L
22
tt
1FM L
22
tt t
N1N
22 2
?
??
??
???
??
?=?+ ? ? ?+? +??
???
??
??
? ??
???
??
??
??
+?? ?+? + ?
??
??
??
??
??
?
??
? ???
??
??+?+ +?? ?
???
??
??
???
?
*
O
O,
O
 (5.6) 
Organizing Equation (5.6)  and grouping the higher order terms, 
 
()()
{
()}
2
nnnn
tt
fFM L N
22
??
??
??
?=?+ ? ? ? +
??
??
??
*
,O (5.7) 
n12/+
? can be approximated using the explicit Euler approximation as 
 67
 
()() ()
{}
2
n12 n n n n
tt
fFM L N
22
+
? ?
??
??
?=?+ ?? ? +
??
??
??
/
,O (5.8) 
Therefore, we obtain the correspondence between 
*
? and 
n12/+
? from Equations (5.7) and 
(5.8) such that 
 
( )
n12 2
t
*/
O
+
? =? + ?  (5.9) 
Substituting Equation (5.9) into Equations (5.4) and (5.5), the original equation set 
becomes a computationally linearized two-step equation system: 
 
( ) ( ) ()( ) ( )
{
()()()}
nn n
11
n
22
fFM L 1 FM L
t
N1N
2
????
?=?+ ? ? ? + ?? ? ?
????
?
??+?? ?
***
*
,
 (5.10) 
 
( ) ( ) ()( ) ( )
{
()()()} ()
n1 n1 n1
11
n1 3
22
fFM L 1 FM L
t
N1N t
2
++ +
+
????
?=?+? ? ?+?? ? ?
????
?
??+?? ? +?
***
*
,
O
 (5.11) 
The weighting factor, 
i
? , distinguishes this derivation above from the original derivation 
of TCS method in Chapter 3 because the value of  
i
? will not affect the 2
nd
-order accuracy 
of the TCS methods. In addition, there is no dependent relationship between 
1
? and 
2
? . In 
other words, as long as 
i
01??? , Equations (5.10) and (5.11) are a 2
nd
-order equivalent 
to Equation (5.2). Although the derivation above uses a single variable, the same 
mathematic principles apply to any quadratic nonlinearity. In ordinary differential 
equations (ODE), the nonlinearity is simply 
2
x  or xy; in partial differential equations 
(PDE), this quadratic nonlinearity can be UU/ x? ?  or ( )HU / x? ? . Therefore, the 
derivation above is applicable to an ODE or ODE system with quadratic nonlinearity and 
typical flow transport equation such as 1D advection-diffusion equation or shallow water 
equations. In addition, the more variables are involved in an equation system, the more 
weighting factors that can be introduced. As a result, the TCS method provides not a 
single kind, but a whole family of TCS discretizations.  
As we illustrated in the previous chapters, all the TCS family methods keep 2
nd
-
order temporal accuracy and they share the same simplicity advantages in computation. 
However, with different combinations of 
i
? , the numerical solutions perform differently. 
In the next section, we will use the 1D non-conservative Burgers? equation as an example 
to explore the different performance of different TCS discretizations.  
 
 68
5.2 APPLICATION OF THE TCS FAMILY METHOD TO THE 1D NON-CONSERVATIVE 
BURGERS? EQUATION 
With the introduction of the weighting factors 
1
?  and 
2
? , the general format of 
the TCS discretized non-conservative Burgers? equation is written as, 
 
()()()()(){}
nn n2 2n
1x 1 x 2x 2x
t
uu uu 1 uu u 1 u
2
****
?
=+ ??? +???? +??? +????  (5.12) 
 
 
()()() ()
n1 n1 n1 2 2 n1
1x 1x 2x 2x
t
uu uu1 uu u1 u
2
****++ + +
?
=+ ??? +???? +???+???? (5.13) 
The above equations are a two-step matrix equation system in the format of Ax=b. Using 
a central difference in spatial derivatives, the LHS coefficient matrix A in the first step is 
tridiagonal such that 
 
1,1 1,2
2,1 2,2 2,3
3,2 3,3 3,4
Q2,N3 Q2,Q2 Q2,Q1
Q1,N 2 Q1,Q1 Q1,Q
Q,Q 1 Q,Q
aa 0
aaa
aaa
aaa
aa a
0a
?? ?? ??
?? ?? ?
?
??
??
=
??
A OOO  (5.14) 
 
where Q is the number of the grid points in space, and the non-zero components of A  for 
{ }i 1,2...Q=  are of the form 
 
n
12
i
nn
i,j 1 i 1 i 1 2
n
12
i
C:ji1
42
1
a1 CC :ji
4
C:1
42
+?
??
?
??? =?
?
?
?
??=+? ? +?? =
?
??
?
??
?
+?? =+
?
?
 (5.15) 
  
Also in the first step,  
 
T
** *
12 Q
u ,u ,...u? ?=
? ?
x  (5.16) 
and 
 
 
 69
 
( )
()
()
()
n
n
n2111 2
1
n
n22
2
n
2Q1n
Q1
nn
2Q1Qn
2
Q
1D
C
u
24
1D
u
2
1D
u
2
1DC
u
24
?
?
???? ?
? ??
+++
??
?? ?
+
=
??
?? ?
+
?? ? ?
? ?
??
+?+
??
b M  (5.17) 
 
In the second step, we also have a tridiagonal coefficient matrix A with different 
tridiagonal elements as 
 () ()
*
12
i
**
i,j 1 i 1 i 1 2
*
12
i
11
C:ji1
42
1
a11 CC 1 :ji
4
11
C:1
42
+?
?? ???
??? =?
?
?
?
??=+ ?? ? +??? =
?
??
?
?? ??
?
+?? =+
?
?
 (5.18) 
 
and 
 
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?=
? ?
x  (5.19) 
 
 
 
( ) ( )
()()
*
*
*112
21
1
*
*
22
2
*
2Q1*
Q1
**
2Q 1Q*2
Q
1C1
D
u
24
D
u
2
D
u
2
D1 C1
u
24 2
?
?
???????
??
++ +
??
??
+
=
??
??
+
?? ?? ???
+? +
??
b M  (5.20) 
 
where C,? and D  in these two steps have the same definitions as in Chapter 4.  
 
 70
5.3 RESULTS AND DISCUSSION 
To explore the properties of different weighting factors, the 1D non-conservative 
Burgers? equation is solved using Equations (5.12) and (5.13) with different 
combinations of 
1
? and 
2
? . Before we analyze the results from various combinations of 
1
? and 
2
? , it is necessary to review their meanings. In the TCS scheme, 
n12
u
/+
 in the 
Burgers? equation is either approximated using the time-centered splitting as 
()()
n12 n n1 2
uuu2t
/
/O
++
=+ +? or u
*
 as 
( )
n12 2
uu t
/*
O
+
= +?. Different 
i
?  indicates 
different approximations of 
n12
u
/+
. When 
i
? =0 or 
i
? =1, 
n1/2
u
+
 is approximated only 
using one kind of the approximations. When 
i
?  lies between 0 and 1, 
n1/2
u
+
 is 
approximated by the combinations of these two approximations. Table 5.1 lists the 
mathematical meaning of different values of 
i
? . 
 
Table 5.1 The mathematical meaning of the value of 
1
? and 
2
? . 
value 
1
?  
2
?  
0 
split the entire gradient term and 
approximate the entire flux term using 
u
*
 
split the entire diffusion 
term 
0and1<< 
split part of the gradient and flux terms 
approximate part of the gradient and flux 
terms using u
*
 
split part of the diffusion 
term and approximate part 
of the diffusion term using 
u
*
 
1 
split the entire flux term and 
approximate the entire gradient term 
using u
*
 
approximate the entire 
diffusion term using u
*
 
 
 
5.3.1 Accuracy of Different Weighting Factors 
The error analysis of different TCS discretizations is evaluated in the same vein as 
in Chapter 4. In figure 5.1, the RMS error is plotted against a range of CFL numbers for 
six different 
1
? values ( { }
1
0131223341,/,/, /,/,?? ). For each 
1
? value, four 
different 
2
?  values ( { }
2
013121,/,/,?? ) are chosen for comparison, resulting 
twenty four different TCS discretizations. 
 The total marching time t=0.6 is used as the basis for error comparison in Figure 
5.1 because the maximum absolute error occurs at t=0.6 (as demonstrated in Chapter 4). 
The CFL number is defined using 
analytical
u, x150/? =  and t? = 06./? , where ? is the 
 71
number of the time steps and { }5, 8, 10, 20, 25, 30, 50?? . Thus, the tested CFL number 
covers the range of  04 CFL 42..??.  It can be observed that every TCS discretization 
has the same 2nd-order temporal accuracy and remains stable in the whole range of the 
CFL numbers. The computation requirement of each TCS discretization is essentially the 
same since they share the same general structure of matrix equations as in Equations 
(5.14) to (5.20). Figure 5.1 clearly shows that for a given value of 
1
? , any choice of 
2
?  
gives an error of the same magnitude. The reverse also holds true, as shown in Figure 5.2.  
The RMS error is plotted for six values of 
2
? , choosing four 
1
?  values for each 
2
?  value. 
Again, given a value for 
2
? , the choice of 
1
?  affects the RMS error less than or equal to 
one order of magnitude.  
Though 
12
23 12/, /?= ?=  gives the lowest RMS error in these two examples, 
the best combination for the smallest error is problem-specific; error arises from 
approximations of different terms. To obtain optimal combinations of 
i
? for the lowest 
error, one needs to examine the nature of the solution and each term at the designated 
marching time. For this Burgers? equation, the solution u, the gradient term ux/??,  and 
the diffusion term 
22
ux/?? ?  all need to be taken into consideration because 
approximations have been introduced to all of these terms.  
 
 
 
 72
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1
 
 
?
2
=0 ?
2
=1/3 ?
2
=1/2 ?
2
=1
2nd order slope
 
Figure  5.1 Temporal accuracy of various combinations of 
1
?  and 
2
? for solution of the 
Burgers? equation at t=0.6, where x150/? = , 005.?=  and t06./?= ? 
( { }5, 8, 10, 20, 25, 30, 50?? ).  
 
 73
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1
 
 
?
1
=0 ?
1
=1/3 ?
1
=1/2 ?
1
=1
2nd order slope
 
Figure  5.2 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=0.6, where x150/? = , 005.?=  and t06./?= ? 
( { }5, 8, 10, 20, 25, 30, 50?? ).  
 
 
 74
5.3.2 Stability of Different Weighting Factors 
Figure 5.1 and 5.2 display the temporal accuracy of the solution of the Burgers? 
equation at t=0.6. All TCS discretizations produce stable solutions for the whole range of 
tested CFL numbers. To further explore the stability characteristics, we investigate the 
performance of the TCS family method at t=0.3, 1 and 2, and these results are shown in 
Figure 5.3, 5.4 and 5.5, respectively. Similar to Figure 5.1, the RMS error is plotted 
against a range of CFL numbers for several choices of  
2
?  given a 
1
?  value. Similar to 
t=0.6, all 24 TCS discretizations stay stable at t=0.3.  
At t=1and t=2, though, we find different results. At t=1 (Figure 5.4), for any given 
1
?  value, the solution of 
2
1?=becomes unstable with the increase of the CFL number, 
but the solutions for any 
2
1?? remain stable. At t=2, the solutions of both 
2
0?= and 
2
1?=become unstable with the increase of the CFL number, but the solutions with any 
2
01<? <  remain stable.  
To understand this stability property associated with the 
2
? value, we need to 
examine the nature of the solution and the meaning of the weighting factor 
2
? . The 
solution of the Burgers? equation evolving with time is plotted in Figure 4.1. It can be 
observed that after t=0.6, damping dominates the solution, which suggests that the 
diffusion effect dominates.  The weighting factor 
2
?  controls the approximation of the 
diffusion term as shown in Table 5.1.  If 
2
0? = , the entire diffusion term, 
2n1/2 2
u/x
+
?? ? , 
is completely approximated using the time-centered split, which means  
n1/2
u
+
in the 
diffusion term is first approximated as 
( )
n12 n n1
uuu2
/
/
++
?+  and then we split the 
diffusion terms into two steps as 
2n 2
u/x??? and 
2n1 2
u/x
+
???. If  
2
1?=, the entire 
diffusion term, 
2n1/2 2
u/x
+
?? ? ,  is completely approximated  using u
*
, which means 
2n1/2 2 2* 2
u/x u/x
+
?? ? ??? ? . Figures 5.1 to 5.5 indicate that the TCS method remains 
stable and the stability is not affected by the values of the weighing factors. If the 
diffusion dominates the solution, the weighting factor associated with the diffusion term 
can affect the stability. Our results suggest that if the weighting factor, 
22
0and 1?? ??, 
the solutions remain stable for the tested CFL numbers. In other words, using only one 
kind of approximation for the diffusion term, the solution may become unstable with the 
increase of the CFL number; using a combination of the two kinds of approximations for 
the diffusion term can stabilize the solution.  
 
 
 75
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1
 
 
?
2
=0 ?
2
=1/3 ?
2
=1/2 ?
2
=1
2nd order slope
 
Figure  5.3 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=0.3, where x150/? = , 005.?=  and t03?= ?./  
( { }5, 8, 10, 20, 25, 30, 50?? ).  
 
 76
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1
 
 
?
2
=0 ?
2
=1/3 ?
2
=1/2 ?
2
=1
2nd order slope
 
Figure  5.4 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=1, where x150? = / , 005?= .  and t1?= ?/  
( { }5, 8, 10, 20, 25, 30, 50?? ).  
 77
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1
 
 
?
2
=0 ?
2
=1/3 ?
2
=1/2 ?
2
=1
2nd order slope
 
Figure 5.5 Temporal accuracy of various combinations of 
2
? and 
1
? for solution of the 
Burgers? equation at t=2, where x150/? = , 005.?=  and t2?= ?/  
( { }5, 8, 10, 20, 25, 30, 50?? ).  
 78
The general form of TCS method is demonstrated using the 1D non-conservative 
Burgers? equation. The same application of the weighting factors can be used for an 
equation system with coupled variables. Hence, a general form of a two-step linearized 
coupled equation system can be created using the same principle. Since the values of the 
weighting factors will not affect the 2
nd
-order temporal accuracy of the scheme, TCS 
method may provide a flexible approach to solve different problems. Furthermore, 
properly chosen weighting factors in the TCS method can computationally decouple a 
system of coupled equations. This decoupling advantage of the TCS method is analyzed 
in the following section.  
 
5.4 COMPUTATIONAL DECOUPLING  
In the previous sections, we demonstrated the effects of weighting factors in the 
TCS method on the accuracy and stability of the solution. In this section, we will 
introduce the decoupling advantage of the TCS method. The two-dimensional Burgers? 
equation provides a convenient example: 
 
22
22
uuu u u
uv
txyx y
???? ?
++=?+?
???? ?
 (5.21) 
 
22
22
vvv v v
uv
txyx y
???? ?
++=?+?
???? ?
 (5.22) 
 
Equations (5.21) and (5.22) are a two-variable two-dimensional system with coupled 
variables ?u? and ?v?. If we use the 2
nd
-order temporal accurate C-N method to discretize 
equations (5.21) and (5.22), the resulting equations become 
 
n n 2n 2n
n1 n n n
22
n1 n1 2 n1 2 n1
n1 n1
22
tuuuu
uu u v
2xyxy
tu vuu
uv
2xyxy
+
++++
++
???? ???
= + ? ? +? +?
??
????
??
?? ???
+? ? +? +?
??
 (5.23) 
 
 
nn2n2n
n1 n n n
22
n1 n1 2 n1 2 n1
n1 n1
22
tvvvv
vv u v
2xyxy
tv vvv
uv
2xyxy
+
++++
++
???? ?? ?
= + ? ? +? +?
??
????
??
?? ?? ?
+? ? +? +?
??
 (5.24) 
Equations (5.23) and (5.24) are a coupled nonlinear equation system. Conventionally, to 
solve a coupled equation system, either a simultaneous solution or a sequential solution 
method can be used (Ferziger and Peric, 1996). However, such an equation system is too 
complex and expensive to solve simultaneously due to its nonlinearity. Therefore, 
sequential solution methods are more appropriate. To solve this equation system 
 79
sequentially, inner iterations are required to solve Equations (5.23) and (5.24) 
individually, and outer iterations are required to satisfy Equations (5.23) and (5.24) 
together.  
Alternatively, we can decouple and linearize Equations (5.23) and (5.24) using 
TCS method without iterations by choosing specific weighting factors. Applying the TCS 
method to Equation (5.21) results in a two-step system with a specific 
i
? weighting factor 
for each term 
 
()()(){
()()()
()()
()()}
nn n
1x 1 x
nn
2y 2 y
22n
3x 3x
4y 4y
t
uu uu 1 uu
2
vu 1 vu
u1 u
u1 u
***
**
*
*
?
=+ ??? +????
+? ? ? + ?? ? ?
+? ?? + ?? ??
+? ?? + ?? ??
 (5.25) 
 
 
()()(){
()()()
()
() }
n1 n1 n1
1x 1x
n1 n1
2y 2y
22n1
3x 3 x
4y 4 y
t
uu uu1 uu
2
vu 1 vu
u1 u
u1 u
***
**
*
*
++ +
++
+
+
?
=+ ??? +????
+? ? ? + ?? ? ?
+??? + ?? ??
+? ?? + ?? ??
 (5.26) 
 
Applying the TCS method to Equation (5.22) results in a two-step system: 
 
()()(){
()()()
()()
()()}
nn n
5x 5 x
nn
6y 6 y
22n
7x 7x
8y 8y
t
vv uv 1 uv
2
vv 1 vv
v1 v
v1 v
***
**
*
*
?
=+ ??? +????
+? ? ? + ?? ? ?
+? ?? + ?? ??
+? ?? + ?? ??
 (5.27) 
 
 
()()(){
()()()
()
() }
n1 n1 n1
5x 5x
n1 n1
6y 6y
22n1
7x 7 x
8y 8 y
t
vv uv1 uv
2
vv 1 vv
v1 v
v1 v
***
**
*
*
++ +
++
+
+
?
=+ ??? +????
+? ? ? + ?? ? ?
+??? + ?? ??
+??? + ?? ??
 (5.28) 
For this two-variable two-dimensional Burgers? equation, eight different weighting 
factors are applied in the general TCS discretized form. Observing Equations (5.21) and 
(5.22), we discover that the coupling only occurs in the nonlinear advection term. 
 80
Therefore, the weighting factors for the diffusion terms will not affect the computational 
decoupling.  
In the following discussion, we omit the diffusion terms and their associated 
weighting factors. We then discover that the coupling occurs between the nonlinear 
advection terms ? vu y/??? and ? uv x/? ? ?. The two weighting factors for these two terms 
are 
2
?  and 
5
?  respectively and these two factors are the key to computationally 
decoupling the equation system. The other two advection terms ? uu x/??? and 
?v v y/??? are not coupled. Therefore, we have more choices for their weighting factors 
1
?  and 
6
? . To completely decouple u and v, Table 5.2 shows the possible combinations 
of the weighting factors for each advection term.  
 
Table 5.2 Weighting factors for computational decoupling the 2D Burgers? equation 
2
?  
5
?  
1
?  
6
?  
0 0 
1 0 
0 1 
0 1 
1 1 
0 0 
1 0 
0 1 
1 0 
1 1 
 
It can be seen from Table 5.2 that two kinds of combinations of 
2
?  and 
5
? can decouple 
the 2D Burgers? equation. For each pair of 
2
?  and 
5
? , four different combinations of 
1
? and 
6
?  can be chosen such as 
16
0? =? = , 
16
1and 0? =?=, 
16
1?=?= , 
16
0and 1?= ?= . To illustrate the procedure of decoupling, we use the values from the 
first row in Table 5.2, applying 
2516
0100,,,? =?=?=?= to Equations (5.25) to 
(5.28) and reordering the equation system to 
 
()
nn n
xy
t
vv uvvv
2
***
?
=+ ????  (5.29) 
 
()
nnn
xy
t
uu uuvu
2
***
?
=+ ????  (5.30) 
 
()
n1 n1 n1
xy
t
u u uu vu
2
** *+++
?
=+ ?? ??  (5.31) 
 81
 
()
n1 n1 n1
xy
t
vv uvvv
2
***+++
?
=+ ????  (5.32) 
Each step from Equation (5.29) to (5.32) is computationally linearized, and each 
step has only one unknown variable. Variables ?u? and ?v? are effectively decoupled in 
each step. The solution process for these chosen weighting factors is as follows: 
1. Solve for v
*
 using 
n
u  and 
n
v   
2. Solve for u
*
 using 
n
u  and v
*
 obtained from step 1 
3. Solve for 
n1
u
+
 using u
*
 obtained from step 2 and v
*
 obtained from step 1 
4. Solve for 
n1
v
+
 using 
n1
u
+
 obtained from step 3 and v
*
 obtained from step 1 
A different combination of weighting factors from Table 5.2 results in a different 
equation system, and the order of the solution procedure is changed appropriately. Using 
the following set of weighting factors in Equations  (5.25) to (5.28): 
2
1?=, 
5
0?= , 
1
0?= , and 
6
0?= , results in the following system: 
 
()
nnn
xy
t
uu uuvu
2
***
?
=+ ????  (5.33) 
 
()
nnn
xy
t
v v uv vv
2
***
?
=+ ????  (5.34) 
 
n1 n1 n1
xy
t
v v uv vv
2
** *+++
?
=+ ?? ??  (5.35) 
 
()
n1 n1 n1
xy
t
uu uuvu
2
** *
?
=+ ?? ? ?  (5.36) 
The corresponding solution process is:  
1. Solve for u
*
 using 
n
u  and 
n
v  
2. Solve for v
*
 using 
n
v  and u
*
 obtained from step 1 
3. Solve for 
n1
v
+
 using v
*
 obtained from step 2 and u
*
 obtained from step 1 
4. Solve for 
n1
u
+
 using 
n1
v
+
 obtained from step 3 and u
*
 obtained from step 1 
This computational decoupling technique can be extended to equation systems 
with more than two variables. The underlying principle for decoupling is that when the 
weighting factor is switched from 0 to 1, the approximation for the variable 
n12/+
?  is 
changed from  
*
?   to either  
n
?  or 
n1+
? . Consequently, we can change the variable in one 
equation from unknown to known.  
 
5.5 SUMMARY 
A general form of the TCS family method is developed in this chapter by 
introducing a weighting factor 
i
?  for each term in a quadratic nonlinear differential 
equation. The 1D non-conservative Burgers? equation is used as an example to test how 
the weighting factors affect the solutions. The value of 
i
?  controls the approximations of 
 82
the advection and diffusion terms.  If 
i
? =0 or
i
? =1, the advection term or the diffusion 
term is approximated only by one kind of approximation; if 
i
?  lies between 0 and 1, the 
advection term or the diffusion term is approximated by the combination of  two 
approximations. Although various TCS discretizations can be generated by changing 
i
?  
the 2
nd
 order temporal accuracy of the TCS method is not affected by the choice of the 
weighting factor. We proved analytically and by examples that the TCS family method is 
2
nd
-order temporal accurate if 
i
01??? .  Different combinations of 
i
? value have been 
tested using 1D Burgers? equation as an example. The results presented here show that 
the optimum combination of the weighting factors for the most accurate solutions is 
problem-specific. However, our results indicate that the optimal values for the weighting 
factors lie between 0 and 1 such that the advection term and diffusion term are 
approximated using the combination of two different approximations. Our examples also 
demonstrate the stability advantage of the TCS method. When the advection effect is 
dominant, the stability is not affected by 
i
?  value and all the tested TCS discretizations 
remain stable for the whole range of CFL numbers. When the diffusion effect is 
dominant, the TCS method stays stable under all 
i
?  values except 0 or 1 for the diffusion 
term. This observation suggests that the solution from the TCS method may be stabilized 
by adjusting the weighting factors. Most of the problems in Environmental Fluid 
Mechanics are dominated by advection. Hence, the TCS family method has a promising 
stability advantage in simulating environmental flows.  
In addition to its accuracy and stability advantages, the TCS method is also 
capable of decoupling coupled equation systems by choosing specific combinations of 
the weighting factors. This advantage is displayed using the two-dimensional Burgers? 
equation as an example. Compared to the conventional decoupling techniques, no 
iteration is required in the TCS method. Therefore, the TCS method is more efficient in 
solving coupled nonlinear equation systems. Moreover, the TCS family method provides 
theoretically unlimited possible discretizations for one problem. This, in turn, provides 
flexibility in discretizing specific problems.  
In the next chapter, we will apply the TCS method to a multi-variable and multi-
dimensional equation system: depth averaged shallow water equations (SWE). The 
advantage of computational linearization and decoupling of the TCS method will be 
further illustrated in the next chapter. Numerical experiments will be carried out to test 
the TCS discretized depth averaged SWE.  
 
 83
Chapter 6 Application of the TCS Method to a 2D Depth Averaged 
Shallow Water Equations (SWE) 
 
The TCS (Time-Centered Split) family of methods are derived and analyzed in 
Chapter 5. To investigate the application of the TCS method in a multi-variable and 
multi-dimension equation system, depth averaged SWE (shallow water equations) is 
solved using the TCS method in this chapter. The computational decoupling 
characteristics of the TCS method are explored using depth averaged SWE as an 
example. The decoupling procedure is discussed and three representative linearized and 
decoupled TCS discretizations are presented in this Chapter. Numerical experiments such 
as a one-dimensional standing wave and a two-dimensional standing wave in a 
rectangular domain are performed to verify the TCS numerical model. In the one-
dimensional standing wave case, the numerical results are compared with the analytical 
solutions; and in the two-dimensional standing wave case, the characteristics and 
performance from three TCS discretizations are compared and discussed.  
 
6.1. THE 2D DEPTH AVERAGED SWE 
The 2D depth averaged SWE are obtained by integrating the 3D incompressible 
Navier-Stokes equations over the water depth with the two following assumptions: 1) 
neglecting the vertical velocity and acceleration; 2) applying a hydrostatic pressure 
distribution (Vreugdenhil 1994). These equations are widely used in hydraulic science 
and engineering (Ancey et al. 2008; Arega et al. 2008; Hunter et al. 2008). The definition 
of ?shallow? requires that the vertical scale of the flow is small compared to the 
horizontal scale. In nature or common engineering practice, many types of flows fall into 
this category such as atmospheric flow (Mohebalhojeh and Dritschel 2007), river flow 
(Arega et al. 2008), and storm surge (Bajo et al. 2007). Although the 2D SWE cannot 
simulate vertical velocity gradients, they are useful for flows where strong turbulence 
provides complete vertical mixing of momentum or for flows dominated by barotropic 
motions. The 2D shallow water equations can be written as: 
 
22
22
UUU UU
UV g
txy xxy
????? ????
++=?+?+
??
??? ???
??
 (6.1) 
 
22
22
VVV VV
UV g
txy yxy
????? ????
++=?+?+
??
??? ???
??
 (6.2) 
 
HHUHV
0
tx y
? ??
+ +=
?? ?
 (6.3) 
 
b
HZ?=+ (6.4) 
where U  and V  are depth-averaged velocities in x  and y  directions; g is the 
gravitational acceleration; H  is the water depth; ?  is the water surface elevation; and 
 84
b
Z is the bottom elevation.  The above statement of the 2D SWE neglects the turbulence 
closure and bottom stress for the test case of the TCS numerical method. 
 
6.2 THE TCS DISCRETIZED SWE   
Applying the midpoint rule between time n and n+1 along with central-
differencing of spatial derivatives on an Arakawa C grid (Arakawa and Lamb, 1977), as 
shown in Figure 6.1, Equations (6.1), (6.2) and (6.3) can be discretized as a set of coupled 
nonlinear equations: 
 
()( )
()
n1/2 n1/2 n1/2 n1/2
i,j 3/2 i,j 1/2 i 1,j 1/2 i 1,j 1/2
n 1 n n 1/2 n 1/2
i,j1/2 i,j1/2 i,j1/2 i,j1/2
n1/2 n1/2 n1/2
i,j 1/2 i,j 1/2 i,j 3/22
n1/2 n
i 1,j 1/2 i,j 1/22
UU U U
UU tU V
2x 2y
U2UU
x
U2U
y
++ + +
+? ?
++
+ +
+++
?+
++
?+ +
?
??
?
=+?? ?
?
??
?
?
?
+?
?
?
+?
?
()
1/2 n 1/2
i1,j1/2
n1/2 n1/2
i,j 1 i,j
U
g
x
+
++
++
+
+
?
????
?
?
?
 (6.5) 
 
( ) ( )
()
n1/2 n1/2 n1/2 n1/2
i 1/2,j 1 i 1/2,j 1 i 3/2,j i 1/2,j
n1 n n1/2 n1/2
i 1/2,j i 1/2,j i 1/2,j i 1/2,j
n1/2 n1/2 n1/2
i 1/2,j 1 i 1/2,j i 1/2,j 12
n1/2
i 1/2,j i 1/2,j2
VV VV
VV tU V
2x 2y
v
V2VV
x
v
V2V
y
++ ++
++ +? + ?
++
+ +
+++
+? + ++
+
?+
?
??
?
=+?? ?
?
??
?
?
+?
?
+?
?
()
n1/2 n1/2
i1/2,j
n1/2 n1/2
i1,j i,j
V
g
y
++
+
++
+
+
?
????
?
?
?
 (6.6) 
  
()
n 1 n n 1/2 n 1/2 n 1/2 n 1/2
i,j i,j i,j 1/2 i,j 1/2 i,j 1/2 i,j 1/2
n1/2 n1/2 n1/2 n1/2
i 1/2,j i 1/2,j i 1/2,j i 1/2,j
1
HH t HU HU
x
1
HV HV
y
++++
++ ??
++ ++
??
?
=+?? ?
?
?
?
?
??
?
?
?
 (6.7) 
 85
  
Figure 6.1 Illustration of Arakawa C grid.  
 
As illustrated in Figure 6.1, Arakawa C grid is one kind of staggered grid. Cartesian 
staggered grid is standard in modeling SWE (Schoenstadt 1980; Stelling and Duinmeijer 
2003). In Equations (6.5) to (6.7), the discrete value of the velocity component U is 
defined at (i, j+1/2); and the value of the velocity component V is defined at (i+1/2, j); 
and the discrete value of surface elevation ?  is defined at (i, j).  In Equation (6.5), V (i, 
j+1/2) is obtained by the linear interpolation of the values: V (i+1/2, j), V (i-1/2, j), V 
(i+1/2, j+1), and V (i-1/2, j+1). In Equation (6.6), U (i+1/2, j) is obtained by the linear 
interpolation of the values: U (i, j+1/2), U (i+1, j+1/2), U (i, j-1/2), and U (i+1, j-1/2). In 
Equation (6.7), H (i, j+1/2) is obtained by the linear interpolation of the values of H (i, j) 
and H (i, j+1). H (I, j-1/2) is obtained by the linear interpolation of the values of H (i, j) 
and H (i, j-1).  
To simplify exposition of the TCS method, we will use the shorthand notation 
x
?  
and 
y
? for spatial derivatives in x and y, resulting in:  
()
n 1 n n 1/2 n 1/2 n 1/2 n 1/2 2 n 1/2 2 n 1/2 n 1/2
xyxyx
t
UU U U V U U U g
2
++++++++
?
??
=+? ? ? ? +?? +? ???
??
(6.8) 
n 1 n n 1/2 n 1/2 n 1/2 n 1/2 2 n 1/2 2 n 1/2 n 1/2
xyxyy
t
VV U V V V V V g
2
++++ ++
?
=+? ? ? ? +?? +? ???  (6.9) 
Velocity, U Velocity, V
Surface elevation, ?  
j
1 2 3 4
1 
2 
3 
i 
U
V
j+1/2
i+1/2 
 86
 
()()
n 1 n n 1/2 n 1/2 n 1/2 n 1/2
xy
t
HU HV
2
++++
?
? ?
?=?+ ?? ??
? ?
 (6.10) 
Equations (6.8) through (6.10) are a nonlinear three-variable coupled time dependent 
system. However, this complex equation system can be computationally linearized and 
decoupled by applying the TCS method. The general form of TCS discretized SWE U-
velocity equation is: 
 
()()(){
()()()
()()
()()
()}
nn n
1x 1 x
nn
2y 2 y
22n
3x 3x
4y 4y
n
5x 5 x
t
UU UU 1 UU
2
VU 1 VU
U1 U
U1 U
g1g
***
**
*
*
*
?
=+??? +????
+? ? ? + ?? ? ?
+? ?? + ?? ??
+? ?? + ?? ??
????? ?? ??
 (6.11) 
 
 
()()(){
()()()
()
()
() }
n1 n1 n1
1x 1x
n1 n1
2y 2y
22n1
3x 3 x
4y 4 y
n1
5x 5 x
t
UU UU1 UU
2
VU 1 VU
U1 U
U1 U
g1g
***
**
*
*
*
++ +
++
+
+
+
?
=+??? +????
+? ? ? + ?? ? ?
+??? + ?? ??
+??? + ?? ??
?? ?? ? ?? ??
 (6.12) 
where the summation of Equations (6.11) and (6.12) is Equation (6.8). Similarly, for the 
other equations, it follows that, 
 
()()(){
()()()
()()
()()
()}
nn n
6x 6 x
nn
7y 7 y
22n
8x 8x
9y 9y
n
10 y 10 y
t
VV UV 1 UV
2
VV 1 VV
V1 V
V1 V
g1g
***
**
*
*
*
?
=+??? +????
+? ? ? + ?? ? ?
+? ?? + ?? ??
+? ?? + ?? ??
??????? ??
 (6.13) 
 
 87
 
()()(){
()()()
()
() }
() }
n1 n1 n1
6x 6x
n1 n1
7y 7y
22n1
8x 8 x
9y 9 y
n1
10 y 10 y
t
VV UV1 UV
2
VV 1 VV
V1 V
V1 V
g1g
***
**
*
*
*
++ +
++
+
+
+
?
=+??? +????
+? ? ? + ?? ? ?
+??? + ?? ??
+? ?? + ?? ??
?? ?? ? ?? ??
 (6.14) 
and the summation of Equations (6.13) and (6.14) is equivalent to Equation (6.9). The 
two steps of the ? -equation are: 
 
 
()()(){
()()()}
*n *n n*
11 x 11 x
*n n*
12 y 12 y
t
HU 1 HU
2
HV 1 HV
?
?=?? ?? + ?? ?
+? ? + ?? ?
 (6.15) 
 
()()(){
()()()}
n1 * * n1 n1 *
11 x 11 x
*n1 n1*
12 y 12 y
t
HU 1 H U
2
HV 1 H V
+++
++
?
?=?? ?? +???
+? ? + ?? ?
 (6.16) 
and the summation of Equation (6.15) and (6.16) is equivalent to Equation (6.10). In the 
above equation system, 
i
? is the weighting factor and { }i12 12,, ,? K for each term. As 
discussed in Chapter 5, properly chosen 
i
? can computationally decouple and linearize 
the equation system. In the next section, we will illustrate the possible decoupled 
discretizations for  Equations (6.1) through (6.4). 
 
6.3 DECOUPLING THE SWE 
Observing Equations (6.1) through (6.4), we discover that the advection terms 
VU/ y??and UV/ x?? couple the two momentum equations for U and V. g / x?? ? and 
HU / x?? couple the U-momentum equation and the continuity equation. g/y?? ? and 
HV / y?? couple the V-momentum equation and the continuity equation. Consequently, 
the weighting factors associated with these terms are crucial in the decoupling process. 
By examining Equations (6.11) through (6.16), the weighting factors are 
256101 12
and????? ?,,, , , . The other weighting factors, including the ones for the viscous 
terms,
34 89
?? ??,, , and the ones for non-coupled advection terms, 
17
??, , are not 
involved in the decoupling process. In other words, the values of the weighting factors, 
13478 9
and????? ?,,,,, , will not affect the decoupling process.  
There are limited combinations of the weighting factors, 
256101 12
and? ???? ?,,, , ,  
that will produce a completely decoupled equation system. Each combination will create 
 88
a unique discretization and a solution order. The principal for the decoupling process is 
that only one dependent variable can appear in each equation and we can change the 
number of the variables by changing the values of the weighting factors. This can be 
better understood by the following example:  
? If Equation (6.11) for 
*
U  is chosen to be solved first, 
2
1? = and 
5
0?=  are 
necessary conditions, because only one dependent variable 
*
U  can be 
included. Thus, the weighting factors for 
**
Vand?  have to be equal to 0. 
It follows that Equation (6.12) for 
n1
U
+
 will include 
n1 n1
Vand
++
? . To 
have a decoupled solution, Equation (6.12) for 
n1
U
+
has to be solved at the 
last. This unique solution sequence holds true for each variable. If 
*
V  or 
*
?  is chosen to be solved first, then 
n1
V
+
 or 
n1+
?  has to be solved last. 
? If Equation (6.13) for 
*
V  is solved right after 
*
U , 
610
0and 0? =?= are 
necessary conditions. In Equation (6.13), 
*
? cannot be included, so 
10
?  has 
to equal to 0, which causes 
n1+
?  to be included in  Equation (6.14) . To 
have a decoupled solution, Equation (6.16) for 
n1+
?  has to be solved 
before Equation (6.14) for 
n1
V
+
. Furthermore, in Equation (6.14), 
n1
U
+
 
cannot be included since  
n1
U
+
 will be obtained at the end, therefore 
6
?  
has to equal to 0. After choosing the second equation, two more weighting 
factors are decided and solution order is decided as 
*
U , 
*
V , 
*
?  then 
n1+
? , 
n1
V
+
 and 
n1
U
+
. 
? Since 
n1+
?  has to be solved before 
n1
U
+
 and 
n1
V
+
, 
11 12
0and 0? =?= are 
necessary conditions because 
n1
U
+
 and 
n1
V
+
cannot be included in 
Equation (6.16).  
? Hence, from the procedure described above, we can see that a solution 
order determines a unique combination of weighting factors. Since we 
have six equations in the system, six solution orders exist. Table 6.1 lists 
all six solution orders and their correspondent combinations of the 
weighting factors.   
 
 
 
 
 
 
 
 89
Table 6.1 Weighting factors and the solution orders for decoupled Equations (6.11) 
through (6.16). 
row 
2
?  
5
?  
6
?  
10
?
11
?
12
? solution order 
a 0 1 1 0 1 0 
*
V, 
*
? , 
*
U,  
n1
U
+
, 
n1+
? , 
n1
V
+
 
b 1 0 0 1 0 1 
*
U, 
*
? , 
*
V,  
n1
V
+
, 
n1+
? , 
n1
U
+
 
c 1 0 0 0 0 0 
*
U, 
*
V, 
*
? ,  
n1+
? , 
n1
V
+
, 
n1
U
+
 
d 0 0 1 0 0 0 
*
V, 
*
U,  
*
? ,  
n1+
? , 
n1
U
+
, 
n1
V
+
 
e 1 1 0 1 1 1 
*
? , 
*
U, 
*
V, 
n1
V
+
, 
n1
U
+
,  
n1+
?  
f 0 1 1 1 1 1 
*
? , 
*
V, 
*
U, 
n1
U
+
, 
n1
V
+
,  
n1+
?  
 
 
6.4 CHARACTERISTICS OF THE TCS DECOUPLED EQUATION SYSTEM 
The main advantage of the TCS method is that it can completely decouple and 
linearize an equation system with quadratic nonlinearity. As shown in the previous 
section, multiple TCS decoupled and linearized discretizations exist. To explore this 
characteristic of the TCS method, we choose values from row a, c and e in Table 6.1 to 
create three different TCS discretizations since the continuity equation is solved in 
different orders in these three rows. 
Choosing weighting factors from row a, Equations (6.11) to (6.16) become a two 
step computationally linearized and decoupled equation system. The implicit steps for 
solution of intermediate variables (
***
U,V,? ) are 
 
()
*n n*n* 2*2* n
xyxyx
t
VV UVVV V V g
2
?
??
=+????+??+? ???
??
 (6.17) 
 
()()
*n *n n*
xy
t
HU HV
2
?
? ?
?=?? ? +?
? ?
 (6.18) 
 
()
*n *n*n 2*2* *
xyxyx
t
UU UUVU U U g
2
?
=+????+??+? ??? (6.19) 
where 
**
b
ZH? ?=.  The implicit steps for the ?n+1? values are sequentially solved as 
 
()
n1 * * n1 * n1 2 * 2 * *
xyxyx
t
UU UUVU U Ug
2
+++
?
??
=+ ?? ?? +??+? ???
??
 (6.20) 
 
()()
n1 * * n1 n1 *
xy
t
HU H V
2
++
?
? ?
?=?? ? +?
? ?
 (6.21) 
 
()
n1 * n1 * n1 * 2 * 2 * n1
xyxyy
t
VV UVVV VVg
2
+++ +
?
=+ ??? ?+??+? ??? (6.22) 
 90
As we stated in section 6.3, weighting factors for the viscous terms and non-
coupled advection terms will not affect the decoupling process, so the values for these 
weighting factors can be anything from 0 to 1.  In Equations (6.17) through (6.22), the 
weighting factors for the viscous terms: 
3489
1? =? =? =? =  and the weighting factors 
for the non-coupled advection terms: 
1
1? =  and 
7
1? = . Thus the original set of the 
coupled nonlinear equations becomes a sequence of linear implicit equations for each 
variable. All variables in Equations (6.17) through (6.22) are solved in order 
as:
** * n1n1 n1
V, ,U ,U , ,andV
++ +
?? .  We call the solution from the above equation system 
?TCS solution 1?.  
In the same vein, if we choose the weighting factors at row c combining 
3489
1?=?=?=?=, 
1
1?=  and 
7
1? = , a different decoupled and linearized 
discretization is obtained. In the first step, intermediate variables (
***
U,V,? ) are solved 
as 
 
()
nnn 22 n
xyxyx
t
UU UUVU U U g
2
*****
?
??
=+????+??+? ???
??
 (6.23) 
 
()
nnn22n
xyxyy
t
VV UVVV V V g
2
*****
?
=+????+??+? ??? (6.24) 
 
 
()()
*n n* n*
xy
t
HU HV
2
?
? ?
?=?? ? +?
? ?
 (6.25) 
and in the second step, 
 
()()
n1 * n1 * n1 *
xy
t
HU HV
2
+++
?
? ?
?=?? ? +?
? ?
 (6.26) 
 
()
n1 n1 n1 2 2 n1
xyxyy
t
VV UVVV VVg
2
** * **++ +
?
??
=+ ?? ?? +??+? ???
??
 (6.27) 
 
 
()
n1 n1 n1 2 2 n1
xyxyx
t
UU UUVU UUg
2
*****+++ +
?
=+ ??? ?+??+? ??? (6.28) 
The solution order of the equation set (6.17) to (6.22) is changed to 
*** n1 n1 n1
U,V, , ,V ,U
+++
?? .  The solution from Equations (6.23) to (6.28) is called ?TCS 
solution 2?.  
Similarly, anther set of equations can be obtained by using the weighting factors 
in row e combining 
3489
1?=?=?=?=, 
1
1? =  and 
7
0? = . In the first step, 
 
()()
*n *n *n
xy
t
HU HV
2
?
? ?
?=?? ? +?
? ?
 (6.29) 
 
()
nnn 22
xyxyx
t
UU UUVU U U g
2
?
??
=+????+??+? ???
??
******
 (6.30) 
 91
 
()
nnn22
xyxyy
t
VV UVVV V V g
2
?
??
=+????+??+? ???
??
******
 (6.31) 
and in the second step, 
 
()
n1 n1 n1 2 2
xyxyy
t
VV UVVV VVg
2
+++
?
=+ ?? ?? +??+? ???
** * ***
 (6.32) 
 
()
n1 n1 n1 2 2
xyxyx
t
UU UUVU UUg
2
+++
?
??
=+ ??? ?+??+? ???
??
******
 (6.33) 
 
()()
n1 * * n1 * n1
xy
t
HU HV
2
++
?
? ?
?=?? ? +?
? ?
 (6.34) 
 
In the above equation system, the solution order is 
*** n1 n1n1
,U ,V ,V ,U ,
+ ++
? ? . We call 
this ?TCS solution 3?.  
Each TCS discretization has a different solution procedure. Table 6.2 provides the 
detailed solution procedures for each TCS method.   
 
Table 6.2 Solution procedures of the three TCS discretizations  
Step TCS solution 1 TCS solution 2 TCS solution 3 
*nnn
VU,V,? ?  
*nnn
UU,V,? ?  
*nnn
U,V,? ?? 
*n*n
U,V,? ?? 
**nn
VU,V,? ?  
*nn*
UU,V,?? 
1 
*n**
UU,V,?? 
***n
U,V,? ?? 
**nn
VU,V,?? 
n1 * * *
UU,V,
+
?? 
n1 * * *
U,V,
+
? ?? 
n1 * * *
VU,V,
+
?? 
n1 n1 * *
U,V,
++
? ?? 
n1 * * n1
VU,V,
+ +
? ?  
n1 * n1 *
UU,V,
++
?? 
2 
n1 n1 * n1
VU,V,
++ +
?? 
n1 * n1 n1
UU,V,
+ ++
? ?  
n1 n1 n1 *
U,V,
+++
? ?? 
 
6.5 NUMERICAL TESTS  
To test the general performance of the TCS method in solving 2D depth averaged 
SWE, a series of numerical experiments are presented in this section. In these numerical 
experiments, we want to test 1) how the method treats the nonlinear terms such as 
UV/ x??, VU/ y??, HU / x?? and HV / y? ? ; 2) how the method decouples the terms 
U, V, and H; and 3) How different discretizations and  different solution procedures 
affect the solutions.  
A one-dimensional standing wave in a closed basin is first tested to compare our 
numerical model with the analytical solution. A two-dimensional standing wave in a 
closed basin is simulated to investigate the different characteristics among three TCS 
 92
solutions. As an initial attempt to simulating an open boundary case, a one-dimensional 
wave traveling through an open boundary system is tested in Appendix A.     
 
6.5.1 One dimensional Standing Waves 
As a first test of the TCS method, a free surface standing wave in a rectangular 
basin has been simulated and is shown schematically in Figure 6.2. Free-slip boundary 
conditions are enforced on all walls. An equally-spaced 20 5?  mesh is applied to the 
computational domain.   
L
H
2a
x
yz
 
Figure 6.2 A standing wave in a rectangular basin 
 
Analytical Solution 
The period of an inviscid free surface wave is (Dean and Dalrymple 1998): 
 ()
12
TtanhkH
g
???
= ? ?
? ?
 (6.35) 
According to linear wave theory, the wave function is a sinusoid for small amplitude 
waves. The surface height (h) above the still water level is: 
 ( ) ( )h(x,t) asin kx sin t= ?  (6.36) 
where the frequency (?) is 2/T? .  For a viscous unconfined wave in deep water, the 
evolution of the wave amplitude can be approximated as (Lamb 1932), 
 
2
2kt
a(t) a(0)e
??
=  (6.37) 
 93
where ?a? is the free surface wave amplitude, ??? is kinematic viscosity and ? k ? is the 
wavenumber.  However, wave damping based on Equation (6.37) is unlikely to be well-
represented in a 2D depth-averaged model as vertical velocity and vertical velocity shears 
are part of the closure term. The Reynolds number is defined as: 
 
Lu
Re =
?
 (6.38) 
where L  is the length of the basin, ? is the kinematic viscosity and u  is the 
characteristic Cartesian fluid velocity based on the wave amplitude (a) and the wave 
frequency (?) : 
 ua= ? (6.39) 
 
Results and Discussion 
Simulations have been run for different cases to examine a variety of model 
characteristics. Three simulations with different initial wave steepness (i.e. ratio of wave 
amplitude to basin length) are presented in Figure 6.3. The three simulations were run in 
an inviscid flow to minimize the viscous damping effects. In Figure 6.3, results of the 
simulations illustrate the nondimensionalized water surface elevation (i.e.
()H/a?? ) 
at xx/2=? , which is the grid cell center closest to the left wall. The important 
parameters for these three test cases are listed in Table 6.3.   
 
Table 6.3 Parameters of simulations of a 1D standing wave 
case H/L a/L L/? 
6.4a 0.05 0.05% 0.5 
6.4b 0.05 0.2% 0.5 
6.4c 0.05 0.5% 0.5 
 
 
 
 
 94
0 1 2 3 4 5 6 7
?1
0
1
t/T
(
?
?H)/a
case a: a/L=0.05%
0 1 2 3 4 5 6 7
?1
0
1
t/T
(
?
?H)/a
case b: a/L=0.2%
0 1 2 3 4 5 6 7
?1
0
1
t/T
(
?
?H)/a
case c: a/L=0.5%
 
 
TCS solution 1 TCS solution 2 TCS solution 3 analytical solution
 
Figure 6.3 Simulations in the inviscid flow with different wave steepness 
 
Figure 6.3 shows results of the three different initial wave steepness. The non-
linearity of the advection term causes nonlinear wave steepening. Hydrostatic 
approximation from the depth averaged SWE enhances the nonlinear steepening because 
the dispersive effects from the non-hydrostatic effects are neglected (Daily and Imberger 
2003). This nonlinear steepening can be observed from Figure 6.3. Furthermore, with the 
increase of the nonlinearity of the wave (i.e. the initial wave steepness), the nonlinear-
steepening effects increase. In case 6.3a, a very small initial wave steepness is introduced 
and the simulation shows little nonlinear effects. When the initial steepness is increased 
in case 6.3b and 6.3c, the nonlinear effects become apparent. In theory, a series of 
solitary wave can be formed when the nonlinear steepening is balanced by the non-
hydrostatic pressure gradient (Miropol'sky 2001). However, In Figure 6.3c, a train of 
solitary waves can be observed from the TCS simulations, which are solutions under 
hydrostatic approximation. This is probably due to the numerical dispersion, which plays 
an opposition effect to the nonlinear steepening. This phenomena has been reported in 
other literatures (Hodges et al. 2006; Wadzuk 2004).  
Figure 6.4 displays the viscous damping effects in our shallow water model. The 
Reynolds number is equal to 0.1, which indicates a very viscous flow. The simulation 
results damped faster than the analytical solution. This is mainly a result of the staggered 
grids used in the numerical model. The staggered grids necessitate averaging from the 
cell surfaces, which creates numerical dissipation (Garcia and Kahawita 1986).   
 
 95
0 1 2 3
?1
0
1
t/T
(
?
?H)/a
Re=0.1, H/L=0.05, a/L=0.05%
 
 
TCS solution 1 TCS solution 2 TCS solution3 analytical solution
 
Figure 6.4 Viscous damping effects in the shallow water model 
 
In Figures 6.3 and 6.4, the results from the three TCS discretizations are 
indistinguishable. To have a better comparison among three TCS solutions, we run a 
simulation of a two dimensional standing wave in a square box in the next section.  
 
6.5.2 Two-dimensional Standing Waves 
In this test case, the ability of the TCS method to simulate a two dimensional 
finite-amplitude free surface wave in a square box is investigated and the performances 
from different TCS discretizations are compared.  
A two-dimensional free surface standing wave in a square box is shown in Figure 
6.5. This initial wave is linearly composed by two one-dimensional (one in x-direction 
and the other one in y direction) identical and orthogonal standing waves, which is also 
shown in Figure 6.5. This superimposed condition creates a two dimensional standing 
wave that is exactly symmetric to the diagonal plane of the square box.  The horizontal 
scale of the square box is 10m 10? m, and the initial still water elevation is 0.5m. The 
initial wave amplitude for the two one-dimensional waves is 0.005m. The computational 
domain is illustrated in Figure 6.6.  The surface elevation is measured from the still water 
surface and nondimensionalized by the wave amplitude which is the summation of the 
one dimensional wave amplitudes. The horizontal length scales are normalized by the 
wave length. A 20 20? mesh is used in this simulation and the free slip condition is 
enforced on all the side walls. The simulation is run in an inviscid flow condition to 
eliminate the viscous effects.  Table 6.4 listed all the important parameters in the 
simulations of this 2D standing wave. 
 
 
 96
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
?1.5
?1
?0.5
0
0.5
1
1.5
y/?
initial suface elevation
x/?
(
?
?H)/2a
(0, 0)
(0, 0.5)
(0.5, 0.5)
(0.5, 0)
 
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
?1
0
1
x/?
initial x?wave
y/?
(
?
?H)/a
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
?1
0
1
x/?
initial y?wave
y/?
(
?
?H)/a
 
Figure 6.5 The initial 2D standing wave and its decomposed x and y direction 1D wave 
 97
10m
10m
0.5m
initial wave amplitude: 0.01m
 
Figure 6.6 Computational domain (not to scale) 
 
Table 6.4 Parameters of simulations of the 2D standing wave 
L ? a H 
x?  
y?  
t?  ? 
10m 20m 0.005m 0.5m 0.5m 0.5m 0.01s 0 
 
 
Although six representative decoupled linearized TCS discretizations and their 
correspondent solution orders are listed in Table 6.1, only three discretizations need to be 
tested here because velocity U and V are given symmetrically. In TCS solution 1, the 
surface elevation 
*
?  and 
n1+
?  are solved in between two velocity equations; In TCS 
solution 2, the surface elevation  
*
?  is solved after the two velocity equations and 
n1+
?  is 
solved before the two velocity equations; In TCS solution 3, 
*
?  is solved before the two 
velocity equations and 
n1+
?  is solved after the two velocity equations.  
To explore how different discretizations and solution orders affect the results. We 
first compare the results from TCS solutions 1, 2 and 3 at different simulation times. The 
simulation time is normalized by the wave period, T.  In Figure 6.7, the surface elevation 
contour calculated using these three TCS forms are plotted for the whole computational 
domain at time 4T. The surface elevation is measured from the still water surface and 
normalized by the amplitude. The horizontal length scale is normalized by the wave 
length, ?. The surface elevation contour from TCS solution 2 and 3 are nearly identical, 
and the result from TCS solution 1 behaves differently than that from solution 2 and 3. 
 98
The result from each solution has a symmetric shape. When the simulation time is 
increased to 15T (results shown in Figure 6.8), TCS solution 2 and 3 still perform 
similarly, but TCS solution 1 has an observable difference. Similar phenomena can be 
better observed in the velocity U and V contours. Figure 6.9 and 6.10 demonstrate the U-
velocity simulated using three TCS methods at 4T and 15 T respectively. Figure 6.11 and 
6.12 display the V-velocity at 4T and 15T respectively. Each velocity component in the 
figures is normalized by its maximum value. The results form TCS solution 2 and 3 are 
almost the same, but the result from TCS solution is different. 
This can be analyzed by the solution procedures of the three tested TCS 
discretizations presented in Table 6.2. First of all, in TCS solution 1, the continuity 
equation (solving ? ) is solved in between two momentum equations (solving U and V); 
but in both TCS solutions 2 and 3, the two momentum equations are solved 
consecutively. The only difference between TCS solution 2 and 3 is that one solve 
momentum equations first, the other solve continuity equation first in one time step. 
However, this difference between solution 2 and 3 could be reduced with the marching of 
the simulation time.  Second, in TCS solution 1, the surface elevation 
*
? is solved using 
the velocity components, 
n
U  and 
*
V  in the first step; and 
n1+
?  is solved using the 
velocity components, 
n1
U
+
 and 
*
V  in the second step.  Thus, in TCS solution 1, the 
surface elevation is obtained using a velocity vector field which is constructed by the 
components from different time steps. Third, we further discover that water depth H is 
not used symmetrically in the continuity equation. In Equation (6.18), 
n
H  is used in the x 
direction, but 
*
H  is used in the y direction to solve 
*
?  in the first step; in Equation (6.21), 
*
H  is used in x direction but 
n1
H
+
 is used in the y direction to solve 
n1+
? .  All of these 
asymmetric solution procedures introduce extra numerical errors since our initial 
condition is symmetric. However, in both TCS solutions 2 and 3, the surface elevation is 
obtained using the velocity components and water depth at the same time step although 
the velocity field and water depth are not exactly symmetric because of the sequential 
solution procedure.  The difference caused by different solution procedure is further 
exemplified by a cross mode analysis in the following.  
 With the given symmetric initial condition, ideally, we should expect this 
standing wave oscillates along one diagonal line and no wave motion should be expected 
in the cross direction. That means, the wave surface oscillates at point (0, 0) and (0.5, 0.5) 
with an amplitude 0.01m ( 2 0.005m? ), and the wave surface should be fixed at point 
(0.5, 0) and (0, 0.5). However, the TCS decoupled equation system is solved in a 
sequential order, which causes each equation to be solved anisotropically. This in turn 
induces numerical error that causes cross mode wave motion. The similar phenomena 
have been observed in Alternate Direction Implicit (ADI) method. In ADI method, this 
cross mode is induced by solving the equation in different coordinate direction as 
observed by Hodges (1997). To examine how cross mode is established in the three TCS 
discretizations, we plot the surface displacement at point (0.5, 0) and (0, 0.5) against 
simulation time in Figure 6.13. The surface displacement is normalized by the wave 
amplitude and the simulation time is normalized by wave period. The results from TCS 
solutions 2 and 3 are indistinguishable. The cross mode at the two monitored points are 
 99
nearly identical in these two TCS solutions and gradually increase with time. At time 15 
T, the surface displacement at those two points is almost 25% of the wave amplitude. In 
contrast, this numerical cross mode from solution 1 behaves rather differently. The 
surface at the two monitored points oscillates with an amplitude and a period the same as 
the initial standing wave. In addition, the oscillation of the cross direction is one half 
phase behind the initial standing wave. To have a better visualization of the cross mode, 
Figure 6.14 and Figure 6.15 display the snapshots of the standing wave at 0.5 T and 15T 
respectively. At 0.5 T, no obvious motion can be observed at the two monitored points 
solved from TCS solutions 2 and 3. With the increase of the simulation time, at time 15 
T, an observable surface displacement can be discovered from both TCS solution 2 and 3. 
In contrast, a maximum oscillation with a displacement of 2a can be observed for TCS 
solution 1 at 0.5 T, but no obvious motion can be found at time 15T. This is probably 
because in TCS solution 1, the surface elevation is solved using asymmetric velocity 
components and water depth components. There is a half time step difference between x 
and y direction in both velocity and water depth components, which causes the wave at 
the cross direction oscillates with a half phase lag.  
 
 
 
 100
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
4 T
 
Figure 6.7 Normalized surface elevation contours in the computational domain at time 
4T.  
 101
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
15 T
 
Figure 6.8 Normalized surface elevation contours in the computational domain at time 
15T. 
 
 102
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
4 T
 
Figure 6.9 Normalized U-velocity contours in the computational domain at time 4T.  
 
 
 
 103
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
15 T
 
Figure 6.10 Normalized U-velocity contours in the computational domain at time 15T.  
 
 104
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
4 T
 
Figure 6.11 Normalized V-velocity contours in the computational domain at time 4T.  
 
 105
x/?
y/
?
TCS solution 1
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 2
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
x/?
y/
?
TCS solution 3
 
 
0 0.5
0
0.5
-1
-0.5
0
0.5
1
15 T
 
Figure 6.12 Normalized V-velocity contours in the computational domain at time 15T.  
 
 
 
 
 
 106
0 1 2 3 4 5 6 7 7.5
?1.5
?1
?0.5
0
0.5
1
1.5
t/T
(
?
?H)/2a
7.5 8 9 10 11 12 13 14 15
?1.5
?1
?0.5
0
0.5
1
1.5
t/T
(
?
?H)/2a
 
 
 
 TCS solution 1 at point (0.5, 0)    TCS solution 1 at point (0, 0.5)   
 TCS solution 2 at point (0.5, 0)     TCS solution 2 at point (0, 0.5) 
 TCS solution 1 at point (0.5, 0)     TCS solution 3 at point (0, 0.5) 
 
Figure 6.13 Surface displacement at points (0.5, 0) and (0, 0.5) simulated using three 
TCS solutions 
 
 
 
 107
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 1
y/?
(
?
-H
)/
2
a
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 2
y/?
(
?
-H
)/2
a
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 3
y/?
(
?
-H
)/
2
a
0.5 T
 
Figure 6.14 Snapshot of the standing wave at 0.5 T. The surface displacements at the two 
monitored points (0.5, 0) and (0, 0.5) are circled.  
 108
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 1
y/?
(
?
-H
)/
2
a
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 2
y/?
(
?
-H
)/
2
a
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
-1
0
1
x/?
TCS solution 3
y/?
(
?
-H
)/
2
a
15 T
  
 
Figure 6.15 Snapshot of the standing wave at 15 T. The surface displacements at the two 
monitored points (0.5, 0) and (0, 0.5) are circled.  
 
 
 
 109
6.6 SUMMARY  
In this chapter, the TCS method is applied to a multi-variable and multi-
dimension equation system: the depth averaged SWE. The TCS method not only 
computationally discretizes the nonlinear term in the equation set, but also decouples the 
variables in the SWE. Three TCS discretizations and are derived for the SWE. Each TCS 
decoupled and linearized equation system has a unique solution order. 
A one-dimensional and two-dimensional standing waves in a closed rectangular 
domain are simulated using all three TCS discretizations. In the one dimensional standing 
wave case, we compare our numerical results with the analytical solutions. The results 
demonstrate that our TCS method well capture the nonlinear effects and display the 
combination effects of nonlinear steepening and hydrostatic approximation. In the 
viscous flow test case, the staggered grids used in the simulation causes numerical 
dissipation. In the two-dimensional standing wave case, we compare the performance 
from three representative TCS solutions: TCS solutions 1, 2 and 3. TCS solutions 2 and 3 
are solved similarly because the surface elevation is obtained using the velocity 
components at the same time step. However, TCS solution 1 performs differently because 
the surface elevation is obtained using the velocity components from different time step.  
The anisotropically solution process also causes numerical cross mode. The cross 
mode analysis shows that TCS solution 2 and 3 has smaller numerical cross mode 
compare to TCS solution 1. However, the cross mode gradually increase with simulation 
time. The solution procedure in TCS solution 1 causes wave oscillate in the cross 
direction. Therefore, when we apply TCS method to decouple the equation system, we 
have to be aware of the numerical errors generated by the decoupling process.   
  
 110
Chapter 7 Conclusions and Recommendations 
 
The main objective for this research is: 
Develop a new implicit solution for unsteady nonlinear advection 
problems. The new numerical algorithm should have the advantage in both 
stability and efficiency while keep the 2
nd
-order temporal accuracy as in 
the most existing shallow water equations (SWE) models. In addition, this 
new method would provide a novel approach in decoupling a coupled 
equation system. 
 
This objective has been accomplished. A new finite difference temporal scheme 
(Time-centered split method) is developed to solve the unsteady nonlinear advection 
problems. This new method provides a general approach for solving any unsteady 
quadratic nonlinearity. In this chapter, we summarize 1) discussions 2) conclusions 3) 
recommendations for future work. 
 
7.1 SUMMARY OF DISCUSSION 
A new computational linearization method, time-centered split (TCS) method, is 
derived from the midpoint rule temporal discretization. The fundamental principle of the 
TCS method is that it splits the quadratic nonlinear term into two steps so that each step 
is computationally linear. The two-step equation system is 2
nd
-oder equivalent to the 
original midpoint rule discretization. The essential theoretical development of the TCS 
method is illustrated in Chapter 3. The derivations in Chapter 3 demonstrate one of the 
most important significances of the TCS method: it is a direct linearization method while 
in an implicit format. The conventional implicit nonlinear solutions including Newton 
method, Picard method and local linearization method require either 1) outer iteration ( 
Newton and Picard methods) or 2) calculations of the Jacobian (Newton and local 
linearization methods). Therefore, compared to the conventional implicit nonlinear 
solutions, the TCS method has the advantage in efficiency and stability.  
In a time marching differential equation, the time-centered split concept can be 
applied to different terms: either part of the quadratic nonlinear term and the linear term. 
Different TCS discretizations can be generated when we apply the split to different terms. 
Based on this characteristic of the TCS method, a general form of the TCS family method 
is created in Chapter 5. A weighting factor 
i
?  (
i
01??? ) is introduced in the TCS 
general form. The weighting factor 
i
?  combines different approximations (splitting on 
different terms) of the nonlinear terms and the linear terms. The value of 
i
?  will not 
affect the 2
nd
-order temporal accuracy of the TCS method. The theoretical proof is shown 
in Chapter 5. This revealed another important significance of the TCS method: unlimited 
TCS discretized forms can be created for one problem. Thus, the TCS method provides 
 111
flexibility when we solve a specific problem because we can choose the discretizations 
based on the specific requirement of the problem.    
The significance of the weighting factors is further exemplified when the TCS 
method is applied to a multi-variable and multi-dimension equation system. Properly 
chosen weighting factors can not only linearize but also decouple the equation system 
without outer iterations. This property theoretically enhances the efficiency advantage of 
the TCS method. The 2D Burgers? equation is used as the example to show the process 
and principle of the decoupling in Chapter 5. The decoupling characteristics of the TCS 
method is analyzed in detail using 2D depth averaged shallow water equations (SWE) in 
Chapter 6. Both cases show that the TCS method can computationally linearize and 
decouple an equation system without outer iterations. Each variable in the TCS linearized 
and decoupled two-step system is solved in a sequential order.  
To verify and explore the basic characteristics of the TCS method, a 1D problem 
is first chosen as the test case. In Chapter 4, three 1D differential equations: 1D 
conservative Burgers? equation, 1D non-conservative Burgers? equation and 1D nonlinear 
ordinary differential equation (ODE) are solved using the TCS method. In all the three 
equations, we apply the time-centered split to different terms: either flux or gradient part 
in the nonlinear advection term in the Burgers? equation, the quadratic nonlinear term of 
the ODE, the diffusion term in the Burgers? equation and the linear term in the ODE. 
Consequently, multiple TCS disretizations are created such as TCSF (split the flux term), 
TCSG (split the gradient term), TCSF-D (split flux and diffusion/linear term) and TCSG-
D (split the gradient and the diffusion/linear term). To compare to the conventional 
implicit linearization method, all three equations are also solved using the implicit Crank-
Nicolson scheme combining with the conventional computational linearization methods: 
Newton method, Picard method and local linearization method. In the two Burgers? 
equation cases, each TCS discretization demonstrates 2
nd
-order temporal accuracy and 
remains stable up to a CFL O (10).   The stability advantage of the TCS method is further 
demonstrated by comparing it to the Runge-Kutta (RK) method. The efficiency of the 
TCS method is examined using operation count at each grid point at each time step. 
Results in the two 1D Burgers? equations show that the TCS method requires fewer 
operations than Newton and Picard methods but has a similar computation expense to the 
local linearization method. The principal advantage of the TCS method over local 
linearization is the relative ease with which the TCS method can be derived and 
implemented as it does not require discrete evaluation of a function Jacobian. The test 
case in 1D nonlinear ODE verified that the TCS method can be applied to any quadratic 
nonlinearity.  
In the 1D ODE test case, two different TCS discretizations collapse into one 
expression, which is also mathematically equivalent to the local linearization. However, 
in both conservative and non-conservative Burgers? equation cases, different TCS 
disretizations have different relative errors although they perform similarly in temporal 
accuracy, stability and efficiency. This difference can also be observed in Chapter 5. 
Because of the introduction of the weighting factors, more TCS discretizations are tested 
for the 1D non-conservative Burgers? equation in Chapter 5. Our results show that 
changing 
i
?  will not affect the overall 2
nd
-order temporal accuracy of all the TCS 
 112
discretizations. However, the relative accuracy of each TCS discretization is different. 
The results show that the most accurate solution occurs when the two weighting factors 
equal to some value between 0 and 1 but not 0 or 1. The similar results can also be 
observed for the stability. The stability is also enhanced when the weighting factors are 
not equal to 0 or1. This can be explained by mathematical meaning of the weighting 
factors. The weighting factors control how we approximate each term in the equation. 
When the weighting factors equal to 0 or1, we only use one kind of approximation: either 
approximate 
n1/2
u
+
 using the time-centered split completely or approximate 
n1/2
u
+
 using 
*
u  completely. Our results suggest that combing both approximations will enhance the 
performance of the TCS method in both accuracy and stability. However, the best 
combinations of the weighting factors for accuracy and stability is problem-specific and 
one needs to analyze each approximated term and the nature of the solution itself at the 
designated marching time. This phenomenon verified that TCS method provides 
flexibility in solving a specific problem because a more accurate or stable solution can be 
obtained by changing the choice of the weighting factors. 
In Chapter 6, the TCS method is used to solve 2D depth averaged SWE. We 
discover that six terms in the SWE couple the entire equation system. These six terms are   
VU/ y??, UV/ x??, HU / x??, HV / y? ? , / x??? and / y???. Their associated 
weighting factors are crucial in the decoupling process. It is revealed that only six 
combinations of these weighting factors can fully decouple the SWE. Each combination 
of the weighting fact has a unique solution order. Three representative TCS solutions are 
chosen to solve the SWE based on the solution order of the continuity equation. In the 
one-dimensional standing wave case, how the TCS method treat the nonlinear term is 
tested. The results show that that the TCS methods well capture the nonlinear effects and 
display the combination effects of numerical dispersion and hydrostatic approximation. 
How the solution order affects the results is analyzed in the two-dimensional standing 
wave case. The results show that TCS solution 1 perform differently than TCS solutions 
2 and 3 because the surface elevation is obtained using the velocity components from 
different time step in TCS solution 1. The fully-decoupled TCS-discretized SWE is 
solved in a sequential order, which causes each equation solved anisotropically. This in 
turn causes the numerical cross modes. The cross mode from TCS solution 2 and 3 
behave similarly: the cross mode gradually increases with time. However, the solution 
procedure in TCS solution 1 causes wave oscillate in the cross direction with the same 
amplitude as the initial standing wave and half phase behind the initial standing wave. 
Hence, we need to take into account the numerical errors generated by the decoupling 
process when we apply the TCS method to decouple an equation system.  
   
 
 
 
 113
7.2 CONCLUSIONS 
The conclusions that may be drawn from the present work are: 
? A new time marching numerical algorithm for solving unsteady nonlinear 
advection problems is proposed using a time-centered split (TCS) 
technique. 
? The principle advantage of the TCS method is it computationally 
linearizes the implicit nonlinear advection without either 1) iterations or 2) 
calculations of the Jacobian.  
? The TCS method is 2
nd
-order temporal accurate and has advantages in 
stability and efficiency. 
? A family of the TCS discretizations can be generated using the same 
principle. This provides flexibility when solving a specific problem using 
the TCS method. 
? The TCS method can fully decouple an equation system. However, 
additional numerical error is introduced through the decoupling process.  
 
 7.3 RECOMMENDATIONS FOR FUTURE WORK 
The TCS method is proposed and verified using 1D transport equation, 1D 
nonlinear ODE, 2D transport equation. The following potential application is 
recommended for future work.   
Application of  the TCS method to a coupled nonlinear ODE system 
This research is motivated by solving unsteady nonlinear advection problems. 1D 
and 2D flow transport problems are tested using the TCS method. However, the new TCS 
method provides a general approach to solve any quadratic nonlinearity. Although a 1D 
nonlinear ODE is used as an example in Chapter 3, application of the TCS method to the 
nonlinear ODE can be further investigated. A coupled quadratic nonlinear ODE system 
can be used as an example to verify the capability of the TCS method in computational 
linearization and decoupling. Moreover, various TCS discretizations can be compared 
and tested by changing the weighting factors.   
 
More test cases for the SWE 
In Chapter 6, initial tests of 2D depth averaged SWE have been conducted. All the 
test cases in 2D standing wave are in an inviscid flow condition. Since the weighting 
factors associated with the viscous terms will not affect the decoupling process, the 
inviscid flow condition is a valid test condition for investigating decoupling process. 
However, the weighting factor associated with the non-coupled advection terms such as 
UU/ x?? and VV/ y??can be any value from 0 to 1. In Chapter 6, we didn?t change the 
weighting factors for these two terms to see how these two will affect the solution. Based 
on the analysis in Chapter 5, the numerical error induced by the decoupling might be 
decreased by choosing different weighting factors for these two terms. An additional test 
for numerical error could be conducted by decreasing the time interval t? .  
Solve problems with discontinued flow conditions 
 114
We often need to solve a transport problem with a discontinuity flow condition, 
for example, hydraulic jump. All the test cases conducted in this research use central 
difference spatial discretizations, which cannot represent the sharp front of a shock wave. 
Test cases can be designed by the TCS temporal discretization combining upwind spatial 
discretization to simulate a shock wave.  
 115
Appendix A Test of a Progressive Wave in an Open Boundary System 
 
Another way to verify our numerical method is to simulate a wave progressing in 
an open boundary system. This test case simulates a free surface oscillation in a 
rectangular channel, which is open at both upstream (inlet) and downstream (outlet 
boundary). This also could be thought of as a step leading to the simulation of a river 
flow. In this section, a progressive wave is introduced by oscillating the inlet boundary of 
an open quiescent water body. A sponge layer is successfully applied as the open 
boundary condition at the outlet of the domain. Our numerical models satisfactorily 
simulate a progressive wave traveling through an open boundary system. As shown in 
Chapter 6, TCS solutions 2 and 3 give similar results. Thus, in this test, only results from 
TCS solution 1 and 2 are presented and compared.  
 
A.1 INITIAL AND INLET BOUNDARY CONDITION 
Initially, the water is at rest with uniform depth: 
  ( )
0
Hx,0 h=  (A.1)   
The inlet water level suddenly raised at t=0 and oscillated as a cosine wave described as 
 
0
H(0,t) h acos( t)= +? (A.2) 
where a is the progressive wave amplitude and ?is the frequency. The initial condition 
and inlet boundary condition is shown in Figure A.1 and A.2. 
 
Figure A.1 Initial water level, H (x, 0), in the rectangular open channel  
 
 
0
h 
x 
H 
 116
 
Figure A.2 Inlet boundary condition H (0, t) 
 
A.2 OUTLET BOUNDARY CONDITION 
One of the challenges to simulate an open boundary system with a finite 
computational domain is that open boundary conditions are required at both the inlet and 
outlet. Developing outlet non-reflective boundary conditions has been the subject of 
extensive research (Blayo and Debreu 2005; Tsynkov 1998).   
In the present work we apply an artificial damping layer (or sponge layer) 
upstream of the outlet boundary and downstream of the ?test section? (i.e. the 
computational domain of interest).  In this approach, an artificial damping function is 
prescribed over a range of grid cells to dissipate the surface wave and its reflections 
before it propagates back into the test section (Durran 1999). Although a sponge layer has 
additional computational costs associated with computations outside of the test section, 
the ratio of the additional cost to the initial cost is generally small (Blayo and Debreu 
2005).  In the present work, the sponge layer damping function uses an increasing 
viscosity from the end of the test section to the outlet (Vinayan 2003). We use an inviscid 
progressive wave for the test case, so the viscosity is set to be zero within test section.  In 
general, the viscosity is represented as a function of position such that
max
0(x)?? ?? .  
If we define the viscosity is a function of x, 
 
max
(x) f(x)? =?  (A.3) 
where, f(x)is a function satisfying the following conditions: 
 
r
s
r
s
f(x ) 0
f(x ) 1
df (x )
0
dx
df (x )
0
dx
=
=
=
=
 (A.4) 
 117
where x
r
 and x
s
 are the upstream and downstream limits of the sponge layer, as shown in 
Figure A.3.  If the artificial viscosity changes rapidly, spurious reflections may also 
occur. A gradual change, as shown in Figure A.3, can be invoked by defining f (x) as   
 
r
sr
xx1
f(x) 1 cos
2xx
?
? ?
???
=??
? ???
??
? ?
 (A.5) 
 
 
Figure A.3 Shape function of viscosity 
 
Using the above definitions, a sponge layer is defined by its length (
sr
xx? ) and its 
maximum viscosity
max
? . In this work, 
sr
xx2? =? and 
max
? is basically a relative big 
number obtained by trial and error. The sponge layer arrangement is shown in Figure 
A.4.  
 
 
Figure A.4 An open boundary rectangular channel with a sponge layer 
 
These initial and boundary conditions generate a traveling wave with height 2a, 
progressing along the channel. To see the damping effects of the sponge layer, we 
monitored the surface height of point A, which is the point in the test section closest to 
 118
the sponge layer, and B, which is the point within the sponge layer closest to the outlet, as 
shown in Figure A.5.   
 
 
Figure A.5 Schematic illustrations of grids 
 
 
0 2 4 6 8 10 12 14 16 18 20
-2
-1
0
1
2
(
?
-H
)/
a
 
a
t
 p
o
i
n
t A
 w
i
t
h
/wi
t
h
o
u
t
 t
h
e
 s
p
o
n
g
e
 l
a
y
e
r
t/T
 
with sponge layer without sponge layer 
Figure A.6 Water level at point A with/without the sponge layer 
 
Figure A.6 shows the evolution of the nondimensionalized water depth (()H/a?? ) at 
point A with and without a sponge layer. Without the sponge layer, when the wave 
reaches the solid wall where zero flux is prescribed, reflection occurs and its amplitude is 
increased to about two times of the original one. However, with a sponge layer condition, 
the wave amplitude stays the same. We also compare the wave at point A and B in Figure 
A.7. The wave inside the sponge layer is damped about 10% of the wave in the test 
section. Thus the reflection effect is minimized.  
 
AB
x
z
test section sponge layer
 119
0 2 4 6 8 10 12 14 16 18 20
-1
0
1
t/T
(
?
-H
)/
a
 
a
t
 p
o
i
n
t A
 
a
n
d
 B
 
 water level at point A, inside the test section 
 water level at point B, inside the sponge layer 
Figure A.7 Water level evolutions inside the test section and sponge layer. 
 
A.3 RESULTS AND DISCUSSION 
Figure A.8 compares the free surface shape in the entire domain (including test 
section and sponge layer) at 7T and 14T for TCS solution 1 and 2 respectively. At 7T, the 
progressive wave has traveled through the domain. When we double the simulation time 
to 14T, the free surface shape remains almost the same. This indicates that we 
successfully simulated a wave progress through an open boundary system, and our 
numerical results model the wave characteristics such as wave speed and period well. 
Both TCS solutions 1 and 2 give us similar results.  
To further compare the two TCS solutions  in this progressive wave test case, we 
compare the free surface shape calculated using both TCS discretizations in the entire 
domain at 7T and 14 T, respectively, in Figure A.9. The results simulated using each TCS 
discretization overlapped at 7T. After a longer simulation time, this similarity between 
two solutions still holds true at 14T.  
The time evolution of free surface at different periods in the entire domain for 
TCS solutions 1 and 2 are shown in Figure A.10 and Figure A.11 respectively. 
 
 
 
 120
0 1
?1
0
1
L/L
(
?
?H)/a
TCS solution 1
 
 
0 1
?1
0
1
L/L
(
?
?H)/a
TCS solution 2
 
 
7T
14T
7T
14T
 
Figure A.8 Comparison of wave shape at time 7T and 14T for TCS solution1 and TCS 
solution 2 respectively 
 
 121
0 1
?1
0
1
L/L
(
?
?H)/a
7T
 
 
TCS solution 1
TCS solution 2
0 1
?1
0
1
L/L
(
?
?H)/a
14T
 
 
TCS solution 1
TCS solution 2
 
Figure A.9 Comparison of wave shape simulated from TCS solution 1 and 2 at time 7T 
and 14T respectively 
 
 
 
 122
1T 2T
4T 6T
7T 14T
 
Figure A.10 Time evolution of a progressive wave at different wave periods from TCS 
solution 1. 
 
 
 
 
 
 
 
 123
1T 2T
4T 6T
7T 14T
 
 
Figure A.11 Time evolution of a progressive wave at different wave periods from TCS 
solution 2. 
 124
References 
Agoshkov, V., Ovchinnikov, E., Quarteroni, A., and Saleri, F. (1994). "Recent 
Developments in the Numerical Simulation of Shallow Water Equations. II. Temporal 
Discretization." Mathematical Models and Methods in Applied Sciences, 4(4), 533-566. 
Amat, S., Busquier, S., and Gutierrez, J. M. (2003). "Geometric constructions of 
iterative functions to solve nonlinear equations." Journal of Computational and Applied 
Mathematics, 157(1), 197-205. 
Anderson, D. A., Tannehill, J. C. and Pletcher, R. H. (1997). Computational Fluid 
Mechanics and Heat Transfer. 2
nd
 ed. Hemisphere Publishing, MacGraw-Hill. 
Andrus, C. W., Long, B. A., and Froehlich. H. A. (1988). ?Woody debris and its 
contribution to pool formation in a coastal stream 50 years after logging.? Can. J. of 
Fisheries and Aquatic Sci., 45, 2080-2086. 
Ancey, C., Iverson, R. M., Rentschler, M., and Denlinger, R. P. (2008). "An exact 
solution for ideal dam-break floods on steep slopes." Water Resources Research, 44(1). 
Arakawa, A., and Lamb, V.R. (1977). ?Computational design of the basic 
dynamical processes of the UCLA general circulation model.? Methods in Computational 
Physics, 17, 173-265. 
Arega, F., Lee, J. H. W., and Tang, H. W. (2008). "Hydraulic jet control for river 
junction design of Yuen Long Bypass Floodway, Hong Kong." Journal of Hydraulic 
Engineering-Asce, 134(1), 23-33. 
Bagatur, T. (2007). "Modified Newton-Raphson solution for dispersion equation 
of transition water waves." Journal of Coastal Research, 23(6), 1588-1592. 
Bajo, M., Zampato, L., Umgiesser, G., Cucco, A., and Canestrelli, P. (2007). "A 
finite element operational model for storm surge predicition in Venice." Estuarine 
Coastal and Shelf Science, 75(1-2), 236-249. 
Baldwin, B.S., Lomax H. (1978). ?Thin-layer approximation and algebraic model 
for separated turbulent flows.? AISS paper, No. 78-257. 
Baldwin, B.S., Barth T.J. (1990). ?A one-equation turbulence transport model for 
high Reynolds number wall-bounded flows.? NASA TM-102847. 
Barros, S. R. M., and Garcia, C. I. (2007). "A global finite-difference semi-
Lagrangian model for the adiabatic primitive equations." Journal of Computational 
Physics, 226(2), 1645-1667. 
Bisson, P.A., et al. (1987). Large woody debris in forested streams in the Pacific 
Northwest, past, present and future; Contribution No.57. in Streamside management, 
forestry and fishery interactions. E. O. Salo and T. W. Cundy. Eds., Coll. Of  Forest 
Resour, Univ. of Washington, Seattle, Wash., 143-190. 
Blayo, E. and Debreu, L. (2005). ?Revisting open boundary conditions from the 
point of view of characteristic variables.? Ocean Modeling, 9, pp. 231-252. 
Blumberg, A. F., and Mellor, G. L. (1987). A description of a three-dimensional 
coastal ocean circulation model, In Three-Dimensional Coastal Ocean Models, 
American Geophysical Union, Washington, DC. 
Bonaventura, L., and Rosatti, G. (2002). "A cascadic conjugate gradient algorithm 
for mass conservative, semi-implicit discretization of the shallow water equations on 
 125
locally refined structured grids." International Journal for Numerical Methods in Fluids, 
40(1-2), 217-230. 
Bourchtein, A., and Bourchtein, L. (2006). "Modified Time Splitting Scheme for 
Shallow Water Equations." Mathematics and Computers in Simulation, 73, 52-64. 
Bourchtein, A., and Bourchtein, L. (2007). "Semi-Lagrangian semi-implicit time-
splitting scheme for the shallow water equations." International Journal for Numerical 
Methods in Fluids, 54(4), 453-471. 
Brown, P. N., and Saad, Y. (1990). "Hybrid Krylov methods for nonlinear 
systems of equations." SIAM J. Sci. Stat. Comput., 11, 450-481. 
Burguete, J., and Garcia-Navarro, P. (2004). "Implicit schemes with large time 
step for non-linear equations: application to river flow hydraulics." International Journal 
for Numerical Methods in Fluids, 46(6), 607-636. 
Burguete, J., Garcia-Navarro, P., and Murillo, J. (2006). "Numerical boundary 
conditions for globally mass conservative methods to solve the shallow-water equations 
and applied to river flow." International Journal for Numerical Methods in Fluids, 51(6), 
585-615. 
Cao, J. W., and Sun, J. C. (2005). "An efficient and effective nonlinear solver in a 
parallel software for large scale petroleum reservoir simulation." International Journal of 
Numerical Analysis and Modeling, 2, 15-27. 
Cao, Z., and Carling, P. A. (2002). "Mathematical modelling of alluvial rivers: 
reality and myth. Part 2: Special issues." Proceedings of the Institution of Civil 
Engineers-Water and Maritime Engineering, 154(4), 297-307. 
Casulli, V. (1990). "Semi-implicit Finite Difference Methods for the Two-
Dimensional Shallow Water Equations." Journal of Computational Physics, 86, 56-74. 
Casulli, V., and Cattani, E. (1994). "Stability, Accuracy and Efficiency of a 
Semiimplicit Method for 3-Dimensional Shallow-Water Flow." Computers & 
Mathematics with Applications, 27(4), 99-112. 
Casulli, V., and Cheng, R. T. (1992). "Semiimplicit Finite-Difference Methods for 
3-Dimensional Shallow-Water Flow." International Journal for Numerical Methods in 
Fluids, 15(6), 629-648. 
Celia, M. A., Ahuja, L. R., and Pinder, G. F. (1987). "Orthogonal collocation and 
alternating-direction procedures for unsaturated flow problems " Adv. Water Resour., 10, 
178-187. 
Chan, T. F., and Jackson, K. R. (1984). "Nonlinearly preconditioned Krylov 
subspace methods for discrete Newton algorithms." SIAM J. Sci. Stat. Comput., 5, 533-
542. 
Chen, S. C., and Peng, S. H. (2006). "Two-dimensional numerical model of two-
layer shallow water equations for confluence simulation." Advances in Water Resources, 
29(11), 1608-1617. 
Cho, Y.-S., and Yoon, S. B. (1998). "A Modified Leap-Frog Scheme for Linear 
Shallow-Water Equations." Costal Engineering Journal, JSCE, 40(2), 191-205. 
Clement, T. P., Wise, W. R., and Molz, F. J. (1994). "A Physically-Based, 2-
Dimensional, Finite-Difference Algorithm for Modeling Variably Saturated Flow." 
Journal of Hydrology, 161(1-4), 71-90. 
 126
Cullen, M. J. P. (2001). "Alternative implementations of the semi-Lagrangian 
semi-implicit schemes in the ECMWF model." Quarterly Journal of the Royal 
Meteorological Society, 127(578), 2787-2802. 
Daily, C., and Imberger, J. (2003). "Modelling solitons under the hydrostatic and 
Boussinesq approximations." International Journal for Numerical Methods in Fluids, 
43(3), 231-252. 
Davis, J. A. (1986). Boundary Layers, Flow Microenvironments and Stream 
Benthos. Limnology in Australia (eds De Deck and Williams), pp. 293-312. CSIRO 
Australia, Melbourne. Dr. W. Junk Publishers, Doedercht. 
Dean, R. G. and R. A. Dalrymple (1998). Water Wave Mechanics for Engineering 
and Scientists. World Scientific Publishing Co. Pte. Ltd., New Jersey. 
Dedong, L., Zhongbo, Y., Zhenchun, H., Chuanguo, Y., and Qin, J. (2007). 
"Groundwater simulation in the Yangtze River basin with a coupled climate-hydrologic 
model." Journal of China University of Geosciences, 18, 155-157. 
Delis, A. I., and Katsaounis, T. (2005). "Numerical solution of the two-
dimensional shallow water equations by the application of relaxation methods." Applied 
Mathematical Modelling, 29(8), 754-783. 
Dettmer, W. G., and Peric, D. (2007). "A fully implicit computational strategy for 
strongly coupled fluid-solid interaction." Archives of Computational Methods in 
Engineering, 14(3), 205-247. 
Dubois, T., Jauberteau, F., Temam, R. M., and Tribbia, J. (2005). "Multilevel 
schemes for the shallow water equations." Journal of Computational Physics, 207, 660-
694. 
Durran, D. R. (1999). Numerical Methods for Wave Equations in Geophysical 
Fluid Dynamics. Springer-Verlag, New York. 
Dutta, D., Alam, J., Umeda, K., Hayashi, M., and Hironaka, S. (2007). "A two-
dimensional hydrodynamic model for flood inundation simulation: a case study in the 
lower Mekong river basin." Hydrological Processes, 21(9), 1223-1237. 
Fennema, R. J., and Chaudhry, M. H. (1990). "Explicit Methods for 2-D Transient 
Free-Surface Flows." Journal of Hydraulic Engineering, 116(8), 1013-1035. 
Ferziger J.H. (1996). Simulation and Modeling of Turbulent Flows: Large eddy 
simulation. M.Y. Hussaini, T. Gatski (eds.), Cambridge Univ. Press, New York. 1996 
Ferziger J. H., Koseff, J.R., Monismith, S.G. (2002). ?Numerical simulation of 
geophysical turbulence.? Computers & Fluids, 31, 557-568. 
Ferziger, J. H., and Peric, M. (1999). Computational Methods for Fluid 
Dynamics, Springer, New York.Fleckenstein, J. H., Niswonger, R. G., and Fogg, G. E. 
(2006). "River-aquifer interactions, geologic heterogeneity, and low-flow management." 
Ground Water, 44(6), 837-852. 
Fu, S. and Hodges, B. R. (2005). ?Grid-scale dependency of subgrid-scale 
structure effects in hydraulic models of rivers and streams.? Proc., 2005 Mechanics and 
Materials Conference (McMat 2005), Louisiana State University, Baton Rouge. 
Fujima, K. and Shigemura, T (2000). ?Determination of Grid Size for Leap-Frog 
Finite Difference Model to Simulate Tsunamis around A Conical Island.? Costal 
Engineering, 42(2), 197-210. 
 127
Garcia, R., and Kahawita, R. A. (1986). "Numerical Solution of the St. Venant 
Equations with the MacCormack Finite-Difference Scheme." International Journal for 
Numerical Methods in Fluids, 6, 259-274. 
Gejadze, I. Y., and Monnier, J. (2007). "On a 2D 'zoom' for the 1D shallow water 
model: Coupling and data assimilation." Computer Methods in Applied Mechanics and 
Engineering, 196(45-48), 4628-4643. 
Gippel, C. J.(1995). ?Environmental hydraulics of large woody debris in streams 
and rivers.? J. Environmental Engineering, 121, 388-395. 
Giri, S., and Shimizu, Y. (2006). "Numerical computation of sand dune migration 
with free surface flow." Water Resources Research, 42(10). 
Gr?tzner, A. et al. (1996). ?The Impact of Sub-grid Scale Sea-ice 
Inhomogeneities on the Performance of the Atmospheric General Circulation Model 
ECHAM3.? Climate Dynamics,12,477-496. 
Grossman, G.D., Rincon, P. A., Farr, M.D., and Ratajczak, R. E. (2002). ?A new 
optimal foraging model predicts habitat use by drift-feeding stream minnows.? Ecol. 
Freshwat. Fish, 11,2-10. 
Hamrick, J. M. (1992). "A three-dimensional environmental fluid dynamics 
computer code: theoretical and computarional aspects." Virginia Institute of Marine 
Science, School of Marine Science, The College of William and Mary, Special Report 
317., Gloucester Point, VA 23062. 
Harmon, M. E., et al. (1986). ?Ecology of coarse woody debris in temperature 
ecosystems.? Adv. in Ecological Res, 15, 133-302. 
Hirose, N., Asai, K., Ikawa, K., and Kawamura, R. (1991). "Euler Flow-Analysis 
of Turbine Powered Simulator and Fanjet Engine." Journal of Propulsion and Power, 
7(6), 1015-1022. 
Hodges, B. R. (1997). "Numerical simulation of nonlinear free-surface waves on a 
turbulent open-channel flow," Ph.D. dissertation, Department of Civil Engineering, 
Stanford University. 
Hodges, B. R. (2000). "Numerical Techniques in CWR-ELCOM (code release 
v.1)." Centre Water Res., Nedlands, Western Australia, Australia. 
Hodges, B. R., B. Laval, Wadzuk, B. (2006). ?Numerical error assessment and 
temporal horizon for internal waves in a hydrostatic model.? Ocean Modeling, 13(1), 44-
64. 
Hunter, N. M., Bates, P. D., Neelz, S., Pender, G., Villanueva, I., Wright, N. G., 
Liang, D., Falconer, R. A., Lin, B., Waller, S., Crossley, A. J., and Mason, D. C. (2008). 
"Benchmarking 2D hydraulic models for urban flooding." Proceedings of the Institution 
of Civil Engineers-Water Management, 161(1), 13-30. 
Iskandarani, M., Levin, J. C., Choi, B. J., and Haidvogel, D. B. (2005). 
"Comparison of advection schemes for high-order h-p finite element and finite volume 
methods." Ocean Modelling, 10(1-2), 233-252. 
Kadalbajoo, M. K., and Awasthi, A. (2006). "A numerical method based on 
Crank-Nicolson Scheme for Burgers' Equation." Applied Mathematics and Computation, 
182, 1430-1442. 
 128
Kar, S. K. (2006). "A Semi-Implicit Runge-Kutta Time-Difference Scheme for 
the Two-Dimensional Shallow Water Equations." Monthly Weather Review, 
134(October), 2916-2926. 
Kazezy?lmaz-Alhana, C. M., Medina, M. A., and Raob, P. (2005). "On numerical 
modeling of overland flow." Applied Mathematics and Computation  166(3), 724-740. 
Kutluay, S., Bahadir, A. R., and Ozdes, A. (1999). "Numerical solution of one-
dimensional Burgers equation: explicit and exact-explicit finite difference methods." 
Journal of Computational and Applied Mathematics, 103, 251-261. 
Keller, E.A., and Swanson, F.J. (1979). ?Effect of large organic material on 
channel form and fluvial processes.? Earth Surface Processes, 4, 361-380. 
Keshari, A. K., and Koo, M. H. (2007). "A numerical model for estimating 
groundwater flux from subsurface temperature profiles." Hydrological Processes, 21(25), 
3440-3448. 
Knoll, D. A., and Keyes, D. E. (2004). "Jacobian-free Newton-Krylov methods: a 
survey of approaches and applications." Journal of Computational Physics, 193(2), 357-
397. 
Kolmogorov, AN. (1942). ?The equations of turbulent motion in an 
incompressible fluid.? Isv. Acad. Sci., USSR, Phys. 6, 56-58. 
Kurbatskii, K. A., and Mankbadi, R. R. (2004). "Review of computational 
aeroacoustics algorithms." International Journal of Computational Fluid Dynamics, 
18(6), 533-546. 
Kwag, S. H. (2000). "Computation of water and air flow with submerged 
hydrofoil by interface capturing method." Ksme International Journal, 14(7), 789-795. 
Kwok, F., and Tchelepi, H. (2007). "Potential-based reduced Newton algorithm 
for nonlinear multiphase flow in porous media." Journal of Computational Physics, 
227(1), 706-727. 
Lomax, H., Pulliam, T. H., and Zingg, D. W. (1999). Fundamentals of 
Computational Fluid Dynamics, Springer, New York. 
Le, V. S., Yamashita, T., Okunishi, T., Shinohara, R., and Miyatake, M. (2006). 
"Characteristics of suspended sediment material transport in the Ishikari Bay in snowmelt 
season." Applied Ocean Research, 28(4), 275-289. 
Lehmann, F., and Ackerer, P. (1998). "Comparison of Iterative Methods for 
Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media." 
Transport in Porous Media, 31, 275-292. 
Li, C. W. (1993). "A simplified Newton iteration method with linear finite 
elements for transient unsaturated flow." Water Resour. Res., 29(4), 965-971. 
Li, G. Y., and Jackson, C. R. (2007). "Simple, accurate, and efficient revisions to 
MacCormack and Saulyev schemes: High Peclet numbers." Applied Mathematics and 
Computation, 186(1), 610-622. 
Liao, C. B., Wu, M. S., and Liang, S. J. (2007). "Numerical simulation of a dam 
break for an actual river terrain environment." Hydrological Processes, 21(4), 447-460. 
Lomax, H., Pulliam, T. H. and Zingg, D.W. (1996) Fundamentals of 
Computational Fluid Dynamics. Springer, New York Inc. 
Luo, Q. (2007). "A distributed surface flow model for watersheds with large water 
bodies and channel loops." Journal of Hydrology, 337(1-2), 172-186. 
 129
MacCormack, R. W. (1969). "The effect of viscosity in hypervelocity impact 
cratering." AIAA-1969-354. 
Marzolf, G. R. (1978). The potential effects of clearing and snagging on stream 
ecosystems. FWS/OBS-78-14, U.S. Dept of Interior. Fish and Wildlife Service, Nat. 
Stream Alteration Team, Washington, D.C. 
McMahon, T. E., and Hartman, G. F. (1989). ?Influence of cover complexity and 
current velocity on winter habitat use by juvenile Coho Salmon (Oncorhynchus kisutch).? 
Can. J. of Fisheries and Aquatic Sci., 46, 1551-1557. 
Mendez-Nunez, L. R., and Carroll, J. J. (1993). "Comparison of Leapfrog, 
Smolarkiewicz, and MacCormack Schemes Applied to Nonlinear Equations." Mon. 
Weather Rev., 121(February), 565-578. 
Miropol'sky, Y. Z. (2001). Dynamics of internal gravity waves in the ocean, 
Kluwer. 
Mohebalhojeh, A. R., and Dritschel, D. G. (2007). "Assessing the numerical 
accuracy of complex spherical shallow-water flows." Monthly Weather Review, 135(11), 
3876-3894. 
Moin, P. (2001). Fundamentals of Engineering Numerical Analysis. Cambridge, 
U.K., Cambridge University Press. 
Mousseau, V. A., Knoll, D. A., and Reisner, J. M. (2002). "An implicit 
nonlinearly consistent method for the two-dimensional shallow-water equations with 
coriolis force." Monthly Weather Review, 130(11), 2611-2625. 
Murillo, J., Garcia-Navarro, P., and Burguete, J. (2008). "Analysis of a second-
order upwind method for the simulation of solute transport in 2D shallow water flow." 
International Journal for Numerical Methods in Fluids, 56(6), 661-686. 
Nguyen, D. K., Shi, Y. E., Wang, S. S. Y., and Nguyen, H. (2006). "2D shallow-
water model using unstructured finite-volumes methods." Journal of Hydraulic 
Engineering-Asce, 132(3), 258-269. 
Niet, A. d., Wubs, F., Scheltinga, A. T. v., and Dijkstra, H. A. (2007). "A tailored 
solver for bifurcation analysis of ocean-climate models." Journal of Computational 
Physics, 227 (1), 654-679. 
Paniconi, C., Aldama, A. A., and Wood, E. F. (1991). "Numerical Evaluation of 
Iterative and Noniterative Methods for the Solution of the Nonlinear Richards Equation." 
Water Resources Research, 27(6), 1147-1163. 
Paniconi, C., and Putti, M. (1994). "A Comparison of Picard and Newton Iteration 
in the Numerical Solution of Multidimensional Varibly Saturated Flow Problems." Water 
Resources Research, 30(12), 3357-3374. 
Peyret, R., and Taylor, T. D. (1983). Computational Methods for Fluid Flow, 
Springer, New York. 
Pozrikidis, C. (2005). Introduction to finite and spectral element methods using 
MATLAB, Chapman & Hall/CRC, Boca Raton. 
Prandtl L. (1925). Uber die ausgebildete turbulenz. ZAMM 5, pp. 331-340. 
Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T. (1992). 
Numerical recipes in FORTRAN : the art of scientific computing, Cambridge University 
Press, New York.Peyret R. and Taylor, T. D. (1983). Computational Method s for Fluid 
Flow. Springer-Verlag, New York Inc.  
 130
Rabeni. C. F., and Jacobson, R. B. (1993). ?The importance of fluvial hydraulics 
to fish-habitat restoration in low-gradient alluvial streams.? Freshwater Biology, 29, 211-
220.  
Reisner, J. M., Mousseau, V. A., Wyszogrodzki, A. A., and Knoll, D. A. (2005). 
"An implicitly balanced hurricane model with physics-based preconditioning." Monthly 
Weather Review, 133(4), 1003-1022. 
Robert, L. W. (1992). Oceanographical Enginereing, Dover Publications, Inc., 
Mineola, New York. 
Rosatti, G., Cesari, D., and Bonaventura, L. (2005). "Semi-implicit, semi-
Lagrangian modelling for environmental problems on staggered Cartesian grids with cut 
cells." Journal of Computational Physics, 204(1), 353-377. 
Rueda, F. J., and Schladow, G. S. (2002). "Quantitative Comparison of Models 
for Barotropic Response of Homogeneous Basins." Journal of Hydraulic Engineering, 
128(2), 201-213. 
Salvetti, MV, Banerjee, S. (1995). ?A priori tests of a new dynamic subgrid-scale 
model for finite- difference large-eddy simulations.? Phys. Fluids, 7(11), 2831-2847. 
Schoenstadt, A. (1980). "A transfer analysis of numerical schemes used to 
simulate geostrophic adjustment." Mon. Weather Rev., 108, 1248-1295. 
Shaw, D., Martz, L. W., and Pietroniro, A. (2005). "Flow routing in large-scale 
models using vector addition." Journal of Hydrology, 307(1-4), 38-47. 
Shields, F.D. Jr., and Nunnally, N.R. (1984). ?Environmental aspects of clearing 
and snagging.? Journal of Environmental Engineering, 110 (1), 152-165. 
Spalart, P.R, Allmaras, S.R. (1992). ?A one-equation turbulence model for 
aerodynamic flow.? AIAA J. 29, 1819-1835. 
Spitaleri, R. M., and Corinaldesi, L. (1997). "A multigrid semi-implicit finite 
difference method for the two-dimensional shallow water equations." International 
Journal for Numerical Methods in Fluids, 25(11), 1229 - 1240. 
Stelling, G. S., and Duinmeijer, S. P. A. (2003). "A staggered conservative 
scheme for every Froude number in rapidly varied shallow water flows." International 
Journal for Numerical Methods in Fluids, 43(12), 1329-1354. 
Steppeler, J. (2006). "HLI, a direct method suitable for partial and fully implicit 
time integration of primitive equation meteorological models." Computers & 
Mathematics with Applications, 52(8-9), 1357-1372. 
Sullivan, K., Lisle, T. E., Dolloff, C. A., Grant, G. E., and Reid, L. M. (1987). 
?Stream channels: the link between forest and fishes.? Streamside management, forestry 
and fishery interactions. E. O. Salo and T. W. Cundy. Eds., Coll. Of Forest Resour., 
Univ. of Washington, Seattle, Wash., 39-97. 
Szymkiewicz, R. (1992). "A Mathematical-Model of Storm-Surge in the Vistula 
Lagoon, Poland." Coastal Engineering, 16(2), 181-203. 
Tannehill, J. C., Anderson, D. A., and Pletcher, R. H. (1997). Computational fluid 
mechanics and heat transfer, Taylor & Francis, Washington, DC. 
Thomee, V. (2001). "From finite differences to finite elements - A short history of 
numerical analysis of partial differential equations." Journal of Computational and 
Applied Mathematics, 128(1-2), 1-54. 
 131
Thouvenin, B., Gonzalez, J. L., Chiffoleau, J. F., Boutier, B., and Le Hir, P. 
(2007). "Modelling Pb and Cd dynamics in the Seine estuary." Hydrobiologia, 588, 109-
124. 
Tsynkov, S. V. (1998).  ?Numerical solution of problems on unbounded domains. 
A review.? Applied Numerical Mathematics, 27, 465-532. 
Turek, S. (1996). "A comparative study of time-stepping techniques for the 
incompressible Navier-Stokes equations: From fully implicit non-linear schemes to semi-
implicit projection methods." International Journal for Numerical Methods in Fluids, 
22(10), 987-1011. 
Vinayan, V. (2003). Boundary-Integral Analysis of Nonlinear Diffraction Forces 
on a Submerged Body. Master thesis, Department of Ocean Engineering, Florida Atlantic 
University. 
Vreugdenhil, C. B. (1994). Numerical Methods For Shallow-Water Flow, Kluwer 
Academic, Dordrecht, The Netherlands. 
Wadzuk, B. M. (2004). "Hydrostatic and non-hydrostatic internerwave models," 
Ph.D dissertation, The Univeristy of Texas at Austin. 
Weerakoon, S. B., Tamai, N., and Kawahara, Y. (2003). "Depth-averaged flow 
computation at a river confluence." Proceedings of the Institution of Civil Engineers-
Water and Maritime Engineering, 156(1), 73-83. 
Webster, R. (2007). "Algebraic multigrid and incompressible fluid flow." 
International Journal for Numerical Methods in Fluids, 53(4), 669-690. 
Wiegel, R. L. (1964). Oceanographical Engineering. Dover Publication, Inc., 
Mineola, New York. 
Wilders, P., Van Stijn, T. L., Stelling, G. S., and Fokkema, G. A. (1988). "A Fully 
Implicit Splitting Method for Accuracy Tidal Computations." International Journal for 
Numerical Methods in Engineering, 26, 2707-2721. 
Wubs, F. W., de Niet, A. C., and Dijkstra, H. A. (2006). "The performance of 
implicit ocean models on B- and C-grids." Journal of Computational Physics, 211(1), 
210-228. 
Yuan, H., and Wu, C. (2004). "An implicit three-dimensional fully non-
hydrostatic model for free-surface flows." International Journal for Numerical Methods 
in Fluids, 44(8), 811-835  
Zhang, Y., Baptista, A. M., and III, E. P. M. (2004). "A cross-scale model for 3D 
baroclinic circulation in estuary?plume?shelf systems: I. Formulation and skill 
assessment." Continental Shelf Research, 24(18), 2187-2214. 
Zhou, J. T., Lin, J. G., and Xie, Z. H. (2007). "A compact explicit difference 
scheme of high accuracy for extended Boussinesq equations." China Ocean Engineering, 
21(3), 507-514. 
Zhou, W. (2002). "An Alternative Leapfrog Scheme for Surface Gravity Wave 
Equations." Journal of Atmospheric and Oceanic Technology, 19(September), 1415-
1423.