Tutorial on how to 'delay' the start of particle advection

In many applications, it is needed to 'delay' the start of particle advection. For example because particles need to be released at different times throughout an experiment.

The way to achieve this in Parcels is by using a combination of pset.execute() and pset.add() calls inside a for-loop.

This tutorial will show how this can be done. We start with importing the relevant modules.

In [1]:
from parcels import FieldSet, ParticleSet, JITParticle
from parcels import AdvectionRK4
import numpy as np
from datetime import timedelta as delta
from scripts import *

Now import a FieldSet (from the Peninsula example, in this case)

In [2]:
fieldset = FieldSet.from_nemo('Peninsula_data/peninsula', allow_time_extrapolation = True)

Define an empty ParticleSet

In [3]:
pset = ParticleSet(fieldset=fieldset, lon=[], lat=[], pclass=JITParticle)

And now also define some other variables

In [4]:
npart = 10  # number of particles to be released
delaytime = delta(hours=1)  # time between each particle release
endtime = delta(hours=24)  # total time of the experiment

We use the same start locations as in examples/example_peninsula.py

In [5]:
x = 3. * (1. / 1.852 / 60)  # 3 km offset from boundary
y = (fieldset.U.lat[0] + x, fieldset.U.lat[-1] - x)  # latitude range, including offsets
lat = np.linspace(y[0], y[1], npart, dtype=np.float32)

Now we define an output file. Note that, since not each particle will have the same number of timesteps outputted (because some will be released later than others), we need to run with type="indexed". See also H3.5 at http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/aphs03.html

In [6]:
output_file = pset.ParticleFile(name="DelayParticle", type="indexed")

Now we run the combination of pset.add() and pset.execute() calls inside a for-loop.

In [7]:
for t in range(npart):
    pset.add(JITParticle(lon=x, lat=lat[t], fieldset=fieldset))
    pset.execute(AdvectionRK4,
                 runtime=delaytime,  # note we execute only for 1 hour
                 dt=delta(minutes=5),
                 interval=delaytime,
                 starttime=delaytime*t,  # note starttime is an hour later for each execute
                 output_file=output_file)
Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gn/T/parcels-501/8d4583bd4e9b01d643c0fddfd29f1080.so

And then we run the particles for the rest of the time until endtime

In [8]:
pset.execute(AdvectionRK4, runtime=endtime-npart*delaytime,
             starttime=delaytime*npart, dt=delta(minutes=5), interval=delta(hours=1),
             output_file=output_file)

And then finally, we can show a movie of the particles. Note that one particle appears every frame for the first ten frames.

In [9]:
plotTrajectoriesFile('DelayParticle.nc', mode='movie2d_notebook')
Out[9]: