# Modified Predictor-Corrector WAF Method for the Shallow Water Equations with Source Terms

Authors: Montri Maleewong
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

A modified
predictor-corrector scheme combining with the
depth gradient method
(DGM) and the weighted average flux (WAF) method has been
presented to solve the one-dimensional shallow water equations
with source terms. Approximate solutions in the predictor step are
obtained by the DGM with piecewise-linear reconstructions in each
cell volume. The source terms can then be calculated directly by
these predicted values at the corresponding half-time step. In the
corrector step, the TVD version of the WAF method is applied to
calculate the numerical fluxes at the same half-time step for each
cell face. The accuracy of numerical solutions is shown by
applying the method to solve various test cases in both steady and
unsteady problems with and without source terms. It
shows that the numerical results are in good agreement with the
existing analytical solutions as well as experimental data in some
test cases.


---

## Body

## 1. Introduction

The shallow water equations have a wide variety of applications in ocean and hydraulic engineering. Some examples are tides in oceans and moving waves in shallow beaches, as well as flood waves in rivers. Due to the nonlinear behavior of the equations, analytical method can only be successful in very special cases of flow. Numerical methods must be used to solve approximately in case of realistic problems. A number of finite volume methods based on Godunov type have been developed to solve the shallow water equations in various forms. Difficulties in the computation arise when some source terms are introduced in the model equations. The source terms are in the form of bed and friction slopes.  A simple and efficient method for solving the equations with source terms is a split-step method [1, 2] in which the nonhomogeneous equations are splitt into a homogeneous equation and a set of ordinary differential equations dealing with only the source and time derivative terms; see [3]. By this method, the scheme is not easy to be improved to achieve high-order accuracy in time discretization because the source terms cannot be computed directly during the semitime integration step. Recently, the operator split-step method related with the Cauchy-Kowalewski method is applied to solve the thermal-peaking wave in open-channel flows, see [4], where the fluxes at interfaces are calculated by the weighted average method.In this study, the predictor-corrector method for time integration is proposed. We present a modified concept by applying both the depth gradient method (DGM) in [5] and the weighted average flux method (WAF) method in [1] to solve numerically the one-dimensional shallow water equations. The approximate solutions in the predictor step which corresponds to the half-time step are obtained by the DGM; see Figure 1. The approximations are performed by the piecewise-linear reconstruction to calculate the conservative variables at cell interfaces in each cell volume. Then, these predicted solutions are employed to calculate the values of source terms at the corresponding half-time step. This is a different point from the usual operator splitting of the WAF method. The problems are then become how we can calculate the numerical fluxes at the corresponding half-time step for each cell face. Various approaches can be applied to resolve this problem. One famous method is the Roe averaging method, see [6], but it is usually suffer from some artificial oscillations near the discontinuity region. The slope limiter concept is needed to remove the oscillations. In this work, the WAF method with the total variation diminishing (TVD) approach [7] is applied to obtain high-order solutions and to remove some oscillations in case of very high gradient problem. This is a modified approach. Consequently, the shallow water equations with source terms can be solved numerically in this work.Figure 1
Computational diagram in thext-plane.The shallow water equations are shown in Section2. The predictor-corrector method with DGM and WAF are provided in Section 3. Accuracy of numerical results for both steady and unsteady flows with and without source terms is demonstrated in Section 4. Conclusions are made in Section 5.


## 2. The Shallow Water Equations

For a rectangular open-channel flow, the one-dimensional unsteady flow under the hydrostatic pressure and small channel slope assumptions can be described in conservative form as(2.1)∂U∂t+∂F∂x=S,
where x and t represent longitudinal distance and time, respectively, and (2.2)U=[huh],F=[huhu2+12gh2],
where h is the fluid depth, u is the horizontal velocity, and g is the gravitation acceleration.The source termS in (2.1) represents bed slope and friction slope at the bottom of the channel(2.3)S=[0gh(S0-Sf)],
where S0 is the slope of the bottom elevation zb, and Sf is the slope of the energy grade line, (2.4)Sf=nM2u2h4/3,
where nM is the Manning’s roughness coefficient.


## 3. Predictor-Corrector WAF with DGM

