Handling discontinuities in numerical ODE methods for Lagrangian oceanography

Abstract. In Lagrangian oceanography, numerical methods for Ordinary Differential Equations (ODEs) are used to model particle transport. In many common applications, the velocity field driving the particle transport is provided as output from ocean models, on a discrete grid of points. Hence, the velocity field must be interpolated. Depending on the choice of interpolation, the velocity field or its derivatives may have discontinuities. These discontinuities have implications for the accuracy of the numerical ODE methods employed. We demonstrate that by using information about the location of the discontinuities, we can take these into account, and improve numerical accuracy over standard integration methods that do not take discontinuities into account. The commonly used combination of the fourth-order Runge-Kutta method and linear interpolation of the velocity field, in fact, only yields second-order accuracy with the standard method. By accounting for discontinuities, we can achieve several orders of magnitude better accuracy with the same timestep. The implementation makes use of a combination of known methods from the field of numerical integration of ODEs. The implementation is quite flexible, agnostic to grid layout and order of interpolation, and contributes only modestly to the code complexity. Hence, the proposed technique for handling discontinuities in interpolated velocity fields could easily be adopted to a range of applications where numerical accuracy or efficiency is of importance. As an example where numerical accuracy is important, we run a backtracking case for particles with known initial conditions, and show that the method with discontinuity handling is to a larger degree able to recover the correct initial positions of the particles, compared to standard fourth-order Runge-Kutta.