35.1.7 Linear numerical stability for Redi and GM

Numerically realizing isoneutral diffusion along steep isoneutral slopes is difficult partly because of the small vertical to horizontal aspect ratio in the ocean and hence the ocean model grid. As isoneutral slopes steepen, the projection of diffusive fluxes onto the vertical become stronger, pushing up against the limits of the linear stability criteria for the diffusion equation. This issue is relevant for discretizing both the small and full diffusion tensors. In particular, the linear numerical constraint from the diffusion equation, as discussed in Cox (1987) and Griff ies et al. (1998), indicates that an explicit numerical scheme with a leap-frog time step will be stable if the grid CFL number satisfies

$${\vert K^{m n}\vert \Delta \, t \over \Delta \, x_{m} \Delta \, x_{n} } \le { 1 \over 4},$$ (35.19)

where K^mn are the components of the diffusion tensor ${\bf K}$, and there is no sum implied in this expression. This stability constraint is also relevant for the Gent-McWilliams scheme when implemented according to the skew-flux approach. Similar stability constraints hold for the advective-flux approach.

Assuming a geophysically relevant vertical to horizontal aspect ratio for the grid ( $\Delta x / \Delta z \le 1/1000$), the two dimensional horizontal sub-matrix is stable when the diffusion equation in the horizontal is stable. In general, satisfying this stability constraint in the horizontal is trivial and so will not be considered further. Solving the vertical K^zz diagonal piece implicitly, as done by Cox (1987), points to the K^xz and K^yz cross terms as setting the most restrictive constraint. From these terms, the diffusion equation using the fluxes from the full tensor will be linearly stable when, for each grid cell,

$${\vert S_{a}\vert \over 1 + S_{x}^{2} + S_{y}^{2}} \le { \Delta a \, \Delta \, z \over 4 A_{I} \Delta \, t } \equiv \delta,$$ (35.20)

where a is either x or y. The small tensor's stability is determined with the 1 + S_x² + S_y² denominator set to unity. For the present discussion, the small slope approximation is made. Linear stability for the full tensor is further discussed in Griff ies et al. (1998). For the small angle tensor, $\delta$ represents the maximum allowable slope which can be used before some prescription must be employed to ensure numerical stability. For many large-scale ocean model configurations, this slope check parameter is roughly 1/100, thus providing for the self-consistency of the discretization of the small slope tensor. The prescriptions for maintaining stability are discussed in Section 34.1.9.

This analysis is based on the conservative assumption that if all components to the diffusion tensor produce linearly stable diffusion, then the scheme is linearly stable. Although conservative, experience has shown that violation of these constraints can result in unacceptably large numerical inaccuracies. These inaccuracies are of special concern since they make it more difficult to realize the balance $\alpha F^{z}_{I}(\theta) = \beta F^{z}_{I}(s)$, thus exposing the solution to the nonlinear instability discussed in Griff ies et al. (1998).