In this section, we present an accurate and efficient method that combines two important concepts which are the weighted average flux method in [1], and the depth gradient method modified from [5], with the predictor-corrector TVD scheme to solve numerically the shallow water equation (2.1) in various test problems.The computational domain is discretized asxi=iΔx, and tn=nΔt, i=1,…N, N is the number of cells, and n is the number of time steps. In Figure 1, the cell volume and time step are given by (3.1)Δx=xi+1/2-xi-1/2,Δt=tn+1-tn.
Left and right interfaces for each cell are xi-(1/2) and xi+(1/2), respectively.The integral averages ofU(x,t) at time t=tn and t=tn+1 over the cell volume are given by(3.2)Uin=1Δx∫xi-1/2xi+1/2U(x,tn)dx,Uin+1=1Δx∫xi-1/2xi+1/2U(x,tn+1)dx.
Time integral averaged flux at interfaces are given by(3.3)Fi-1/2=1Δt∫tntn+1F(U(xi-1/2,t))dt,Fi+1/2=1Δt∫tntn+1F(U(xi+1/2,t))dt.
The splitting scheme cooperating with the source term can be expressed by(3.4)Uin+1=Uin-ΔtΔx[Fi+1/2n-Fi-1/2n]+ΔtS(Ui(s)).
The source term is approximated at Ui(s). There are two obvious choices which are Ui(s)=Uin and Ui(s)=Uin+1, or a linear combination of these two values, (3.5)Ui(s)=(1-β)Uin+βUin+1,
where β is a weighted constants in which 0≤β≤1.The scheme (3.4) is just the first-order method in time. We can improve the scheme by applying the predictor-corrector method. There are two time levels with Δt/2. The computations can be separated into two time steps as(3.6)Ui*=Uin-Δt2Δx[Fi+1/2DGM-Fi-1/2DGM]+ΔtS(Uin),(3.7)Uin+1=Ui*-Δt2Δx[Fi+1/2TVD-Fi-1/2TVD]+ΔtS(Ui*),
where Fi±(1/2)DGM and Fi±(1/2)TVD are the numerical fluxes at interfaces obtained by the depth gradient method and the TVD version of WAF method, respectively.Solving (3.6) provides the values of the conservative variables Ui* at t=tn+Δt/2. The values of fluxes at cell interfaces on the RHS of (3.6) are approximated by the method of DGM described in the next section. Then, we use these predicted solutions Ui* to calculate the source terms on the RHS of (3.7). The values of fluxes at cell interfaces on the RHS of (3.7) are approximated by the WAF method which is in the process of the half-time step; see Figure 1. Finally, the corrected solutions at t=tn+Δt can be obtained. So, the time integration process is performed repeatly until the final time reached. All of computational details of the present method are summarized in the next sections.


### 3.1. Predictor Step with DGM

We apply the DGM for solvingUi* in (3.6). This scheme is an accurate reconstruction of conservative variables at cell interfaces. Then, numerical fluxes at cell faces can be calculated directly using the values of predicted-conservative variables. The values of conservative variables within a cell are calculated using a piecewise-linear reconstruction, that is, for any function ϕ within a cell i, the approximation is given by(3.8)ϕ=ϕi+(x-xi)δ(ϕi),
where δ is the gradient of ϕ given by (3.9)δϕi=G(ϕi+1-ϕixi+1-xi,ϕi-ϕi-1xi-xi-1),
where G is a slope limiter used to avoid spurious oscillations at the cell interfaces [5]. In our work, we apply G as the minmod limiter which is in the form of(3.10)G(a,b)=max⁡(0,min⁡(a,b)).
Thus, the values of ϕ on the left and right of the considering cell interface (i-1/2) are (3.11)ϕi-1/2L=ϕi-1+12Δxi-1δϕi-1,ϕi-1/2R=ϕi-12Δxiδϕi.
The values of ϕi+1/2L and ϕi+1/2R can be calculated in the same way.We apply the conservative variablesh and u as the function ϕ in (3.8). The conservative variables can then be approximated at the cell interfaces for each cell. Using (2.2), we can calculate U and F at the cell faces. After substituting these values into (3.6), we can obtain Ui* such that the source term is calculated directly by the values of h and u in each cell. This completes the calculations of the predictor step by the depth gradient method.


### 3.2. Corrector Step with WAF

