In case of two layer-fluid, two propagating wave modes are there. One is at the free surface and another is at interface. At interface propagating wave number is larger than the free surface wave number. In case of linear water wave potential theory, the free surface elevation is much affected by interface elevation in case of shallow water. These two wave numbers are the real roots of a quadratic dispersion relation. There exist two velocity potentials and one can be found with the help of other.

It is worth thinking about the derivation of the two layer shallow water equations. The the basic idea (in the simplest case) is to do depth averaging of the Navier-Stokes equations ignoring viscous and surface tension effects, and assuming a hydrostatic assumption between layers. This results in a system of coupled system of PDES for the heights h_1 and h-2 , and the depth averaged velocities U_1 and U_2. The equations are conditionally hyperbolic, and when the system loses hyperbolicity you get the Kelvin Helmholtz instability. I believe this is the basic idea. Of course if the system is hyperbolic you can have shocks, waves of permanent form , and tsunami-like structures.

Thanks Brian Higgins