The model uses the analytical solution to the linearized 3D internal wave equation to advect particles. This is actually more computationally time consuming than simply solving the equations directly, but it eliminates numerical diffusivity and allows for a restricted subdomain in the vertical.
The first example shows a single mode solution, in this case the first vertical (baroclinic) mode, and first (k=1, l=1) horizontal modes. The stratification is taken to be constant.
This second example shows a full spectrum of waves meant to mimic typical ocean conditions. Because the domain and stratification aren’t particularly realistic, the spectrum is set to only 0.1 of the Garrett-Munk reference level.
The interface at the top and bottom of the ocean doesn’t move, while the density interfaces in the interior wiggle around quite a bit. As a consequence, the particles at the surface only move horizontally (as might be expected), while the particles in the interior move in all three dimensions.
Code for the above two examples is available from GitHub.