Algorithmic Improvements To Finding Approximately Neutral Surfaces

Interior oceanic motions occur predominantly along, rather than across, the neutral tangent plane. These planes do not link together to form well-defined surfaces, so oceanographers use approximately neutral surfaces. To date, the most accurate such surface is the omega-surface, but its practical utility was limited because its numerical implementation was slow and sometimes unstable. This work upgrades the speed, robustness, and utility of omega-surfaces. First, we switch from solving an overdetermined matrix problem by minimal least squares, to solving an exactly determined matrix problem, obtained either by the normal equations (multiplication by the matrix's transpose) or by discretizing Poisson's equation derived from the original optimization problem by the calculus of variations. This reduces the computational complexity, roughly from O(N1.6) to O(N1.2), where N is the number of grid points in the surface. Second, we update the surface's vertical position by solving a nonlinear equation in each water column, rather than assuming the stratification is vertically uniform. This reduces the number of iterations required for convergence by an order of magnitude and eliminates the need for a damping factor that stabilized the original software. Additionally, we add "wetting" capacity, whereby incrops and outcrops are reincorporated into the surface should they become neutrally linked as iterations proceed. The new algorithm computes an omega-surface in a 1,440 by 720 gridded ocean in roughly 15 s, down from roughly 11 h for the original software. We also provide two simple methods to label an omega-surface with a (neutral) density value.