Numerical solution of the two-layer shallow water equations with bottom topography

We present a simple, robust numerical method for solving the two-layer shallow water equations with arbitrary bottom topography. Using the technique of operator splitting, we write the equations as a pair of hyperbolic systems with readily computed characteristics, and apply third-order-upwind differences to the resulting wave equations. To prevent the thickness of either layer from vanishing, we modify the dynamics, inserting an artificial form of potential energy that becomes very large as the layer becomes very thin. Compared to high-order Riemann schemes with flux or slope limiters, our method is formally more accurate, probably less dissipative, and certainly more efficient. However, because we do not exactly conserve momentum and mass, bores move at the wrong speed unless we add explicit, momentum-conserving viscosity. Numerical solutions demonstrate the accuracy and stability of the method. Solutions corresponding to two-layer, wind-driven ocean flow require no explicit viscosity or hyperviscosity of any kind; the implicit hyperdiffusion associated with third-order-upwind differencing effectively absorbs the enstrophy cascade to small scales.