Print ISSN: 0976-2876 Online ISSN: 2250-0138

Available online at: http://www.ijsr.in

Indian Journal of Scientific Research
DOI:10.32606/IJSR.V14.I1.00005

Received: 28-05-2023 Accepted: 29-07-2023 Publication: 31-08-2023

Indian J.Sci.Res. 14 (1): 29-36, 2023 Original Research Article

USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH

INDIAN OCEAN
KARAN CHAKRAVARTHY1
Department of Environmental Sciences, American International School, Chennai, India

ABSTRACT
As oceans get progressively more polluted and garbage patches are discovered across the globe, tracking the
movement of waste across the oceans becomes increasingly important. Modeling this movement makes identifying ocean gyres
(and possible garbage patches) easier and could serve as an educational tool to demonstrate the effects of discarding trash in
the ocean. To create such a model, we represented ocean currents as a vector field, and the fourth-order Runge Kutta method
is used along with three different interpolation methods as bilinear, biquadratic, and bicubic interpolation. Bicubic
interpolation was chosen as the reference solution, and the error associated with bilinear and biquadratic interpolation was
compared. Quantitative results show that there is no clear superiority of the biquadratic model over the bilinear model.
Qualitative results show that the model performs as expected based on studies of ocean currents in the North Indian Ocean.
This paper demonstrates how to implement such a model, a method for evaluation, and recommendations for further study.

KEYWORDS: Ocean, Biquadratic Model, Bilinear

Our oceans have progressively become more and DATASET AND REPRESENTATION
more polluted with plastics, and much of that plastic
We decided to model the ocean currents as a
comes from rivers. (Ritchie H., 2021) Moreover, with the
two-dimensional vector field and plastics moving through
discovery of the Great Pacific garbage patch, the Indian
the ocean as a particle’s trajectory through this field. The
Ocean garbage patch, and others across the oceans, ocean
vector field is a two-dimensional array of vectors sampled
current analysis has shown that trash entering the oceans
at equally spaced points across the Earth’s surface.
tends to be caught within the ocean gyres and remain at
sea. Studies are also showing that much of the trash We used data from NASA’s Ocean Surface
expelled by rivers wash up on beaches (Egger M., 2023). Current Analyses Real-time (OSCAR) global surface
current database. OSCAR data provides current velocities
While there exists a handful of studies to
for the ocean’s mixed layer (also known as the surface
understand how currents affect the movement of plastics
layer). This topmost layer of the ocean is of almost
across the oceans, most focus on the development of
uniform density and lies above the pycnocline (Figure 1).
theoretical models that track the motion of particles in a
OSCAR data is derived from satellite measurements of
vector field. Part of the challenge comes from the fact
sea surface heights, ocean winds, and sea surface
that there is little ground truth available to validate any
temperatures, and the dataset is updated approximately
model that is developed. Nonetheless, as the research
every five days. For the work described in this paper, the
develops and further empirical evidence is gathered, the
dataset used was from January 1, 2023.
background work done to develop and evaluate these
models will be valuable.
The purpose of this work is to create a model
using numerical methods to track the movement of
plastics across the oceans, with a focus on the North
Indian Ocean. The fourth-order Runge Kutta model was
used to determine the trajectory of a point within a vector
field of ocean currents. Three interpolation methods –
bilinear, biquadratic, and bicubic – were implemented
and evaluated, and the results presented.

Figure 1: Diagram of the ocean layers

________________________________
1
Corresponding author
 CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN

Four fields were extracted from the OSCAR ⃗

dataset: u, v, latitude, longitude. The u and v fields
represent the two components of the vector current: the ⃗
zonal current component (east-west currents along the
lines of latitude) and the meridional current component ⃗
(north-south currents along the lines of longitude)
respectively with units of meters per second. The latitude ⃗
field is a set of Float64 values representing the latitudes
from -80N to +80N divided into 480 values, with each ⃗
value separated evenly by 0.333 degrees. The longitude
⃗
field represents the longitudes from 20E to 420E divided
into 1200 values, with each value similarly separated The final position vector k is created as a
evenly by 0.333 degrees. Together, this data results in a weighted average of the four vectors , , ,
grid of 480 x 1200 with each grid point associated with a according to the following equation.
latitude, longitude, and 〈u, v〉 vector. Each vector
represents the surface currents across approximately 32 ⃗ ⃗ ⃗ ⃗ ⃗
square kilometers.
With the final position vector calculated, the next point in
Instead of using latitudinal and longitudinal the trajectory is given by
values to calculate the trajectory of a particle, the latitude
and longitude values were mapped to a 2D array of ⃗
index values with values ranging from 1 to 480 and y This is depicted pictorially as shown in Figure 2.
values ranging from 1 to 1200. The latitude and longitude
data were overlaid in such a way that =
represented a latitude of 80N and a longitude of 20E
while = represented a latitude of -80N
and a longitude of 420E.

