Assessment of numerical integration schemes and boundary conditions for Lagrangian tracking of marine debris

Abstract

A growing environmental issue today is the amount plastic pollution in the ocean. To investigate how it might affect marine life, an important aspect is to study how plastics are transported by ocean currents and where it ends up. This can be done by use of Lagrangian simulation of plastic particles in computer ocean models. To improve the quality of the output from such simulations the number of simulated particles is continuously increased, and grid resolutions are refined. This puts increasing pressure on the simulation programs to perform the computations as efficiently as possible.

In the main part of this project the numerical time stepping schemes used in a Lagrangian ocean simulator, and their influence on the efficiency and accuracy of the simulations, were studied. The currently used 4th order Runge-Kutta method (RK4) was compared to the Runge-Kutta-Fehlberg method (RKF45) to assess the potential performance gain by switching to a more advanced scheme. The Explicit Euler method was included in initial tests to illustrate the benefit of using a higher order method. Error testing was performed on an idealized test case based on the Stommel equations, and time testing was performed on model data from the Agulhas region off the coast of South Africa. RKF45 was found to produce less accurate results than RK4 in very similar computation time. However, it used less function evaluations, which means it might still be useful in the future if more advanced interpolation schemes are introduced.

Another numerical aspect was targeted in the secondary part of the project, namely the effect of the boundary conditions between land and sea on the particle trajectories. Simple conditions, such as the currently used no-normal flow partial slip condition, include slowing particles as they approach land to prevent them from accidentally beaching. To allow them to flow freely even along shorelines, a free-slip condition was introduced and evaluated. The two different methods were compared using the Agulhas data. With partial slip, many particles beached along the shorelines, but with free-slip not a single one of the 300 particles used got stuck.

Popular Abstract

The plastic pollution in the ocean is a growing issue for the marine environment. To determine how and where the pollution affects marine life, an important step is to map the spread of this pollution. This project concerns the numerical aspects of a computer program developed to simulate the distribution of plastics in the ocean.

Each year millions of tons of plastics enter our oceans around the world. With time the larger plastic debris is broken down into smaller pieces by the sunlight and grinding movements of the waves. These smaller pieces can easily get ingested by marine animals, and thus harm them. To evaluate the extent of this threat to marine life, it is important to map the spread of plastics in the ocean. This can be done with the help of computer simulations of the ocean currents.

One method of simulating the flow of objects in the ocean is by tracking individual particles, in this case pieces of plastic. In an ocean model, the ocean currents are given as velocity fields that define the currents at discrete points in space. The velocity of a particle is calculated by interpolating the velocities at its nearest grid points, and then the particle is transported by taking discrete steps of a predefined time length. In practice, millions of particles are tracked simultaneously. This way of using particles to simulate transportation in a fluid is known as Lagrangian simulation.

PARCELS is a computer program that is currently being developed to not only track plastics in the oceans, but also other particles such as oil droplets or tuna. In the main part of this project, different methods for stepping a particle forward in time are compared, to evaluate whether PARCELS would benefit from using a more advanced scheme than it currently does. The current method is known as the 4th order Runge-Kutta method (RK4). To take a step of length dt in time, it uses the velocity at not only the starting point, but also at three intermediate points. This produces a 4th order approximation, i.e. it reproduces polynomials up to the 4th degree exactly.

The method that I implemented in this project is known as the Runge-Kutta-Fehlberg method (RKF45). In addition to the 4th order step, it also calculates a 5th order step, and takes the difference between these two as an approximation of the local error. By putting a tolerance on this error the step size can then be varied to take longer steps where the tracked path is nice, and shorter steps only where necessary.

In the secondary part of the project, the boundary conditions that determine how a particle behaves close to land were considered. The currently used conditions are no-normal flow partial slip. This effectively means that the velocity perpendicular to land goes to 0 as the particle approaches the boundary, and velocities along a beach are slowed, but not quite to 0. To allow particles to flow near a beach without getting slowed, I implemented a free-slip condition. By mirroring the closest ocean velocities onto land when the particle is close, the particle can be transported without getting slowed even close to boundaries.

To compare the methods used for time stepping, two different test cases were used. The first one used the Stommel equations to simulate a periodic trajectory in a rectangular region. The velocities along the western edge of this region were very large compared to the rest. This test case was used to get a correlation between the number of steps used and the error produced by RK4 and RKF45.

The second test case used real simulated currents from the Agulhas region, off the coast of South Africa. Here the correlation found in the Stommel test was used in an attempt to produce runs with RK4 and RKF45 that gave rise to similar errors, and compare the computation time of these runs.

The boundary conditions were also evaluated using the Agulhas data. 300 particles were released along the coastline, using first the partial slip condition, and then the free-slip condition. The number of particles getting stuck along the shore in both runs were counted.

In the Stommel test, RKF45 was found to produce more accurate results in far fewer steps than RK4. However, the correlation found did not translate well to the Agulhas region. The timed runs took 934.4 s in RK4 and 935.2 s in RKF45, i.e. almost exactly the same. However, the errors produced by the RKF45 run in the Agulhas region were found to be larger than those produced by RK4, meaning RK4 produced more accurate results in the same computation time.

The part of the time stepping that was expected to be most computationally expensive was the function evaluations, i.e. interpolating the grid points to calculate the velocity of the particle. RKF45 was indeed found to use fewer function evaluations than RK4, but was still slower, which implies that other overhead computation was relatively more expensive in RKF45. However, the current method for calculating velocities is by using bilinear interpolation of the grid points. This method is cheap, but might not be very accurate. Thus as the development of PARCELS continues, it might be switched to a more advanced method which would make function evaluations more expensive, and possibly also make RKF45 perform better than RK4.

The free-slip condition proved to be very successful in the testing in the Agulhas region. With partial slip, 95 particles, i.e. more than a third, got stuck along the shoreline. With free-slip, however, not a single particle got stuck. Theoretically, there are still ways for particles to get stuck even with the free-slip conditions, but this is highly unlikely and will probably not cause problems in simulations that use millions of particles.

To conclude, even though RKF45 used less function evaluations than RK4, they ran at very similar computation times. RKF45 produced larger errors in this test, but if function evaluations are made more advanced and computationally expensive in the future, it might still prove useful. The new free-slip boundary conditions prevented particles from getting stuck along the boundaries between land and ocean, and were thus very successful.

• author: Kronborg, Joel
• supervisor: Eskil Hansen; Erik van Sebille
• organization: Mathematics (Faculty of Engineering)
• course: FMN820 20161
• year: 2016
• type: H2 - Master's Degree (Two Years)
• subject: Mathematics and Statistics
• keywords: Lagrangian tracking, PARCELS, Runge-Kutta, time stepping, boundary conditions
• publication/series: Master's Theses in Mathematical Sciences
• report number: LUTFNA-3037-2016
• ISSN: 1404-6342
• other publication id: 2016:E36
• language: English
• id: 8888535
• date added to LUP: 2016-11-16 10:37:05
• date last changed: 2016-11-16 10:37:05