In the corrector step, we need to know the approximate values of fluxes at interfaces at time steptn+Δt/2. Then we can march the numerical solutions in time to the next time step tn+1. In this work, we apply the WAF method to obtain the numerical fluxes at time step tn+Δt/2. This is a second-order accurate method in space; see [1]. Main concepts of this method are summarized in this subsection.In each cell volume, the Riemann problem must be solved locally in order to obtain numerical fluxes at interfaces. We approximate these fluxes by the HLL (Harten, Lax, and van Leer) approach; see [8]. Note that another related approach is HLLC approach which needs to divide the region at interface into four regions. The HLL numerical flux can be written as(3.12)Fi+1/2HLL={FLifSL≥0,SRFL-SLFR+SRSL(UR-UL)SR-SLifSL≤0≤SR,FRifSR≤0.
Here, SL and SR are the signal velocities in the solutions of the Riemann problem with known data UL≡Uin and UR≡Ui+1n, and the corresponding fluxes FL≡Fn(UL) and FR≡Fn(UR). The signal velocities must be estimated separately for two cases of wet and dry beds. More detail derivations can be found in [9]. Summarized details of the approximations are given in the next subsections.


#### 3.2.1. Wet-Bed Approximation

Assuming that both the left and right waves are rarefaction waves, the wave speed estimatesSL and SR are(3.13)SL=uL-aLqL,SR=uR+aRqR,
where qK with K is L or R given by (3.14)qK={12[(h*+hK)h*hK2]ifh*>hK,1ifh*≤hK,
where h* is an approximation of the exact solution for h in the interface area or usually called as star region. The approximations are given by(3.15)h*=1g[12(aL+aR)+14(uL-uR)]2,u*=12(uL+uR)+aL-aR,
where aK=ghK for K is L or R.


#### 3.2.2. Dry-Bed Approximation

The computations of wave speed to the exact dry front are given by(3.16)SL={uR-2aRifhL=0,usualestimateifhL>0,SR={uL+2aLifhR=0,usualestimateifhR>0.
Usual estimate means that the computations of wave speed are obtained from (3.13).


#### 3.2.3. TVD Version of WAF Method

The basic idea of the weighted average flux method is described in this subsection. This method provides numerical fluxes at interfaces for the next half-time step. The method is performed by applying cell integration along the interface of the half-time step solution,(3.17)Fi+(1/2)WAF=1Δx∫-Δx/2Δx/2F(Ui+(1/2)(x,Δt2))dx.
In the form of wave structure, the integral (3.17) can be written as the summation(3.18)Fi+(1/2)WAF=∑k=1NC+1wkFi+(1/2)(k),
where the weights wk are(3.19)wk=12(ck-ck-1)withc0=-1,cNC+1=1,(3.20)ck=ΔtSkΔx
is the Courant number of wave k.For this method,Sk is the speed of wave k and NC is the number of conservation laws, or the number of waves in the solution of the Riemann problem.Substituting the weights (3.19) into (3.18), an alternative expression can be obtained,(3.21)Fi+1/2WAF=12(Fi+Fi+1)-12∑k=1NCckΔFi+1/2(k),
with (3.22)ΔFi+1/2(k)=Fi+1/2(k+1)-Fi+1/2(k)
representing the flux jump across the wave k.For a one-dimensional flow, we haveNC=2. So, we need to know Fi+(1/2)(1), Fi+(1/2)(2), and Fi+(1/2)(3) in order to compute weighted average flux in (3.21). Left and right areas can be represented by k=1 and k=3, respectively. We have already known the values of conservative variables in these areas from the initial conditions, except that in the interface area which corresponds to k=2. In this work, the computation by (3.12) is applied to calculate the numerical flux at the interface region.Following [1], a high-resolution method can be obtained by enforcing a total variation diminishing (TVD) constraint to the weighted numerical flux at the interface. The resulting TVD version of the WAF method can be expressed by(3.23)Fi+1/2TVD=12(Fi+Fi+1)-12∑k=1Nsgn⁡(ck)ϕi+1/2(k)ΔFi+1/2(k),
where ϕi+1/2(k) is a WAF limiter function. There are several choices of this limiter function. In this work, we apply the minmod limiter function(3.24)ϕi+1/2(k)=ϕi+1/2(r(k),|ck|)={1,r(k)≤0,1-(1-|ck|)r(k),0≤r(k)≤1,|ck|,r(k)≥1.
The flow parameter r(k) is (3.25)r(k)={Δqi-1/2(k)Δqi+1/2(k)ifck≥0,Δqi+3/2(k)Δqi+1/2(k)ifck<0,
where q is a suitable variable. For shallow water flows, q is set to be h. Note that the value of Δqi-1/2(k) denotes the jump in q across wave k in the Riemann problem with data (Ui-1n,Uin). Hence, the approximation of flux at interface is finally obtained by (3.23). This completes the weighted average flux method for computing flux at interfaces for each cell volume.


#### 3.2.4. Calculation of Time Step

For any celli with length of Δxi, time step is defined by(3.26)Δti=CnΔxi(Si+1/2L-Si-1/2R),
where Si+1/2L and Si-1/2R are the left and right wave speed in the Riemann problem for cell i, and Cn is the Courant number coefficient.For easy programming in the case of uniform cell length, the time step can be defined by(3.27)Δt=CnΔxSmax⁡,
where Smax=max{Si+1/2L-Si-1/2R}, for all i.


#### 3.2.5. Time Integration in the Corrector Step

Numerical solutions for the next time step can be obtained by using (3.7) as follows. The numerical fluxes at interfaces are calculated by the flux slope limiter (3.23). The source term S involving the bed slope and friction effects is actually the function of the conservative variables h and u. We have already known the values of h and u at the half-time step from the predictor calculations. Hence, we can calculate all terms on the RHS of (3.7). Consequently, numerical solutions for the next time step are obtained. This completes the overall steps of the predictor-corrector method with the TVD version of the WAF method.


## 3.1. Predictor Step with DGM

We apply the DGM for solvingUi* in (3.6). This scheme is an accurate reconstruction of conservative variables at cell interfaces. Then, numerical fluxes at cell faces can be calculated directly using the values of predicted-conservative variables. The values of conservative variables within a cell are calculated using a piecewise-linear reconstruction, that is, for any function ϕ within a cell i, the approximation is given by(3.8)ϕ=ϕi+(x-xi)δ(ϕi),
where δ is the gradient of ϕ given by (3.9)δϕi=G(ϕi+1-ϕixi+1-xi,ϕi-ϕi-1xi-xi-1),
where G is a slope limiter used to avoid spurious oscillations at the cell interfaces [5]. In our work, we apply G as the minmod limiter which is in the form of(3.10)G(a,b)=max⁡(0,min⁡(a,b)).
Thus, the values of ϕ on the left and right of the considering cell interface (i-1/2) are (3.11)ϕi-1/2L=ϕi-1+12Δxi-1δϕi-1,ϕi-1/2R=ϕi-12Δxiδϕi.
The values of ϕi+1/2L and ϕi+1/2R can be calculated in the same way.We apply the conservative variablesh and u as the function ϕ in (3.8). The conservative variables can then be approximated at the cell interfaces for each cell. Using (2.2), we can calculate U and F at the cell faces. After substituting these values into (3.6), we can obtain Ui* such that the source term is calculated directly by the values of h and u in each cell. This completes the calculations of the predictor step by the depth gradient method.


## 3.2. Corrector Step with WAF

In the corrector step, we need to know the approximate values of fluxes at interfaces at time steptn+Δt/2. Then we can march the numerical solutions in time to the next time step tn+1. In this work, we apply the WAF method to obtain the numerical fluxes at time step tn+Δt/2. This is a second-order accurate method in space; see [1]. Main concepts of this method are summarized in this subsection.In each cell volume, the Riemann problem must be solved locally in order to obtain numerical fluxes at interfaces. We approximate these fluxes by the HLL (Harten, Lax, and van Leer) approach; see [8]. Note that another related approach is HLLC approach which needs to divide the region at interface into four regions. The HLL numerical flux can be written as(3.12)Fi+1/2HLL={FLifSL≥0,SRFL-SLFR+SRSL(UR-UL)SR-SLifSL≤0≤SR,FRifSR≤0.
Here, SL and SR are the signal velocities in the solutions of the Riemann problem with known data UL≡Uin and UR≡Ui+1n, and the corresponding fluxes FL≡Fn(UL) and FR≡Fn(UR). The signal velocities must be estimated separately for two cases of wet and dry beds. More detail derivations can be found in [9]. Summarized details of the approximations are given in the next subsections.


### 3.2.1. Wet-Bed Approximation

Assuming that both the left and right waves are rarefaction waves, the wave speed estimatesSL and SR are(3.13)SL=uL-aLqL,SR=uR+aRqR,
where qK with K is L or R given by (3.14)qK={12[(h*+hK)h*hK2]ifh*>hK,1ifh*≤hK,
where h* is an approximation of the exact solution for h in the interface area or usually called as star region. The approximations are given by(3.15)h*=1g[12(aL+aR)+14(uL-uR)]2,u*=12(uL+uR)+aL-aR,
where aK=ghK for K is L or R.


### 3.2.2. Dry-Bed Approximation

The computations of wave speed to the exact dry front are given by(3.16)SL={uR-2aRifhL=0,usualestimateifhL>0,SR={uL+2aLifhR=0,usualestimateifhR>0.
Usual estimate means that the computations of wave speed are obtained from (3.13).


### 3.2.3. TVD Version of WAF Method

The basic idea of the weighted average flux method is described in this subsection. This method provides numerical fluxes at interfaces for the next half-time step. The method is performed by applying cell integration along the interface of the half-time step solution,(3.17)Fi+(1/2)WAF=1Δx∫-Δx/2Δx/2F(Ui+(1/2)(x,Δt2))dx.
In the form of wave structure, the integral (3.17) can be written as the summation(3.18)Fi+(1/2)WAF=∑k=1NC+1wkFi+(1/2)(k),
where the weights wk are(3.19)wk=12(ck-ck-1)withc0=-1,cNC+1=1,(3.20)ck=ΔtSkΔx
is the Courant number of wave k.For this method,Sk is the speed of wave k and NC is the number of conservation laws, or the number of waves in the solution of the Riemann problem.Substituting the weights (3.19) into (3.18), an alternative expression can be obtained,(3.21)Fi+1/2WAF=12(Fi+Fi+1)-12∑k=1NCckΔFi+1/2(k),
with (3.22)ΔFi+1/2(k)=Fi+1/2(k+1)-Fi+1/2(k)
representing the flux jump across the wave k.For a one-dimensional flow, we haveNC=2. So, we need to know Fi+(1/2)(1), Fi+(1/2)(2), and Fi+(1/2)(3) in order to compute weighted average flux in (3.21). Left and right areas can be represented by k=1 and k=3, respectively. We have already known the values of conservative variables in these areas from the initial conditions, except that in the interface area which corresponds to k=2. In this work, the computation by (3.12) is applied to calculate the numerical flux at the interface region.Following [1], a high-resolution method can be obtained by enforcing a total variation diminishing (TVD) constraint to the weighted numerical flux at the interface. The resulting TVD version of the WAF method can be expressed by(3.23)Fi+1/2TVD=12(Fi+Fi+1)-12∑k=1Nsgn⁡(ck)ϕi+1/2(k)ΔFi+1/2(k),
where ϕi+1/2(k) is a WAF limiter function. There are several choices of this limiter function. In this work, we apply the minmod limiter function(3.24)ϕi+1/2(k)=ϕi+1/2(r(k),|ck|)={1,r(k)≤0,1-(1-|ck|)r(k),0≤r(k)≤1,|ck|,r(k)≥1.
The flow parameter r(k) is (3.25)r(k)={Δqi-1/2(k)Δqi+1/2(k)ifck≥0,Δqi+3/2(k)Δqi+1/2(k)ifck<0,
where q is a suitable variable. For shallow water flows, q is set to be h. Note that the value of Δqi-1/2(k) denotes the jump in q across wave k in the Riemann problem with data (Ui-1n,Uin). Hence, the approximation of flux at interface is finally obtained by (3.23). This completes the weighted average flux method for computing flux at interfaces for each cell volume.


### 3.2.4. Calculation of Time Step

For any celli with length of Δxi, time step is defined by(3.26)Δti=CnΔxi(Si+1/2L-Si-1/2R),
where Si+1/2L and Si-1/2R are the left and right wave speed in the Riemann problem for cell i, and Cn is the Courant number coefficient.For easy programming in the case of uniform cell length, the time step can be defined by(3.27)Δt=CnΔxSmax⁡,
where Smax=max{Si+1/2L-Si-1/2R}, for all i.


### 3.2.5. Time Integration in the Corrector Step

Numerical solutions for the next time step can be obtained by using (3.7) as follows. The numerical fluxes at interfaces are calculated by the flux slope limiter (3.23). The source term S involving the bed slope and friction effects is actually the function of the conservative variables h and u. We have already known the values of h and u at the half-time step from the predictor calculations. Hence, we can calculate all terms on the RHS of (3.7). Consequently, numerical solutions for the next time step are obtained. This completes the overall steps of the predictor-corrector method with the TVD version of the WAF method.


## 3.2.1. Wet-Bed Approximation

Assuming that both the left and right waves are rarefaction waves, the wave speed estimatesSL and SR are(3.13)SL=uL-aLqL,SR=uR+aRqR,
where qK with K is L or R given by (3.14)qK={12[(h*+hK)h*hK2]ifh*>hK,1ifh*≤hK,
where h* is an approximation of the exact solution for h in the interface area or usually called as star region. The approximations are given by(3.15)h*=1g[12(aL+aR)+14(uL-uR)]2,u*=12(uL+uR)+aL-aR,
where aK=ghK for K is L or R.


## 3.2.2. Dry-Bed Approximation

The computations of wave speed to the exact dry front are given by(3.16)SL={uR-2aRifhL=0,usualestimateifhL>0,SR={uL+2aLifhR=0,usualestimateifhR>0.
Usual estimate means that the computations of wave speed are obtained from (3.13).


## 3.2.3. TVD Version of WAF Method

The basic idea of the weighted average flux method is described in this subsection. This method provides numerical fluxes at interfaces for the next half-time step. The method is performed by applying cell integration along the interface of the half-time step solution,(3.17)Fi+(1/2)WAF=1Δx∫-Δx/2Δx/2F(Ui+(1/2)(x,Δt2))dx.
In the form of wave structure, the integral (3.17) can be written as the summation(3.18)Fi+(1/2)WAF=∑k=1NC+1wkFi+(1/2)(k),
where the weights wk are(3.19)wk=12(ck-ck-1)withc0=-1,cNC+1=1,(3.20)ck=ΔtSkΔx
is the Courant number of wave k.For this method,Sk is the speed of wave k and NC is the number of conservation laws, or the number of waves in the solution of the Riemann problem.Substituting the weights (3.19) into (3.18), an alternative expression can be obtained,(3.21)Fi+1/2WAF=12(Fi+Fi+1)-12∑k=1NCckΔFi+1/2(k),
with (3.22)ΔFi+1/2(k)=Fi+1/2(k+1)-Fi+1/2(k)
representing the flux jump across the wave k.For a one-dimensional flow, we haveNC=2. So, we need to know Fi+(1/2)(1), Fi+(1/2)(2), and Fi+(1/2)(3) in order to compute weighted average flux in (3.21). Left and right areas can be represented by k=1 and k=3, respectively. We have already known the values of conservative variables in these areas from the initial conditions, except that in the interface area which corresponds to k=2. In this work, the computation by (3.12) is applied to calculate the numerical flux at the interface region.Following [1], a high-resolution method can be obtained by enforcing a total variation diminishing (TVD) constraint to the weighted numerical flux at the interface. The resulting TVD version of the WAF method can be expressed by(3.23)Fi+1/2TVD=12(Fi+Fi+1)-12∑k=1Nsgn⁡(ck)ϕi+1/2(k)ΔFi+1/2(k),
where ϕi+1/2(k) is a WAF limiter function. There are several choices of this limiter function. In this work, we apply the minmod limiter function(3.24)ϕi+1/2(k)=ϕi+1/2(r(k),|ck|)={1,r(k)≤0,1-(1-|ck|)r(k),0≤r(k)≤1,|ck|,r(k)≥1.
The flow parameter r(k) is (3.25)r(k)={Δqi-1/2(k)Δqi+1/2(k)ifck≥0,Δqi+3/2(k)Δqi+1/2(k)ifck<0,
where q is a suitable variable. For shallow water flows, q is set to be h. Note that the value of Δqi-1/2(k) denotes the jump in q across wave k in the Riemann problem with data (Ui-1n,Uin). Hence, the approximation of flux at interface is finally obtained by (3.23). This completes the weighted average flux method for computing flux at interfaces for each cell volume.


## 3.2.4. Calculation of Time Step

For any celli with length of Δxi, time step is defined by(3.26)Δti=CnΔxi(Si+1/2L-Si-1/2R),
where Si+1/2L and Si-1/2R are the left and right wave speed in the Riemann problem for cell i, and Cn is the Courant number coefficient.For easy programming in the case of uniform cell length, the time step can be defined by(3.27)Δt=CnΔxSmax⁡,
where Smax=max{Si+1/2L-Si-1/2R}, for all i.


## 3.2.5. Time Integration in the Corrector Step

Numerical solutions for the next time step can be obtained by using (3.7) as follows. The numerical fluxes at interfaces are calculated by the flux slope limiter (3.23). The source term S involving the bed slope and friction effects is actually the function of the conservative variables h and u. We have already known the values of h and u at the half-time step from the predictor calculations. Hence, we can calculate all terms on the RHS of (3.7). Consequently, numerical solutions for the next time step are obtained. This completes the overall steps of the predictor-corrector method with the TVD version of the WAF method.


## 4. Numerical Results

In this section, the presented method is applied to solve some benchmark problems in both steady and unsteady flows. The accuracy of numerical solutions are demonstrated in various test cases with and without friction effects, as well as the problems of wet and dry beds. In most of all our numerical investigations, we setCn=0.95 for calculating time step size. This value is relatively large comparing with unity that shows the efficiency of the present method.


### 4.1. Steady Flows over a Bump

A steady one-dimensional flow in a channel 25-m-long with a bump at the bottom is defined by(4.1)zb(x)={0.2-0.05(x-10)2,if8<x<12,0,otherwise.This flow is a frictionless test case. The steady flow can be subcritical, trancritical, or supercritical flow with or without a steady shock depending on the initial and boundary conditions. Some analytical solutions in various test cases are proposed by Goutal and Maurel [10]. In our investigations, the simulations run initially by appropriate data until the numerical solutions dramatically converge to some expected solutions.


#### 4.1.1. Transcritical Flow with a Shock

In this test case, the upstream boundary is specified by a discharge per unit width ofq=0.18 m2/s, while the downstream boundary is imposed by a fixed water depth h=0.33 m. The flow domain is discretized by 1,000 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 2 which shows the accuracy of the proposed method.Figure 2
Steady transcritical flow over a bump with a shock.


#### 4.1.2. Transcritical Flow without a Shock

The upstream boundary condition is imposed by a discharge per unit width ofq=1.53 m2/s. No boundary condition is specified at the downstream flow. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 3. They are in good agreement.Figure 3
Steady transcritical flow over a bump without a shock.


#### 4.1.3. Subcritical Flow

A discharge per unit width ofq=4.42 m2/s is imposed at the upstream boundary, while a water depth h=2 m is specified at the downstream boundary. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 4 which again shows very good agreement.Figure 4
Steady subcritical flow over a bump.


### 4.2. Steady Flow over a Variable Bed

A 1000 m long rectangular channel has a discharge of 20m3/s. The flow is subcritical at inflow and is subcritical at outflow. The outflow water depth is fixed to be 0.748409 m. The Manning coefficient for the channel is 0.03. The bed slope is given by(4.2)S0(x)=(1-4gŷ(x)3)ŷ′(x)+0.36(2ŷ(x)+10)4/3(10ŷ(x))10/3,
where (4.3)ŷ(x)=(4g)1/3[1+12exp⁡(-16(x1000-12)2)].The analytical solution to this problem is given byh≡ŷ proposed previously by MacDonald et al. [11]. The flow domain is discretized by 1000 grid cells. The comparison of numerical result and analytical solution is shown in Figure 5. This shows the accuracy of the presented solution. It also reveals that the present method can be used to simulate steady flow in case of variable bed with friction roughness.Figure 5
Steady flow over a variable bed with friction.


### 4.3. Quasistationary Flow

A quasistationary test case was proposed previously by LeVeque [12]. This test case is used to demonstrate the ability of the proposed numerical scheme for computing some small perturbations of initial wave traveling along the channel. The computational domain is 0<x<1. The bed profile is given by(4.4)zb(x)={0.25(cos⁡(π(x-0.5)0.1)+1),if|x-0.5|<0.1,0,otherwise.
The initial conditions are the stationary solution with u=0 and(4.5)h(x,0)={1+ϵ,if0.1<x<0.2,1,otherwise.
where ϵ are some small perturbations.At the early stage, the initial disturbance splits into two waves propagating on the left and right at the characteristic speeds±gh. Many numerical methods may encounter the difficulties of capturing the correct wave speeds of disturbances. The comparisons of the presented solutions with those of LeVeque [12] at t=0.7 for ϵ=0.2 and ϵ=0.01 are shown in Figures 6 and 7. They are in good agreement with the same wave speeds. They also show that the presented scheme can provide accurate solutions in comparing with those are obtained by the high-resolution Godunov method based on balancing the source terms and flux gradients [12].Figure 6
Quasistationary flow over a bump forϵ=0.2 at t=0.7.Figure 7
Quasistationary flow over a bump forϵ=0.01 at t=0.7.


### 4.4. Unsteady Dam-Break Flow in Adverse Slope Channel

The laboratory experiment to induce shock formation of unsteady dam-break flow is proposed by Aureli et al. [13]. This experiment shows reverse flow on wetting and drying conditions. The experiment is conducted in a rectangular channel with length of 7 m, width of 1 m, and Manning n coefficients of 0.025. The initial water depth is 0.292 m in the reservoir and a dry bed in the tailwater. A zero-discharge is specified at the upstream boundary, while a zero-gradient is imposed at the downstream boundary. The bed topography is given by(4.6)zb(x)={0,if0≤x<xc,0.1(x-xc),ifxc≤x≤7,
in which xc=3.5 m.The flow domain is discretized by 700 grid cells. The Courant number is 0.95. The comparisons between numerical results and measured water depth profiles at the four different locations (1.4, 2.25, 3.4, 4.5 m) are shown in Figures8, 9, 10, and 11. These results show very good agreements indicating that the presented scheme can capture not only shock and reverse flows but also wetting and drying wave fronts in transient case.Figure 8
Dam break flow in adverse sloping channel atx=1.4 m.Figure 9
Dam break flow in adverse sloping channel atx=2.25 m.Figure 10
Dam break flow in adverse sloping channel atx=3.4 m.Figure 11
Dam break flow in adverse sloping channel atx=4.5 m.


## 4.1. Steady Flows over a Bump

A steady one-dimensional flow in a channel 25-m-long with a bump at the bottom is defined by(4.1)zb(x)={0.2-0.05(x-10)2,if8<x<12,0,otherwise.This flow is a frictionless test case. The steady flow can be subcritical, trancritical, or supercritical flow with or without a steady shock depending on the initial and boundary conditions. Some analytical solutions in various test cases are proposed by Goutal and Maurel [10]. In our investigations, the simulations run initially by appropriate data until the numerical solutions dramatically converge to some expected solutions.


### 4.1.1. Transcritical Flow with a Shock

In this test case, the upstream boundary is specified by a discharge per unit width ofq=0.18 m2/s, while the downstream boundary is imposed by a fixed water depth h=0.33 m. The flow domain is discretized by 1,000 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 2 which shows the accuracy of the proposed method.Figure 2
Steady transcritical flow over a bump with a shock.


### 4.1.2. Transcritical Flow without a Shock

The upstream boundary condition is imposed by a discharge per unit width ofq=1.53 m2/s. No boundary condition is specified at the downstream flow. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 3. They are in good agreement.Figure 3
Steady transcritical flow over a bump without a shock.


### 4.1.3. Subcritical Flow

A discharge per unit width ofq=4.42 m2/s is imposed at the upstream boundary, while a water depth h=2 m is specified at the downstream boundary. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 4 which again shows very good agreement.Figure 4
Steady subcritical flow over a bump.


## 4.1.1. Transcritical Flow with a Shock

In this test case, the upstream boundary is specified by a discharge per unit width ofq=0.18 m2/s, while the downstream boundary is imposed by a fixed water depth h=0.33 m. The flow domain is discretized by 1,000 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 2 which shows the accuracy of the proposed method.Figure 2
Steady transcritical flow over a bump with a shock.


## 4.1.2. Transcritical Flow without a Shock

The upstream boundary condition is imposed by a discharge per unit width ofq=1.53 m2/s. No boundary condition is specified at the downstream flow. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 3. They are in good agreement.Figure 3
Steady transcritical flow over a bump without a shock.


## 4.1.3. Subcritical Flow

A discharge per unit width ofq=4.42 m2/s is imposed at the upstream boundary, while a water depth h=2 m is specified at the downstream boundary. The flow domain is discretized by 100 grid cells. The free surface profiles for analytical and numerical solutions are shown in Figure 4 which again shows very good agreement.Figure 4
Steady subcritical flow over a bump.


## 4.2. Steady Flow over a Variable Bed

A 1000 m long rectangular channel has a discharge of 20m3/s. The flow is subcritical at inflow and is subcritical at outflow. The outflow water depth is fixed to be 0.748409 m. The Manning coefficient for the channel is 0.03. The bed slope is given by(4.2)S0(x)=(1-4gŷ(x)3)ŷ′(x)+0.36(2ŷ(x)+10)4/3(10ŷ(x))10/3,
where (4.3)ŷ(x)=(4g)1/3[1+12exp⁡(-16(x1000-12)2)].The analytical solution to this problem is given byh≡ŷ proposed previously by MacDonald et al. [11]. The flow domain is discretized by 1000 grid cells. The comparison of numerical result and analytical solution is shown in Figure 5. This shows the accuracy of the presented solution. It also reveals that the present method can be used to simulate steady flow in case of variable bed with friction roughness.Figure 5
Steady flow over a variable bed with friction.


## 4.3. Quasistationary Flow

A quasistationary test case was proposed previously by LeVeque [12]. This test case is used to demonstrate the ability of the proposed numerical scheme for computing some small perturbations of initial wave traveling along the channel. The computational domain is 0<x<1. The bed profile is given by(4.4)zb(x)={0.25(cos⁡(π(x-0.5)0.1)+1),if|x-0.5|<0.1,0,otherwise.
The initial conditions are the stationary solution with u=0 and(4.5)h(x,0)={1+ϵ,if0.1<x<0.2,1,otherwise.
where ϵ are some small perturbations.At the early stage, the initial disturbance splits into two waves propagating on the left and right at the characteristic speeds±gh. Many numerical methods may encounter the difficulties of capturing the correct wave speeds of disturbances. The comparisons of the presented solutions with those of LeVeque [12] at t=0.7 for ϵ=0.2 and ϵ=0.01 are shown in Figures 6 and 7. They are in good agreement with the same wave speeds. They also show that the presented scheme can provide accurate solutions in comparing with those are obtained by the high-resolution Godunov method based on balancing the source terms and flux gradients [12].Figure 6
Quasistationary flow over a bump forϵ=0.2 at t=0.7.Figure 7
Quasistationary flow over a bump forϵ=0.01 at t=0.7.


## 4.4. Unsteady Dam-Break Flow in Adverse Slope Channel

The laboratory experiment to induce shock formation of unsteady dam-break flow is proposed by Aureli et al. [13]. This experiment shows reverse flow on wetting and drying conditions. The experiment is conducted in a rectangular channel with length of 7 m, width of 1 m, and Manning n coefficients of 0.025. The initial water depth is 0.292 m in the reservoir and a dry bed in the tailwater. A zero-discharge is specified at the upstream boundary, while a zero-gradient is imposed at the downstream boundary. The bed topography is given by(4.6)zb(x)={0,if0≤x<xc,0.1(x-xc),ifxc≤x≤7,
in which xc=3.5 m.The flow domain is discretized by 700 grid cells. The Courant number is 0.95. The comparisons between numerical results and measured water depth profiles at the four different locations (1.4, 2.25, 3.4, 4.5 m) are shown in Figures8, 9, 10, and 11. These results show very good agreements indicating that the presented scheme can capture not only shock and reverse flows but also wetting and drying wave fronts in transient case.Figure 8
Dam break flow in adverse sloping channel atx=1.4 m.Figure 9
Dam break flow in adverse sloping channel atx=2.25 m.Figure 10
Dam break flow in adverse sloping channel atx=3.4 m.Figure 11
Dam break flow in adverse sloping channel atx=4.5 m.


## 5. Conclusions

A modified predictor-corrector method for solving the one-dimensional shallow water equations with and without source terms is presented in this paper. This is a modified approach that combines two important concepts which are the depth gradient method (DGM) and the weighted average flux (WAF) method. The DGM is used to approximate the conservative variables at the interfaces in each cell volume for the half-time step. This is the predictor step. So, the source terms can be calculated directly at the same half-time step. By using the WAF method, fluxes at interfaces can be approximated and these values correspond to the half-time step as well. Combining these information together, we can march the numerical solutions in time to reach the next time step of calculations. This step is called the corrector method. For fixing a Courant number value, time step size can be obtained. We apply the presented scheme to solve both steady and unsteady problems. The considering steady problems are flows over a bump with and without shock formation behind the obstacle. It is found that our numerical results are in good agreement with the analytical results in all test cases. For unsteady problem, we have run simulations in order to compare with the existing experimental data of dam-break flow. This problem shows reverse flow with shock on wetting and drying conditions. It is found that our numerical solutions are in good agreement with the experimental data at all observed locations. Hence, the presented scheme can be applied to capture shock formation in shallow water flow for both steady and unsteady cases with or without source terms.



---

*Source: 178491-2011-06-12.xml*