METHODS AND IMPLEMENTATION

Runge Kutta Method
After consideration of various numerical
methods to trace a particle’s path through a vector field, Figure 2: Diagram of weighted vectors in the RK4
we chose the fourth order Runge Kutta algorithm. The algorithm
fourth-order Runge Kutta algorithm is a commonly used
predictor-corrector method and uses a weighted average Bilinear Interpolation
of four nearby vectors to determine the next step in a Often, a vector field in is described as the
particle’s path. gradient of a function so a vector can be
Given a starting point , the first vector ⃗ is determined at any point in the field. However, because
the vector at point . The point is then calculated as ocean currents are represented by a sampling of points
across a vector field, it is necessary to interpolate the
one half of the step size along the vector ⃗ from the
vectors that lie between the discrete samples. To do so,
point . The second vector ⃗ is the vector at , and we considered three methods to interpolate vectors
is then calculated as one half the step size along the between the discrete samples: bilinear interpolation,
vector ⃗ again from the point . The third vector ⃗ is biquadratic interpolation, and bicubic interpolation. These
the vector at , and is calculated as the step size along methods are variations of linear, quadratic, and cubic
the vector ⃗ from point . Finally, the fourth vector ⃗ interpolation applied to a 2D space.
is the vector at point The equations for each of these Bilinear interpolation is the simplest of these
vectors are given below where is the vector at the three and is commonly used in computer graphics and
point and is a constant step size. image processing to estimate values between grid points
⃗ or other applications where the dataset is relatively
uniform. It involves determining the value of an unknown
point based on the values of the 4 neighboring points.

30 Indian J.Sci.Res. 14 (1): 29-36, 2023

 CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN

Assume a point to be interpolated with discrete points. This technique involves fitting a smooth
location lies within a rectilinear cell of four points curve between data points using a piecewise-defined
( , ), ( , ), ( , ), and ( , ) each of which is quadratic polynomial. Unlike bilinear interpolation, the
associated with a vector as shown in Figures 3 and 4. biquadratic interpolation method takes nine data points
surrounding a central point instead of four. This allows
for a greater representation of the data as more data points
are included when determining the interpolation.
With linear interpolation, one point is chosen on
either side of the point to be interpolated. With quadratic
interpolation, three points are considered. However, the
odd number of points leaves two points on one side and
one point on the other, and the interpolation would not be
Figure 3: Visual representation to find percentage centered about the point to be interpolated.
distance from closest corner ( , ) to To address this, we shifted the interpolation by
half a unit. This was done by creating two virtual control
points half way between the first and second points and
the second and third points. The point to be interpolated
is designed to be between these two virtual control points.
The values for each of these virtual control points are
determined by linearly interpolating between the two
corresponding points. When applied to a 2D space, this
results in 4 virtual control points as shown in Figure 5.

Figure 4: Rectilinear cell of four vectors with calculate

vector at point

The vector ⃗ at point is calculated as follows:

⃗ ⃗ ⃗

where

⃗ ⃗ ⃗ Figure 5: Biquadratic grid of points used to

interpolate the data
⃗ ⃗ ⃗
Virtual control points shown in blue and
interpolated points shown in red.
Biquadratic Interpolation
With the linear interpolation included as
While bilinear interpolation can provide a described above, the general equation used for quadratic
suitable representation of the vector field, biquadratic interpolation across three points , , and is given as
interpolation provides a smoother curve between the follows.
( ) ( )

where and are normalized to be between 0

and 1 between the virtual control points.
Once the values for , , and are
When extended to biquadratic interpolation, the determined, the quadratic interpolation is performed once
quadratic function is first interpolated three times along more along the y-axis to determine the value for , the
the x-axis as follows. point to be interpolated.

Indian J.Sci.Res. 14 (1): 29-36, 2023 31

 CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN

