📖 Parcels conceptual workflow#

Parcels is a set of Python classes and methods to create customisable particle tracking simulations using gridded output from (ocean) circulation models.

Here, we will explain the most important classes and functions. This overview can be useful to start understanding the different components we use in Parcels, and to structure the code in a simulation script.

A Parcels simulation is generally built up from four different components:

  1. FieldSet. The input dataset of gridded fields (e.g. ocean current velocity, temperature) in which virtual particles are defined.

  2. ParticleSet. The dataset of virtual particles. These always contain t, z, y, and x, for which initial values must be defined. The ParticleSet may also contain other, custom variables.

  3. Kernels. Kernels perform some specific operation on the particles every time step (e.g. advect the particles with the three-dimensional flow; or interpolate the temperature field to the particle location).

  4. Execute. Execute the simulation. The core method which integrates the operations defined in Kernels for a given runtime and timestep, and writes output to a ParticleFile.

We discuss each component in more detail below. The subsections titled “Learn how to” link to more detailed how-to guide notebooks and more detailed explanations of Parcels functionality are included under “Read more about” subsections. The full list of classes and methods is in the API reference. If you want to learn by doing, check out the quickstart tutorial to start creating your first Parcels simulation.

Parcels concepts diagram

Parcels concepts diagram with key classes in blue boxes#

1. FieldSet#

Parcels provides a framework to simulate particles within a set of fields, such as flow velocities and temperature. To start a Parcels simulation we must define this dataset with the parcels.FieldSet class.

The input dataset from which to create a parcels.FieldSet can be an xarray.Dataset with output from a hydrodynamic model or reanalysis. Such a dataset usually contains a number of gridded variables (e.g. "U"), which in Parcels become parcels.Field objects. A set of parcels.Field objects is stored in a parcels.FieldSet in an analoguous way to how xarray.DataArray objects combine to make an xarray.Dataset.

For several common input datasets, such as the Copernicus Marine Service analysis products, Parcels has a specific method to read and parse the data correctly:

