1.  On strategies for model control to manage instabilities

All numerical models, and RSM is no exception, have a way of becoming unstable under a certain set of circumstances. Thus, it seems appropriate, at the start, to provide a general discussion on strategies for model control to manage instabilities. A good physically based mathematical model is based on generally accepted partial differential equations describing the relevant physical processes. RSM uses 1-D and 2-D formulations of watershed, channel, reservoir, and groundwater flow, coupling them as appropriate to better represent the physical reality at the chosen level of abstraction.

All numerical models suffer from problems of stability and convergence. Stability is related to roundoff errors; convergence to discretization errors (O'Brien et al., 1950). A model run on a computer of infinite word length would theoretically be free from roundoff errors; therefore, stable. However, such a computer does not exist. The computers in use today typically have a 32-bit word length, that is, each rational number is represented by a collection of 32 zeros (0) and ones (1), with an accuracy of approximately seven (7) significant digits. In practice, however, this accuracy is not enough; in the longer runs, roundoff errors propagate beyond the stated accuracy, eventually rendering the solution unstable.

Convergence, which is akin to accuracy (in the sense of convergence to the analytical solution), is determined by the size of the discretization, i.e., the values of the discrete space and time steps, which are chosen by the person performing the modeling. In theory, the steps should be small enough to reduce the (n+1)th-order errors of an n-th order scheme to insignificant amounts. This is normally obtained by a careful choice of the discrete steps in order to achieve good spatial and temporal resolutions. The temptation may be to choose very small discrete steps; however, generally this is not the answer. Decreasing the discrete steps increases the number of computations required to reach a solution, thereby increasing the chance for round-off errors to propagate, not to mention the increased computer time required to get a solution.

In practice, the control of numerical instability is seen to be a careful balancing act: How to build a scheme that has enough numerical diffusion to handle the high-frequency perturbations that are responsible for the instability, while at the same time making sure that the solution itself is not being substantially affected by the artificially introduced numerical diffusion. This dilemma is at the crux of all numerical modeling.

The laws of mass and momentum conservation, which underpin all physical-process modeling of unsteady flows, may be combined, through appropriate linearization, into a single second-order, convection-diffusion equation (Hayami, 1951). In one extreme, when the diffusion term vanishes, the equation becomes hyperbolic; in the other extreme, when the convection term vanishes, the equation becomes parabolic.

Numerical models of hyperbolic systems are subject to the Courant law, which expresses the ratio of physical celerity (c) to numerical celerity (Δx/Δt), also referred to as the grid ratio. On the other hand, numerical models of parabolic systems are subject to what has sometimes been referred to (for lack of a better name) as the cell Reynolds number law, which expresses the ratio of physical diffusivity (ν) to numerical, or grid, diffusivity [(Δx)²/Δt]. Both Courant and cell Reynolds numbers control the properties of numerical models of unsteady flow; their values should be calculated a priori (Ponce et al., 2001).

The properties of numerical schemes may be analyzed using various tools of advanced mathematics. A time-tested approach uses Fourier analysis to develop amplitude and phase portraits following the pioneering work of Leendeertse (1967). Significant strides along these lines have already been accomplished by SFWMD scientists. Further clarification of various concepts appears to be in order at this juncture.

In hyperbolic systems, an assessment of numerical accuracy (i.e., convergence) focuses on the spatial resolution L /Δx, where L is the predominant wavelength of the perturbation and Δx is the chosen space step. Generally, numerical models of hyperbolic systems are shown to be more accurate when the grid size follows the characteristic lines, i.e., for a Courant number C = 1, wherein the physical celerity c matches the grid ratio Δx/Δt. In theory, selecting a sufficiently high spatial resolution, say, L /Δx ≥ 100 and a Courant number C = 1 should suffice. In practice, however, a certain scheme may lack enough numerical diffusion to confront the high-frequency perturbations that are likely to appear in well-balanced schemes; thus, additional filtering (numerical diffusion) is normally required to render the system workable.

For instance, there is a wealth of accumulated experience on the numerical properties of the well-known Preissmann scheme, wherein stability and convergence are determined by the spatial resolution L /Δx, the Courant number C, and the weighting factor θ (Ponce et. al., 1978). The latter is required to control nonlinear instabilities which tend to plague the computation as the scheme approaches second order. Values of the weighting factor in the range 0.55 ≤ θ ≤ 1 are recommended, with values near the lower limit approaching convergence (to second order) at the expense of stability, and values near the upper limit approaching stability at the expense of convergence.

An excellent example of the use of Fourier analysis in numerical modeling of flood flows is that of the Muskingum-Cunge model, a diffusion wave model that is based on the matching of physical and numerical diffusivities (Cunge, 1969). A review of the amplitude and phase portraits of the Muskingum-Cunge model, including an online calculator, has recently been accomplished by Vuppalapati and Ponce (2016).