Bicubic Interpolation
Bicubic provides an even higher degree of
smoothness between grid points, using a cubic function to
fit a curve between the discrete points. Like bilinear
interpolation, bicubic interpolation uses an even number
of points so the creation of virtual control points is not
needed.
Bicubic interpolation works by taking a group of
16 points in a 4 x 4 grid of data points surrounding the
point to be interpolated. Like both of the other
Figure 6: Grid of data points used for bicubic
interpolation methods, the interpolation is first run across
interpolation. Interpolated coordinate points shown
the x-axis to calculate the points , , , and , and
in red
then run once along the y-axis to determine the value at
. As earlier, and are normalized to be between 0 The equation for cubic interpolation across the x-axis is
and 1. This is illustrated in Figure 6. defined below.

( ) ( )

This equation is applied four times across the distance is first calculated between corresponding
points along the x-axis to determine the values for points points along the two trajectories as follows:
at and .
√ ( )

where P and Q are points in the 2D space.

The values for a given trajectory are then
combined as the square root of the sum of squares. This is
known as the lock-step Euclidean distance and is
The value at is finally determined through described as follows:
one more application of the equation for cubic
interpolation, but this time across the four interpolated √∑
points and along the y-axis.
where M and N are two vectors of points in the
2D space and ranges from 1 to the number of elements
Implementation in the vector. In our case, M corresponds to the trajectory
With the ocean currents represented as a 2D associated with the bicubic model and N corresponds to
vector field, a model was implemented in the Wolfram the trajectory associated with either the bilinear or
language using the fourth order Runge-Kutta algorithm biquadratic model. This approach allows us to quickly
and the three interpolation methods described above. compare the bilinear and biquadratic trajectories.
Inputs to the model are the starting point, step size, and This process was repeated for four 5 x 5 grids of
number of iterations. Each run of the model created three GPS coordinates in the North Indian Ocean the Arabian
trajectories — one each for the linear, quadratic and cubic Sea, the Bay of Bengal, northeast of Seychelles and
interpolation methods. For all runs, a step size of 0.5 was southwest of Sri Lanka as shown in Figure 7, and the total
used with 800 iterations. errors across the trajectories were compared.
Because there is no known ground truth with In addition, 16 GPS coordinates were evaluated
which to evaluate the models, the bicubic interpolation qualitatively by observing the trajectories. Eight of these
was used as a reference by which to evaluate the linear were chosen along the coast of India to be near river
and quadratic models as it represents the most outlets into the sea, and eight were spread across the Bay
mathematically accurate solution of the three. To of Bengal and the Arabian Sea. The processing times for
calculate the error between the linear and quadratic each of these runs was also recorded.
trajectories with the cubic trajectory, the Euclidean

32 Indian J.Sci.Res. 14 (1): 29-36, 2023

 CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN

Figure 7: The four regions of 5 x 5 grids of Figure 8: Bilinear and biquadratic error for each
coordinates on which the runs were tested trajectory along with a best fit line
RESULTS AND DISCUSSION Table 2 shows the time taken to calculate several
Figure 8 illustrates a graph of the total error trajectories with each of the interpolation methods. The
associated with both the bilinear and biquadratic models results are consistent across the runs, regardless of the
when compared to bicubic for all of the points tested. For starting points of the trajectories. Bilinear and biquadratic
a given trajectory, the bilinear error is represented on the interpolation are similar in the absolute times taken with
x-axis, and the biquadratic error is represented on the y- bilinear interpolation being slightly faster than
axis. A determination of the line of best fit yields the line biquadratic. However, bicubic interpolation results in
and is also overlayed on the plot. Because approximately a 10x increase in computation time over
the slope of the line is nearly one, points under the line both bilinear and biquadratic.
generally indicate samples where the biquadratic model Table 1: Comparison of bilinear and biquadratic
performed with less error than the bilinear model, and error across trajectories by region
points above the line indicate the contrary.
Bilinear error > Biquadratic error >
As can be seen from the figure, the points appear Region
Biquadratic error Bilinear error
to cluster around the line , suggesting the I 9 16
biquadratic model performs similar in accuracy to the II 13 12
bilinear model, if not slightly better. Table 1 also shows III 12 13
the count of trajectories for each region where the bilinear IV 9 16
error is greater than the biquadratic error and vice versa.
The biquadratic model performs better in 57% of the
trajectories.
Table 2: Time taken in seconds to generate the trajectories for the bilinear, biquadratic, and bicubic methods
Latitude Longitude Bilinear Biquadratic Bicubic
19 72.4 0.578843 0.612795 5.88402
13 80.5 0.566436 0.593552 5.88056
15.4 73.7 0.573438 0.607441 5.83616
12.9 74.6 0.570757 0.613081 5.87552
9.95 76.2 0.574994 0.600569 5.8592
8.64 78.15 0.58338 0.602804 5.93379
11.36 79.85 0.566767 0.60121 5.84195
15.66 80.98 0.572612 0.599231 5.88178
20.08 86.77 0.563314 0.609077 5.92286
21.4 91.02 0.562457 0.593427 5.86548
13.9 84.1 0.580862 0.609787 5.90423
13.9 88.1 0.580396 0.620646 5.99402
13.9 92.1 0.582617 0.619004 5.95934
19 70 0.583237 0.629257 5.87745
14 70 0.581369 0.612338 5.89972
9 70 0.581801 0.615364 5.94343