ds_fields = xr.open_mfdataset("insert_copernicus_data_files.nc")
fields = {"U": ds_fields["uo"], "V": ds_fields["vo"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

In some cases, we might want to combine fields from different sources in the same parcels.FieldSet, such as ocean currents from one dataset and Stokes drift from another. This is possible in Parcels by creating multiple parcels.FieldSet objects and combining them into a single parcels.FieldSet:

ds_fields2 = xr.open_dataset("insert_stokes_data_files.nc")
fields2 = {"Ustokes": ds_fields2["ustokes"], "Vstokes": ds_fields2["vstokes"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields2)
fieldset += parcels.FieldSet.from_sgrid_conventions(ds_fset, vector_fields={"UVstokes": ["Ustokes", "Vstokes"]})

Grid#

Each parcels.Field is defined on a grid. With Parcels, we can simulate particles in fields on both structured (parcels.XGrid) and unstructured (parcels.UxGrid) grids. The grid is defined by the coordinates of grid cell nodes, edges, and faces. parcels.XGrid objects are based on Xarray Datasets with attached SGRID metadata, while parcels.UxGrid objects are based on uxarray.Grid objects.

The user doesn’t need to manually construct or manage grids - this is done internally by Parcels. When using parcels.FieldSet.from_sgrid_conventions() or parcels.FieldSet.from_ugrid_conventions() the correct grid object is constructed for your data.

📖 Read more about grids

Interpolation#

To find the value of a parcels.Field at any particle location, Parcels interpolates the gridded field. Depending on the variable, grid, and required accuracy, different interpolation methods may be appropriate. Parcels comes with a number of built-in parcels.interpolators.

📖 Read more about interpolation

🖥️ Learn how to use Parcels interpolators

2. ParticleSet#

Once the environment has a parcels.FieldSet object, you can start defining your particles in a parcels.ParticleSet object. This object requires:

  1. The parcels.FieldSet object in which the particles will be released.

  2. The type of parcels.Particle: A default Particle or a custom Particle-type with additional parcels.Variables (see the custom kernel example).

  3. Initial conditions for each parcels.Variable defined in the parcels.Particle, most notably the release coordinates of t, z, y and x.

t = np.array([0])
z = np.array([0])
y = np.array([0])
x = np.array([0])

# Create a ParticleSet
pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, t=t, z=z, y=y, x=x)

🖥️ Learn more about how to create ParticleSets

3. Kernels#

A kernel function is a little snippet of code, which is applied to the particles in the parcels.ParticleSet, for every time step during a simulation. Kernels define the computation or numerical integration done by Parcels, and can represent many processes such as advection, ageing, growth, or simply the sampling of a field.

Advection of a particle by the flow, the change in position \(\mathbf{x}(t) = (x(t), y(t))\) at time \(t\), can be described by the equation:

\[ \begin{aligned} \frac{\text{d}\mathbf{x}(t)}{\text{d}t} = \mathbf{v}(\mathbf{x}(t),t), \end{aligned} \]

where \(\mathbf{v}(\mathbf{x},t) = (u(\mathbf{x},t), v(\mathbf{x},t))\) describes the ocean velocity field at position \(\mathbf{x}\) at time \(t\).

In Parcels, we can write a kernel function which integrates this equation at each timestep particles.dt. To do so, we need the ocean velocity field fieldset.UV at the particles location, and compute the change in position, particles.dx and particles.dy.

@parcels.validate_kernel
def AdvectionEE(particles, fieldset):
    """Advection of particles using Explicit Euler (aka Euler Forward) integration."""
    (u1, v1) = fieldset.UV[particles]
    particles.dx += u1 * particles.dt
    particles.dy += v1 * particles.dt

Note that the (optional) @parcels.validate_kernel decorator can be used to ensure that the kernel function is correctly defined and follows recommendations.

Basic kernels are included in Parcels to compute advection and diffusion. The standard advection kernel is parcels.kernels.AdvectionRK2, a second-order Runge-Kutta integrator of the advection function.

Warning

It is advised not to update the particle coordinates (particles.t, particles.z, particles.y, or particles.x) directly within a Kernel, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use a change in the coordinates: particles.dy, particles.dx and/or particles.dz. Read the kernel loop tutorial to understand why.

We can write custom kernels to add many different types of ‘behaviour’ to the particles. To do so, we write a function with two arguments: particles and fieldset. We can then write any computation as a function of any variables defined in the Particle and any Field defined in the FieldSet. Kernels can then be combined by creating a list of the kernels. The kernels are executed in order:

# Create a custom Particle object with an "age" variable
CustomParticle =  parcels.Particle.add_variable(
    parcels.Variable("age", initial=0)
)

# Create a custom kernel which keeps track of the particle age
def Age(particles, fieldset):
    particles.age += particles.dt

# define all kernels to be executed on particles using an (ordered) list
kernels = [Age, parcels.kernels.AdvectionRK2]

Note

Every Kernel must be a function with the following (and only those) arguments: (particles, fieldset)

Warning

We have to be careful with kernels that sample velocities on “spherical” grids (so with longitude and latitude in degrees). Parcels can automatically convert velocities from m s-1 to degrees s-1, but only when using VectorFields. This guide describes how to use velocities on a “spherical” grid in Parcels.

📖 Read more about the Kernel loop

🖥️ Learn how to write Kernels

4. Execute#

The execution of the simulation is done using the method parcels.ParticleSet.execute(), given the FieldSet, ParticleSet, and Kernels defined in the previous steps. This method requires the following arguments:

  1. The kernels to be executed.

  2. The runtime defining how long the execution loop runs. Alternatively, you may define the endtime at which the execution loop stops.

  3. The timestep dt at which to execute the kernels.

  4. (Optional) The ParticleFile object to write the output to.

dt = np.timedelta64(5, "m")
runtime = np.timedelta64(1, "D")

# Run the simulation
pset.execute(kernels=kernels, dt=dt, runtime=runtime)

Output#

To analyse the particle data generated in the simulation, we need to define a parcels.ParticleFile and add it as an argument to parcels.ParticleSet.execute(). The output will be written in a parquet format, which can be opened as a polars.DataFrame. The dataset will contain the particle data with at least t, z, y and x, for each particle at timesteps defined by the outputdt argument.

There are many ways to analyze particle output, and although we provide a short tutorial to get started, we recommend writing your own analysis code and checking out related Lagrangian analysis projects in our community page.

🖥️ Learn how to run a simulation