FINITE-DIFFERENCE SCHEME

Although sophisticated, higher-order finite-difference schemes exist for the solution of hyperbolic conservation laws (Centrella and Wilson, 1984; Hawley et al., 1984; Harten et al., 1987), a first-order upwind finite-difference scheme appears generally accurate enough for geophysical applications. We implemented a first-order method described by Engquist and Osher (1980).

The first-order forward and backward finite differences in space are, respectively,

$$\begin{array} {ll} \Delta_+\u\ &=\ \u_{j+1} - \u_j; \\ \ \\ \Delta_-\u &=\ \u_j - \u_{j-1}, \end{array}$$ (8)

where $\u_j$ is the discrete representation of $\u(\th)$ on the grid $\th_j = j\Delta \th$.

The basic upwind scheme for equation (7) is  

$${{\u^{n+1}_j-\u^n_j} \over {\Delta \r}} \ =\ {1\over{\Delta ... ...'(\u^n_j)\ {\rm and}\ F'(\u^n_{j-1}) \gt 0, \end{array}\right.$$ (9)

with n the discrete-radius index of the grid, $\r^n = n\Delta \r$.The direction of the finite-difference operator depends on $F'(\u)$:when the flux increases from left to right, flow is to the left; when it decreases, flow is to the right (see Figure 1). Note that an equivalent scheme can be set up in Cartesian coordinates, where $\r$ is replaced by z, and $\th$ by x.

This specification is incomplete, of course. The various upwind schemes differ in the way in which the intermediate cases are handled; i.e. when the sign of $F'(\u)$ changes amongst the three points of the stencil. The scheme of Engquist and Osher handles the intermediate cases via the formula:

$$\u^{n+1}_j= \u^n_j\ +\ {{\Delta \r}\over{\Delta \th}}\ \Bigl(\Delta_+F_-(\u^n_j)\ +\ \Delta_-F_+(\u^n_j) \Bigr).$$ (10)

Here  

$$\begin{array} {ll} F_+(u)\ &=\ F({\rm max}(\u,{\overline{\u... ... \\ F_-(u) \ &=\ F({\rm min}(\u,{\overline{\u}})), \end{array}$$ (11)

with ${\overline{\u}}$ the stagnation point: $F'({\overline{\u}}) = 0$. ${\overline{\u}}= 0$ for the function F under consideration here. Engquist and Osher's scheme reduces to the standard first-order scheme (9) when $\u$ has a consistent sign, and is stable near gradient discontinuities.

Note that the specification given above is still not complete: the slowness is position-dependent, therefore so is the flux. We have suppressed this dependence in the notation, but it must be respected. Some obvious ways of evaluating the flux lead to inconsistent schemes. In general it seems sufficient to verify that the elementary upwind differences above are recovered when the sign of $\u$ is consistent. For example, taking the slowness locally constant at its value at $\u^n_j$ seems to work.

Because an explicit scheme is used, each step (fixed by the input slowness grid) may require several partial steps to satisfy the Courant-Friedrichs-Lewy stability condition. The appropriate partial steps are determined adaptively, and the full step is always taken if stably possible. Note the contrast with Vidale's code, which uses an implicit upwind scheme.