Fast lagrangian particle tracking in unstructured ocean model grids

Abstract

Lagrangian particle tracking, based on currents derived from hydrodynamic models, is an important tool in quantifying bio-physical transports in the ocean. Particle tracking in the unstructured grids typically used in coastal regions is computationally slow, limiting the number of particles and ranges of behaviours that can be modeled. Techniques used in a new offline particle tracker “OceanTracker” are shown to be two orders of magnitude faster than those used in an existing ocean particle tracker for unstructured grids when run on a single computer core. More significantly, its computational speed can exceed that achieved when particle tracking on a regular grid. The techniques for unstructured grids make it possible to routinely calculate the trajectories of millions of particles. This large number of particles allows much better estimates of dispersion and transport statistics, particularly when the probability of connection is low but the consequences are significant, e.g. the spread of invasive species. It also enables wider exploration of parameter sensitivity and particles’ bio-physical behaviours to provide more robust results. The speed increases result largely from exploiting history and reuse within the spatial interpolation of the hydrodynamic model’s output. Using multiple computer cores further increased the speed to track a given number of particles.

Appendix.: OceanTracker structure (version 0.1)

To make it more accessible, OceanTracker (OT) is written in the modern language Python. Key computational tasks used Numba.

The OT framework uses an explicit family tree of Python classes to connect major parts of the code, Fig. 11. This allows the user to easily insert different versions of each family member into the tree, customized to their usage. The OceanTrackerSim class sets up the particle tracking, which can be on one or more computer cores. Solver does the time integration to calculate the particle trajectories. Spatial and temporal interpolation is managed by the FieldGroup which contains the hindcast file reader and an appropriate spatial interpolator. The file reader fills memory buffers with multiple time steps of the hindcast, which are then used for particle tracking. At each time step of the Solver, the Mover class returns the water velocity at each particle’s location, plus any modifications to the velocity due to particle behaviours, such as a fall velocity. The Solver integrates this velocity using RK1, RK2, or RK4 numerical integration.

Fig. 11

Structure of parent-child relationships of classes within Lagrangian code OceanTracker v0.1

Data for individual particles are held within a Particle Group, which holds several types of particle properties. Group properties are relevant to the whole set of particles, e.g. the current time. Individual particle properties which derive from spatial interpolation of hindcast fields are held as Field Particle properties, e.g. particle location, velocity, water temperature, and water depth. Finally, Custom properties are those added by a user, which can be calculated from the other properties; for example, the decaying concentration of a compound within each particle, where the decay rate depends on the water temperature experienced by the particle. Finally, the Particle Writer underpins the Particle Group to write output as required by user. This grouping of particle properties into groups allows many aspects of property creation and use to be automated, e.g. requesting a temperature to be loaded from the hindcast file causes it to be added to a list of particles properties which can then be automatically read, interpolated to particle locations and written to output files, without any changes to the code. All hindcast variables are classified in types of 2D/3D, time-dependent/time-independent, vector/scalars, to further aid this automation.

OT’s explicit family tree allows members to easily obtain information from other members, e.g. the Mover requests the interpolated water velocity from its sibling “fields” using text strings (written in Python as v = self.sibling(“fields”).interp(“water_velocity”) ); similarly a custom particle property could request the interpolated water temperature using code self.greataunt(“fields”).interp(“water_temperature”).

All family members are driven by parameters held in Python dictionaries, which can have specified default values. The parameters contain particle release locations, integration time step etc., appropriate to each family member. This was done so that the OT could be entirely driven by a set of parameter dictionaries, making it easier to support web-based requests for calculated trajectories. The parameter dictionary approach also means that, as users modify or extend a family member through code inheritance, previous parameters are acquired by its descendants. The most likely member to be extended are the Mover, in order to model how a user’s bio-physical particle moves in response to its environment, or the Reader to customize it for the user’s file format.

Another feature to make adaption and extensions easier to implement is to use internal property names as text strings (e.g.“water_temp”), which are mapped to file variable names by the user. These internal names are then used as keys of a Python dictionary to access and update variables stored as matrices using high level code, e.g. particle.set_property(“x”, x_new) to update all particle locations with new vectors. These properties hold matrix data of dimensions determined at start-up from the dimensions of the hindcast data or model parameters. This abstraction of particle properties to dictionary keys also enables automation of operations like reading data files and writing output, which simply loop through the set of keys of the property dictionaries.

Some other features of OT include the ability for particles to be stranded/re-floated by the tide and be re-suspended based on critical shear velocity. A basic lateral boundary condition is implemented in this first version of the code. Particles which leave the hydrodynamic model’s domain are flagged as bad. There is an option to return these bad particles to their last good position, where they are allowed to move at the next time step.

The numerical accuracy of OT’s RK4 time integration was tested by releasing particles in a synthetic hindcast with the 90k nodes spaced at 100-m intervals over a 30-km-wide domain. Particles were released at a 10-km radius from the centre of the eddy, where peak flows were 1 m/s. The numerical drift from the starting radius was 0.1 m per month, for particle tracking done at 5-min time steps.