Indian J.Sci.Res. 14 (1): 29-36, 2023 33

34
Bilinear Model Biquadratic Model Bicubic Model

CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN
Starting Point: (12.5, 85)

Starting Point: (14, 63)

Indian J.Sci.Res. 14 (1): 29-36, 2023

Figure 9: Trajectories associated with the 3 interpolation methods with a GPS starting point (13.9, 88.1)
Indian J.Sci.Res. 14 (1): 29-36, 2023

Bilinear Model Biquadratic Model Bicubic Model

Starting Point: (-7, 56)

CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN
Starting Point: (-7, 85)

Figure 9: Trajectories associated with the 3 interpolation methods with a GPS starting point (13.9, 88.1)
35
 CHAKRAVARTHY: USING NUMERICAL METHODS TO TRACK OCEAN PLASTICS IN THE NORTH INDIAN OCEAN

These results indicate there is no significant computation times for each model were recorded, and
benefit in a model with biquadratic interpolation over with the bicubic model designated as the reference
bilinear interpolation for this application under the solution, the error was calculated for both the bilinear and
current conditions. biquadratic models. In addition, trajectories created by
each of the models were assessed qualitatively.
In addition to the quantitative evaluation above,
the trajectories were evaluated visually. All trajectories Across a set of 100 trajectories taken from
are not shown here for conciseness, but Figure 9 shows starting points in the Northern Indian Ocean, the error
the trajectories associated with the first GPS coordinate associated with the bilinear model appears to be similar to
for each region. In most cases, the trajectories starting that of the biquadratic model, indicating there is no clear
near the shore result in endpoints along a nearby superiority of the biquadratic model as might have been
coastline. This points to the hypothesis that most plastic expected for this scenario. The time taken to execute the
entering the ocean from points along the Indian coastline biquadratic model is marginally longer than that of the
result in beaching along the Indian coastline or in fewer bilinear model, so there also there is no clear advantage
cases along other countries in the North Indian Ocean (eg. of one model over the other with regards to computation
Sri Lanka, Maldives) (van der Mheen et al., 2020). time.
Trajectories originating further from land were There are numerous avenues for further study
taken to be at least 90 km from the nearest point on the including variation of the step size and more granular
Indian coastline. This ensured that no points in the first ocean current data, as both of these could show more
iteration of the cubic interpolation calculation would be 0. differentiation in performance between bilinear and
While this set produced more variation in the trajectories biquadratic models.
than the set of trajectories starting near land, all
REFERENCES
trajectories starting in the Bay of Bengal ended in the Bay
of Bengal, while all trajectories starting in the Arabian Egger M., 2023. ―The Other Source: Where Does Plastic
Sea ended in the Arabian Sea. This again agrees with the in the Great Pacific Garbage Patch Come
hypothesis that most debris originating in the Northern From?‖ The Ocean Cleanup, 17 July 2023,
Indian Ocean remains in the Northern Indian Ocean. www.theoceancleanup.com/updates/the-other-
source-where-does-plastic-in-the-great-pacific-
CONCLUSION AND AVENUES FOR
garbage-patch-come-from/.
FURTHER STUDY
Ritchie H., 2021. ―Where Does the Plastic in Our Oceans
The aim of this paper is to explore the use of Come From?‖ Our World in Data, 1 May 2021,
numerical methods to understand the movement of www.ourworldindata.org/ocean-plastics.
plastics in the ocean. With the ocean modeled as a vector
field, the fourth order Runge-Kutta algorithm was used to van der Mheen M., van Sebille E. and Pattiaratchi C.,
track the movement of a particle across that vector field, 2020. Beaching patterns of plastic debris along
with bilinear, biquadratic, and bicubic interpolation used the Indian Ocean rim. Ocean Sci., 16: 1317–
to interpolate points between discrete samples. The 1336, https://doi.org/10.5194/os-16-1317-2020.

36 Indian J.Sci.Res. 14 (1): 29-36, 2023