LAPS v1.0.0: Lagrangian Advection of Particles at Sea, a Matlab program to simulate the displacement of particles in the ocean

We develop a Matlab program named LAPS (Lagrangian Advection of Particles at Sea) to simulate the advection of suspended particles in the global ocean with a minimal user effort to install, set and run the simulations. LAPS uses the 3D sea current velocity fields provided by ECCO2 to track the fate of suspended particles injected in the ocean, at specific places and times, during a period of time. LAPS runs with a short configuration file set by the user and returns the distribution of the particles at the end of the advection. A continuous tracking option is also available to record the complete trajectory of the 5 particles throughout the entire period of advection. The effect of water waves, or Stokes drift, which alter sea surface current velocities, can also be taken into account. The principle and usage of the program is detailed and then applied to three case studies. The first two cases studies are applied to suspended sediment transport. We show how LAPS simulations can be used to investigate the spatio-temporal distribution of fine particles observed by satellites in the upper ocean. We also estimated sediment deposit areas on the seafloor as a function of sediment grain sizes. The third case study simulates the dispersion of 10 microplastic particles during a tropical cyclone, and shows how the Stokes drift, which is significant during storm events, alters the particles trajectories compared to the case where the Stokes drift is neglected.


Introduction
The quantitative analysis of particle transport by ocean currents relies on in-situ and remote observations, and on Lagrangian simulation of particle transport (van Sebille et al., 2018). The latter allows hypothesis testing and statistical analysis when 15 observations are missing. Indeed, remote sensing observations are global and are focused on sea surface. By contrast, in-situ ocean observations, from ships or drifters, can probe deeper water but remain local. Thus, Lagrangian simulation allows to bridge the gap between different spatial scales of observations. For instance, sediments eroded from the continent are advected until they finally settle in deep oceanic basins. The finest particles remain in suspension for a long time, making wide structures that are visible by remote sensing satellites (Martinez 20 et al., 2015). Then for greater depths, the Lagrangian simulation of particle transport allows evaluating where they eventually accumulate and build sediment deposits (Mouyen et al., 2018). The ocean also spreads living organisms over the globe and Lagrangian simulation can explain how biological dissemination toward the most remote areas is possible (Fraser et al., the program puts an emphasis on particles transport at depth, first because ECCO2 provides 3D velocity fields but also because particles can be set to sink (e.g. due to their mass) with respect to the bathymetry used in ECCO2.
In this article, we will first describe the method and the complete procedure to use LAPS. We then show three case studies applied to sediment transport and movement of microplastic particles in the ocean. They shall emphasize the main simulation features of this program. The velocities of the current come from ECCO2 (Menemenlis et al., 2008) for the global oceanic circulation and a reanalysis of the WAVEWATCH III model for the Stokes drift (Rascle and Ardhuin, 2013). These two datasets must be downloaded by 45 the user before running LAPS, or at least ECCO2 if the Stokes drift is not relevant for the purpose of the simulation (below a few meters depth for instance, as we described in a next section.) The main principle of the program is thus to attribute each particle a velocity according to their location in the velocity grids. The initial position of each particle is provided as an input file by the user. The position of each particle is updated every 15 minutes (simulation time) during the entire simulation period, according to the current velocity at the location of the particle. The injection of particle in the ocean can be done either once, 50 at the beginning of the simulation, or every day until a given date.
The ECCO2 velocity fields used in LAPS have a time resolution of 3 days and a spatial resolution of 0.25 • . Dealing with such a large grid size can represent a computational challenge for Lagrangian models. Therefore, LAPS is designed to efficiently manage the ECCO2 grid by only uploading the spatial and temporal data of interest. When the simulation goes beyond three days, it loads a new ECC02 field and continues the advection. To work with smaller velocity field matrices and decrease the 55 computation time, the program loads only a sub-region of the global ECCO2 velocity field. This sub-region fits the extreme by users willing to do so. At the end of simulation, LAPS stores the final positions of particles and, optionally their trajectories, in folders defined by the user in the configuration file. The date of the simulation and a few particle properties, mostly defining how they can sink, are also defined in the configuration, which in total requires 13 simple inputs from the user.

Oceanic currents
The oceanic currents used in LAPS are the 3-dimensions velocity fields obtained from ECCO2 (Estimating the Circulation and 70 Climate of the Ocean 2). These products assimilate data such as sea surface height from Topex/Poseidon, Jason-1 and Jason-2 missions and ocean bottom pressure from the Gravity Recovery and Climate Experiment (GRACE (Tapley et al., 2004)) mission into a global circulation model (MITgcm (Marshall et al., 1997)) through a Green functions approach that allows to efficiently adjust the model parameters, initial conditions, and boundary conditions (Menemenlis et al., 2005b, a). They are available for the 1992-present time period at a time resolution of three days and a spatial resolution of 0.25 degree. This velocity 75 field is provided for up to 50 levels of water depths, ranging in thickness from 10 m near the surface to approximately 450 m at a maximum model depth of 6150 m from the sea surface down to the sea floor.

Wave action -Stokes drift
Sea surface waves can drag particles in the direction of the wave propagation. This process, also known as the Stokes drift, has a significant impact on particle advection at sea, especially during storms (Fraser et al., 2018;Dobler et al., 2019). Stokes drift 80 velocities derived from wave hindcasts based on the WAVEWATCH III® model (WW3, Tolman, 2009) and forced by winds (Rascle and Ardhuin, 2013) provided by Ifremer at a three-hours sampling can be added to the ECCO2 velocities. Adding the Stokes drift velocity provides a more accurate sea surface velocity field, and can be especially relevant to account for storm events on the particle transport (Fig. 1).
The spatial resolution of the Stokes drift velocity fields is 0.5 degree and its time resolution is 3 hours. LAPS provides a 85 function that converts these velocity fields to the same spatial resolution as ECCO2, prior to any advection simulation. Thus, it is only necessary to locate the particles in the ECCO2 field, then its latitude and longitude indices can be re-used for the Stokes drift velocity field. This saves computation time because locating the particles is computationally demanding. The Stokes drift velocity is maximal at the sea surface and rapidly decreases with depth. LAPS approximates the effect of the depth on the Stokes drift velocity, which is given in Equation 1 (Lynch et al., 2015, eq. 10.52) where U w (z) and U w0 are the Stokes drift velocities at depth z and at the sea surface, respectively, U * is the frictional velocity, κ is the von Karman constant and z 0 is the roughness length (0.5 mm ≤ z 0 ≤ 1.5 mm). With the typical approximations U w0 = 0.03W 10 and U * = 0.0012W 10 (Lynch et al., 2015), where W 10 is the wind velocity 10 m above the sea level, U w (z) can be simplified as Equation 2 95 Thus, we can approximate the Stokes drift velocity at depth z from its surface value, the von Karman constant and the roughness length, which are set by default to κ = 0.4 and z 0 = 1 mm, respectively. Fig. 2 shows the rapid decrease of U w (z)/U w0 i.e. the term in brackets in Equation 2, as the depth increases, for three different values of z 0 . For z 0 = 1 mm, the default value set in LAPS, the Stokes drift velocity at 10 m depth decreases by a factor 10 compared to its value at the sea surface, and by a 100 factor 100 at 25 m depth, consistent with the fact that the Stokes drift mainly impact particles located in shallow waters. This is especially relevant for buoyant particles such as micro-plastic debris and very fine sediment particles (1 µm) that sink at a slow rate.
where U ECCO2 (z) is the depth-dependent velocity of ECCO2 and c(z) is the factor that account for the decrease of the Stokes drift surface velocity ( Fig. 2

Settling of particles
Non-buoyant suspended particles, such as sediment particles, will sink at the Stokes settling velocity, expressed in Equation 4 v s = D 2 g ρ p − ρ s 18µ (4) where g is the mean gravitational acceleration, ρ s = 1030 kg m −3 and µ = 1.4 10 −3 Pa s are the seawater density and dynamic viscosity, respectively. This formalism presupposes that the sinking particles are spherical with ρ p and D the particle density and diameter, respectively, which are set by the user. The settling velocity is added to the vertical component of the oceanic currents and is the main factor governing how long a particle can remain in suspension at sea (Table 1).
Microplastic particles are assumed as fully buoyant at the beginning of each advection simulation. Nevertheless, they can 115 start sinking with a given probability after they spend a certain amount of time at sea, about a few tens of days, as a result of degradation and biofouling (Fazey and Ryan, 2016;Liubartseva et al., 2018). Both this age threshold and the probability of sinking have default values that can be modified by the user. The microplastic particles that are old enough to sink are then chosen randomly with respect to that probability. The sinking velocity of these particles is set to 0.016 m s −1 by default following the study of Kaiser et al. (2017) and can also be changed by the user. The diversity of biofouling processes, plastic 120 types and local parameters at sea do not allow a definitive estimation of such values. That is why it is defined as a user parameter that can be easily modified to any other estimated value.
3 Workflow and algorithm

Workflow
The main task for the user is to write a configuration file that sets up the simulation parameters (an example of such configu-125 ration file is provided with the program). Aside from paths and the simulation dates, the user must only specify another 5 or 6 parameters: 1. if the Stokes drift should be enabled 2. if the particles should be tracked at regular time intervals, as opposed to only the final positions Sediment particles will sink at the Stokes settling velocity, which requires to specify the diameter and density of the sediment particle. Choosing micro-plastic debris means that particles either stay at the surface or may sink after a userdefined time at sea, with also user-defined probability and velocity.
See the online documentation and the example configuration file for setting such parameters.
The user should then download the ECCO2 velocity for the zonal (ECCO2, 2021a), meridional (ECCO2, 2021b) and vertical  Table 2. At the end of the simulation, the results are saved in the folder defined by the user in the configuration file. 145

Algorithm
A simplified algorithm is presented in Algorithm 1. It takes ECCO2 or a combination of ECCO2 and Stokes drift (see section 2.3) as well as particles initial positions as inputs. The main while-loop runs until the final time of simulation set by the user is reached. It performs the injection of particles and updates the particle position for each time step of 15 mins, according to current velocity where the particle is located. The latter step thus implies to locate each particle in the current 3-D velocity field, which is done using the Matlab function histc.
Algorithm 1 Algorithm advection. T is the current time in the simulation, T F is the final time of simulation, T SI is the final time of particle injection. dT is the time increment of the simulation (15 min).
Require: Configuration file, ECCO2 3-D velocity field and initial locations of the particles The validity of the simulation relies on the fact that each particle is properly located in the 3-D velocity field and that it moves according to that local current velocity at its position. In addition, since the velocity field changes at regular time interval, the program must track the simulation time and load a new velocity field whenever the end of such a time interval is reached. We validate that these constraints are respected by creating synthetic ECCO2 and Stokes drift velocity fields, with simple velocity 155 components: 1 m s −1 in pure east, west south and north directions split on four equal longitudinal sections. One particle is set in each of this section and the simulation is run over 5 days of transport, a time long enough to necessitate loading a new velocity field (ECCO2 has a 3-days resolution). We ensured that the simulated particles trajectories and distances agree with the expected results. The synthetic velocity fields, the simulation results and their analysis are archived together with the program LAPS. 10 times smaller than at the surface according to Fig. 2. Fig. 3    One limitation in this comparison is that, because of its size (about 40 cm diameter) the buoy of the drifter, even drogued, is much more sensitive to winds than a particle floating just at the surface or below. Nevertheless, this comparison helps to provide an uncertainty of the particle positions. The standard deviation of the residuals is 1060 km for the East component

Spatiotemporal distribution of suspended material in the ocean
About 19 Gt yr −1 of suspended sediment are discharged by rivers into the oceans (Milliman and Farnsworth, 2011) and take part in various processes, such as building delta fans or fuelling the nutrient cycle that sustains marine life (Darby et al., 2016).
Therefore, simulating the oceanic transport of sediment discharged by rivers can provide useful constraints on the spatio-190 temporal distribution of sediment, which is difficult to assess in-situ due to the wide areas and depths ranges of the ocean.
In this section, we compare the spatio-temporal distribution of fine particles ( The simulation cannot account for all the processes and the variety of material that compose the gelbstof. In addition, the actual particles input leading to MODIS-Aqua observations covers larger areas and time periods than the one we simulated in LAPS. Rather we intend to compare the spatiotemporal pattern of slowly sinking particles out of the basin with one of the world largest sediment discharge (1.43 Gt yr −1 , (Milliman and Farnsworth, 2011)). Indeed, we observe common features: 1. Particles are transported eastward in January and April, and north-westward In July and October 210 2. South of the Amazon estuary, we note a relative increase of particle content in April and relative decrease in October.
Conversely, LAPS can also be used to highlight inconsistencies between observations and simulations in order to better identify sediment sources and grain size. In addition, optical observations are restricted to the sea surface and the transport of sediment that sank deeper shall take advantage from such simulations. LAPS can thus be used to investigate sediment transport for various grain size distributions. Indeed, as already shown in Table 1, sediment grain sizes control their settling velocity and, 215 therefore, significantly alter how sediment distribute in the water column.

Deposition zones of suspended sediment
Here we simulate the area of deposition of the suspended sediment delivered by the rivers Magdalena (Colombia) and Godavari (India), which have significant sediment loads of 140 and 170 Mt yr −1 , respectively. We use the "Sediment mode" of LAPS to evaluate how the sediment deposition area might be controlled by the grain size of the sediment, ranging from 25 to 95 µm, 220 with 10 µm increase. Indeed, for a given density, the size of grain is responsible for most of the variation of particle settling velocities. The results in Fig. 6 show that, overall, fine grain size extend further than coarser ones, simply because they need more time to settle and can thus travel over longer distances. However, the shape of the sedimentation zone is a rather complex combination of the sea-floor depth (as given in ECCO2), the grain size and the sea current velocities, from which we expect to 230 obtain realistic estimate of contemporary sedimentation zones.

Ocean plastic pollution
The pollution of the marine environment with macro and microplastic debris is rising and has become a global problem (Borrelle et al., 2020). It is estimated that more than 5 trillion plastic particles weighing 250,000 tons float on the ocean surface (Eriksen et al., 2014)  ). Therefore, modelling the transport and sinking of plastic particles is of great interest as it allows to estimate the location of pollution hotspots at the sea surface as well as in the deep sea. Additionally, in recent years, it has become more and more evident that rivers play a major role in polluting the oceans with plastic (Jambeck, 2015;Lebreton et al., 2017). It was estimated that every year 1.15-2.41 million tons of plastic waste are carried by rivers into the ocean . The fate of riverine plastic waste is not well known. Modelling can help to identify the transport ways and the locations where plastic 245 might sink to the ocean floor or is beached on the coast. These results can in turn help to plan mitigation or cleaning efforts. We use LAPS to run two different scenarios, scenario 1 shows the release of plastic particles by rivers in East-Asia and in scenario 2 we demonstrate the impact of a typhoon on the particle trajectory.

Scenario 1
In this scenario, we use the Tamsui River in Taiwan, and the Pearl River and Yangtze River in China as examples for riverine 250 plastic output. All three rivers carry significant amounts of micro and macroplastic waste into the ocean (Lin et al., 2018;Wong et al., 2020;Schneider et al., 2021;Zhang et al., 2021) and belong to the most polluting rivers worldwide Lebreton et al., 2017).  Karkanorachaki et al. (2021).

260
The resulting trajectories are shown in Fig. 7a, b and c, respectively, and it can be seen that Stokes drift has a significant influence on the particle trajectories. Especially for the particles that were released by the Pearl River. Without Stokes drift, the majority of particles traveled from the Pearl River mouth in a north-eastern direction. Some particles beached at the north coast of Taiwan, others continued further northwards. However, when applying Stokes drift in the model, all particles released by the Pearl River stay in a gyre close to the river mouth. For particles that were released by the Tamsui River and the  China. Ten particles were placed randomly in a small area in the eastern part of the Bashi Channel between Taiwan and the Philippines. Injection of particles was done at two different times. First injection was on May 1st 2019 and particles traveled for three months through the undisturbed western Pacific (Fig. 8a). Second injection was on July 1st 2019 and particles crossed the path of typhoon Lekima (Fig. 8b). For both time periods, the trajectories with and without Stokes drift effect were modelled.
For the time period May to July, it can be seen that Stokes drift has a significant impact on the location of the endpoints. The

290
Two features are considered for future development. The first one is the ability for storms and currents to re-mobilize sediment particles that settled at shallow depths (Stastna and Lamb, 2008). At the moment, a workaround is to update the particles input file with points located in areas prone to sediment re-suspension, and to keep the injection running over the requested time.
This will act as if sediment particles were re-suspended from these locations.
The second feature is the effect of tidal currents, which can also impact the sediment redistribution in coastal areas. The 295 significance of this effect depends on the area considered and on the timescale at which the sediment transport is considered, since the strongest tidal currents have a 6-to 12 hours period (to be compared with the 3-days resolution of ECCO2 velocity fields used for LAPS). The last FES2014b products (Carrère et al., 2015) should be suitable for that purpose.

Conclusion
We described a Matlab program for Lagrangian advection that allows to simulate the 3D displacement of particles in the 300 ocean. The utilisation of LAPS is made simple by only requesting to download publicly available current velocity data and setting a short configuration file for basic parameters (such as the time of the simulation, the paths to working directories, particle sinking or particle tracking). By comparing LAPS results with actual drifters trajectories over one year of transport, we evaluate particle position uncertainties to 1060 km for the East component and 560 km for the North component. We run example of applications that show how LAPS can handle particle grain size and Stokes drift in order to simulate particles 305 trajectories in various cases. We suggest that LAPS is suitable for investigating sediment transport, plastic littering and other suspended particles.
Code and data availability. The exact version of the model used to produce the results used in this paper is archived on Zenodo (Mouyen, 2021) under the MIT licence, as are input data and scripts to run the model and produce the plots for all the simulations presented in this paper .