University of Ghana  http://ugspace.ug.edu.gh
Tracking pollutants using Lagrangian Coherent
Structures.
By
Amengor Cyril Makafui
(10638272)
June 2018
THIS THESIS IS SUBMITTED TO THE UNIVERSITY OF GHANA,
LEGON IN PARTIAL FULFILMENT OF THE REQUIREMENTS
FOR THE AWARD OF MASTER OF PHILOSOPHY
MATHEMATICS DEGREE.
University of Ghana  http://ugspace.ug.edu.gh
DECLARATION
This thesis was written in the Department of Mathematics, University of Ghana,
Legon from August 2017 to July 2018 in partial fullment of the requirements for the
award of Master of Philosophy degree in Mathematics under the supervision of Dr.
Joseph K. Ansong and Dr Anton Asare-Tuah of the University of Ghana.
I hereby declare that except where due acknowledgement is made, this work has never
been presented wholly or in part for the award of a degree at the University of Ghana
or any other University.
Signature:..................................................
Student: Amengor Cyril Makafui
cyril@aims.edu.gh
Signature: ...................................................
Dr. Joseph K. Ansong
Signature: ...................................................
Dr. Anton Asare-Tuah
i
University of Ghana  http://ugspace.ug.edu.gh
ABSTRACT
In steady ows, the notion of boundaries separating dynamically distinct regions is
not ambiguous. This is because the invariant manifolds of time-independent ows and
the critical points of time-periodic ows provide adequate information to determine
the behaviour of the solutions of these systems. However, for time dependent systems,
it is strenuous to determine the nature of their solutions due to their dependence on
time. Nevertheless, it was observed that just like steady ows, most time-dependent
systems have boundaries that prevent cross-mixing of dynamically distinct regions.
They are known as Lagrangian Coherent Structures(LCSs) and they are embedded
in time-dependent ows as robust structures that determine the ow pattern of uid
particles.
This project investigates LCSs and also employs a numerical method to compute the
Finite Time Lyapunov Exponent to detect these structures. Initially, the coherent
structures are dened as hyperbolic material lines that separate dynamically distinct
regions in an unsteady ow. Then, the LCSs are classied into attracting and re-
pelling structures based on their inuence on the time-dependent ow. Subsequently,
the LCSs are also dened as a second derivative ridges of the FTLE elds. This de-
nition is perceptible from the numerical computations of the double-gyre model where
the coherent structures are extracted as ridges of the computed FTLE elds. Further-
more, we employ the Finite time Lyapunov Exponent model to carry out numerical
simulations on satellite observed surface velocities along the coast of Ghana. The aim
of this realistic application is to determine the Lagrangian Coherent Structures that
are formed in geophysical ows. Finally, based on these results, we hypothesize the
implications of a crude oil spill along the coast of Ghana. It was realized that in the
event of a spill, the oil is likely to be conned to the coast temporarily due to the
concentration of repelling LCSs. Also, for a longer time interval the oil spill is likely
to be advected from the coastline.
ii
University of Ghana  http://ugspace.ug.edu.gh
DEDICATION
I dedicate this project work to my family and to the late Prof. Francis K. A Allotey
(AIMS Ghana, Founder).
iii
University of Ghana  http://ugspace.ug.edu.gh
ACKNOWLEDGEMENTS
My greatest thanks is to God Almighty for granting me good health and replenishing
my strength daily through the course of this project. Secondly, in line with the
famous quote of Sir Isaac Newton which says "If I have seen further it is by standing
on the shoulders of Giants", I would express my gratitude to Dr Joseph K. Ansong
for allowing me to stand on his shoulders. It has been a great learning experience.
Also to Dr. Anton Asare-Tuah, thanks for all the encouraging words.
To Dr. Ralph Twum; Head of Mathematics Department and Lecturers thanks so
much for giving me the platform to aspire for academic excellence. Also my sincere
gratitude to Mr. Bennett Foli for the help you oered to me in understanding and
working on the satellite data. To all my colleagues Some, Hannah and Esther I am
very grateful for all your ideas and help you rendered to me throughout the course of
this project. Finally not forgetting my God given friend Dzifa, thank you for all your
encouraging words and support despite your busy schedule.
iv
University of Ghana  http://ugspace.ug.edu.gh
Contents
Declaration i
Abstract ii
Dedication iii
Acknowledgements iv
1 Introduction 1
1.1 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Literature Review 5
3 Dynamical Systems 9
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Linear and Non-linear systems . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Phase Portrait of dynamical systems . . . . . . . . . . . . . . . . . . 11
3.3.1 Real and distinct eigenvalues . . . . . . . . . . . . . . . . . . . 12
3.3.2 Repeated eigenvalue . . . . . . . . . . . . . . . . . . . . . . . 12
3.3.3 Complex eigenvalues . . . . . . . . . . . . . . . . . . . . . . . 12
3.4 Time independent systems . . . . . . . . . . . . . . . . . . . . . . . . 13
3.5 Time-Periodic systems . . . . . . . . . . . . . . . . . . . . . . . . . . 19
v
University of Ghana  http://ugspace.ug.edu.gh
3.6 Time-dependent systems . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.7 Lagrangian Coherent Structures(LCSs) . . . . . . . . . . . . . . . . . 26
3.7.1 LCSs as material lines and surfaces . . . . . . . . . . . . . . . 27
3.7.3 Classifying LCSs as attracting or repelling material lines . . . 28
3.8 Finite Time Lyapunov Exponent(FTLE) . . . . . . . . . . . . . . . . 30
3.8.1 Maximum Stretch . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.8.2 Dening the LCSs as the second derivative of the FTLE . . . 33
4 Numerical Simulations and Analysis 34
4.1 Double gyre model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.1 Computational method for the Double gyre . . . . . . . . . . 36
4.1.2 Time-independent double gyre . . . . . . . . . . . . . . . . . . 36
4.1.3 Time-dependent double gyre . . . . . . . . . . . . . . . . . . 38
4.2 Numerical Computations of the FTLE on Satellite data . . . . . . . . 40
4.2.1 Velocity elds of the Atlantic Ocean . . . . . . . . . . . . . . . 41
4.2.3 FTLE plots for Numerical Computations on the Atlantic Ocean 43
4.3 LCSs in the Gulf of Guinea . . . . . . . . . . . . . . . . . . . . . . . 44
4.4 Discussion of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.5 Implications for Oil spills on the Gulf of Guinea . . . . . . . . . . . . 51
5 Conclusion 53
References 57
University of Ghana  http://ugspace.ug.edu.gh
Chapter 1
Introduction
Since time immemorial human beings have devised various techniques for solving
complex-life problems, however, in spite of all the technological advancements, one
major challenge remaining is the problem of pollution. Pollution is dened as the re-
lease of unwanted substances into the natural environment that cause adverse changes.
Also [3] describes pollution as the addition of any substance or any form of energy
such as radioactivity into the environment at a rate faster than it can be dispersed,
diluted, decomposed, recycled or stored in some harmless form. Pollution exist in sev-
eral forms such as land, water and air just to name a few. In August 2006, toxic waste
produced by a Singaporean oil company called Tragura Beheer BV, was ooaded
to an Ivorian waste handling company. The company then discarded the waste at
the port of Abidjan. Subsequently, this resulted in the contamination of their water
bodies and also caused a lot of infections due to the presence of Hydrogen Sulphide
in the waste. The aftermath of this event resulted in 17 deaths and about 100, 000
inhabitant needed critical medical attention [28]. Later on in April 2010, airborne
volcanic ash caused pandemonium in European airspaces. This natural disaster led to
the grounding of several commercial airlines because the atmosphere was not suitable
for air plane operations. These airline companies incurred huge debts. Finally in that
same month, the world recorded the greatest oil spill in history. An oil drilling rig
in the Gulf of Mexico exploded leading to the enormous spillage of oil on the Pacic
Ocean. Consequently it had adverse eects on aquatic life. Numerous aquatic plants
and organisms died. According to [4], 78% of the oil spill still remain in the Gulf of
Mexico.
Due to these three globally signicant events, it was of paramount importance to
propose a solution to help control and detect pollutants. Scientists from diverse elds
1
University of Ghana  http://ugspace.ug.edu.gh
suggested dierent approaches to help track the movement of these pollutants in u-
ids. One such method is to run complex three-dimensional ocean models such as
the Modular Ocean Model(MOM), the Hybrid Coordinate Ocean Model(HYCOM)
or the Massachusetts Institute of Technology General Circulation Model(MITgcm).
However, these global and regional ocean models are quite expensive to run numeri-
cally. They also require relatively more simulation time before results can be obtained
to help alleviate the impact of a disaster. Due to this reason, one of the proposed
methods to help track pollutants is to employ dynamical system theory. This theory
is a eld in mathematics used to describe complex systems, usually by employing
dierential equations [14]. A dynamical system is dened as a system that describes
the time dependence of a point in geometrical space [5]. There are three main groups,
namely: the time independent, the time- periodic and the time dependent systems.
For time-independent and periodic ows the behaviour of the system can be deter-
mined for any given point. The phase portraits of time-independent ows provides
additional details about the behaviour of their solutions. There are four main phase
portraits generated for these steady ows namely the node, focus, center and the
saddle point. The importance of studying these plots is to identify the concept of
boundaries in the ow which separate dynamically distinct regions. This is mostly
visible in saddle plots where the invariant manifolds act as the boundary within the
ow. Furthermore for time-dependent ows, due to their reliance on time, the be-
haviour of the systems is unpredictable. Nevertheless from numerous research work
done on time-dependent ows, it has emerged that there are structures in the ow
that act as boundaries thereby shaping the advection of particles. These barriers are
known as Lagrangian Coherent Structures (LCSs).
LCSs are time-evolving structures that shape trajectory patterns in complex dynam-
ical systems [21]. The Lagrangian term describes the movement of the particles
through uids. In 2D ows, LCSs appear as material lines (continuous, smooth)
curves that determine the direction of particles as they are advected in the ow. The
LCSs detection methods provide predictive capabilities, by identifying robust struc-
tures that are too strong to be destroyed by short-term uctuations in the ow eld.
These surfaces are classied into elliptic, parabolic and hyperbolic barriers depending
on the type of inuence they exhibit on neighbouring uid particles [8]. However this
project focuses on the hyperbolic barriers. These barriers behave as the invariant
manifolds of the steady ows by attracting or repelling other material elements in
the ow. Some examples of dynamical systems that employ the concept of LCSs are
Chlorophyll patterns in the ocean[12], surface drifters[16], clouds of volcanic ash[23],
2
University of Ghana  http://ugspace.ug.edu.gh
oil spills[17] and spores in the atmosphere[26]. One approach that may be used for
detecting these structures in unsteady ows is by running regional ocean models to
simulate the advection of tracers in the ow. The regional ocean models are based
on the Eulerian approach of viewing the ow eld. With this method, it focuses
on specic locations or regions in space as time elapses. As mentioned before it is
computationally expensive to run regional models, and in addition minor errors with
the initial conditions would result in a faulty result.
One of the most important diagnostic techniques for detecting LCSs in unsteady uid
ows is by computing the Finite time Lyapunov Exponent (FTLE). FTLE has proven
to be an eective tool in analyzing unsteady ows. Hitherto numerous computational
techniques that have been proposed to study LCSs have not been as ecient as
the FTLE. The FTLE is a scalar value which determines the amount of stretching
about the trajectory point of uid particles. It can also be dened as a scalar value
that determines the separation of minute particles in an unsteady ow over a nite
time interval [27]. Additionally, it varies as a function of time and space. Although
FTLE's are derived from Eulerian methods, it is considered as Lagrangian because it
considers individual uid particles in the ow. They are able to predict the succeeding
positions of uid particles based on the history of the ow. Therefore, the study of
LCSs and FTLE is not just for the sake of academic curiosity but they have yielded
ground-breaking results in tracking pollutants[22].
1.1 Research Objectives
The main objectives of this project include understanding the concept of dynamical
systems in relation to autonomous, time-periodic and time dependent systems. Our
specic objectives are:
1. To carry out a detailed study of LCSs and how they apply to uid dynamics.
2. To compute the FTLE eld in simple uid systems in order to detect LCSs in
steady and unsteady ows.
3. Perform numerical computations of FTLE from satellite observed surface veloc-
ity elds along the coast of Ghana.
Section 1.3: Outline
3
University of Ghana  http://ugspace.ug.edu.gh
Chapter 1 gives a brief introduction of some key terminologies used in the project.
Chapter 2 discusses the ndings of some research papers which would be vital for this
project. Chapter 3 explores the properties of the dynamical systems in detail, such as
the equations governing LCSs and the FTLE. In Chapter 4, we carry out numerical
simulations for the FTLE and the last chapter is the conclusion.
4
University of Ghana  http://ugspace.ug.edu.gh
Chapter 2
Literature Review
LCSs are dened by [7] as time-evolving material surfaces or lines that shape tra-
jectory patterns in complex dynamical systems such as geophysical ows and multi-
body problem such as the movement of celestial objects. Detecting these surfaces
in unsteady ows is challenging. Nevertheless, [7] proposed an algorithm that study
unsteady velocity data using the variational theory of hyperbolic LCSs in 2D ows.
Firstly, the paper shows that the LCSs are the most attracting or repelling material
lines in a geophysical ow. It considered the advection of a material surface denoted
as M(t) within a given time interval. In addition the normal repulsion rate was de-
termined for M(t) along the trajectory patterns. From the outcome, it was observed
that when the repulsion rate is greater than 1, then the material line is a repelling
LCSs and vice-versa. These coherent structures are collectively called the hyperbolic
LCSs. Subsequently, [7] carried out numerical simulations to determine strainlines in
ow. Strainlines are dened as curves tangent to the eigenvector eld of the Cauchy
green deformation tensor. The equation for the deformation tensor is given as
( )
Ct
∗
0+T (x ) = ∇F t0+T (x ) ∇F t0+Tt 0 t 0 t (x0).0 0 0
Where F t0+Tt (x0) represents the ow map from initial position x0 from initial time t0 0
to nal time t0 +T and ∇F t0+Tt (x0) denotes the Jacobian of the ow map. Addition-0
ally, the star represents the matrix transposition. The Cauchy deformation tensor is
symmetric and positive denite hence it admits two real positive eigenvalues and or-
thogonal real eigenvectors [7]. Also the eigenvectors are dened on a subspace called
the strain direction. One of the results of [7] is that strainlines are used to locate and
advect hyperbolic LCS over a given interval. However, some challenges were faced
5
University of Ghana  http://ugspace.ug.edu.gh
in the computations of the strainlines namely: sensitivity of strain direction and de-
generate points. For the deformation tensor, the eigenvectors are very sensitive to
computational errors. These errors are most likely near hyperbolic LCSs due to the
large deviations in close particles. This results in huge mistakes in the computation of
the Cauchy deformation tensor. Subsequently, degenerate points are produced when
the eigenvalues of the deformation tensor are equal. The tensor eld then becomes
a scalar multiple of the identity matrix therefore the eigenvectors at these points are
not well dened [7]. Finally, the algorithm was tested on 2D turbulence with the
velocity eld derived as the numerical result of the Navier-Stokes equation [7]. The
main ndings of this paper are that LCSs are attracting and repelling material lines
and the tangents of M(t) can be employed to determine the advection of LCSs using
a higher-order Hermite interpolation.
[11] dened LCSs as regions of qualitatively dierent tracer dynamics. Tracers are
particles embedded in geophysical ows that determine the ow patterns in the uids.
The two main types are active and passive tracers. For the former, they change the
ow pattern by transforming the properties of uids and the latter do not have any
eect on the ow. Therefore the latter is essential in determining coherent structures
in the ow. [11] investigated the role of the LCSs in the mixing of 2D turbulence.
Similar to [7] LCSs boundaries are dened as stable and unstable material lines, how-
ever, for turbulent ows due to the lack of regular time dependence the material lines
are tracked for all given time intervals to determine their stability or instability[11].
Since it's not possible, [11] proposed the nite-time stability and instability method
for detecting stable and unstable material lines. This approach is derived from the
concept of uniform normal hyperbolicity for dynamical systems. Therefore [11] de-
ned LCSs boundaries as material lines with locally the longest or shortest stability
or instability points in turbulent ows. Lastly, an algorithm was developed and tested
on a barotropic turbulence and was compared with the Lagrangian Statistics. It was
observed that the algorithm was able to detect mesoscale structures whilst the Sta-
tistical method was not. Moreover, the algorithm can also be used to isolate votex
cores and regions with distinct hyperbolicity levels [11].
[21] described FTLE as a measure of sensitivity of a uid particle's future behaviour
to its initial position in a ow eld. Initially, [21] investigated the two main per-
spectives of uid ow namely the Eulerian and Lagrangian approaches. The eulerian
method looks at the uid properties at a xed point in space and time. Conversely,
the Lagrangian considers the individual elements of the uid for varying space and
time. Therefore it is able to track changing velocity of individual uid elements as
6
University of Ghana  http://ugspace.ug.edu.gh
they are advected in the ow [21]. This makes the Lagrangian method the best tool
for studying material advection in geophysical ows. In understanding how the La-
grangian method is relevant in uid transport, [21] investigated the double-gyre. The
double-gyre is a model that is generated by a set of equations and it consist of two
oppositely swirling vortices whose strengths and locations vary with time [21]. In this
paper, [21] studied the advection of a circular drop of a dye in a 2D time periodic ow
in the interval [t0, t1]. It was observed at t1 that the dye is stretched and advected
throughout the periodic ow. Also, the dye is separated into two distinct regions
represented as red and blue curves by a boundary in the ow which is depicted as a
white line. This boundary represents the strongest repelling LCSs and it is responsi-
ble for the separation of the circular blob of dye. Therefore, the two distinct regions
can be considered as representing the crude oil and water. Due to the success of this
experiment, the LCS approach was employed to study the satellite data for the oil
spill on the deep-water horizon o the Gulf of Mexico. As a result, it was discovered
that the attracting LCSs determine the advection of the Oil spill by hindcasting.
Hindcasting is a process of predicting a previous event based on the data gathered
earlier. Subsequently, [21] realized that future events can be determined based on the
current information by forecasting. This method is able to locate major transport
barriers in unsteady ows such as eddies and vortices. However, the satellite data
must be comprehensive, reliable and up-to-date [21].
Finally, [24] introduced a dierent but a well founded approach to obtaining LCSs
in uid ows. The paper dened LCSs as ridges of FTLE elds. The ridges are
dened as special gradient lines of the eld that are transverse in the direction of
the minimum curvature [24]. First and foremost, two mathematical denitions were
proposed for a LCSs. The rst denition considered the LCSs as the curvature ridges
of the FTLE eld and it was derived from the principal curvature and the direction
of material surfaces in dierential geometry [18]. Also the LCSs were dened as the
second derivative ridge of an FTLE eld and it is relies on the Hessian of the FTLE.
The Hessian is given as [24]
∑ d2σTt (x)= 0 .
dx2
where σTt (x) represents the FTLE values at x. In addition, [24] determined the La-0
grangian properties of the FTLE. [24] realised that the FTLE values are Lagrangian
for large integration times, however, that does not necessarily mean the LCSs obtained
from these computations are Lagrangian. This is because the LCSs are obtained from
7
University of Ghana  http://ugspace.ug.edu.gh
the computations of higher derivatives of the FTLE elds [24]. Lastly, numerical com-
putations were carried out on three examples namely the double gyre model and the
Very high frequency(VHF) data along the coast of Florida. In the rst example, the
paper analyzed the time independent and dependent gyre models. In both models,
the low and high FTLE elds act as the attracting and repelling LCSs. Also for the
subsequent example the ux across the LCSs was determined to be less than 0.05%
[24]. Therefore, uid particles are not likely to cross the LCSs in the ow. In con-
clusion, LCSs are ridges in FTLE elds and also LCSs and FTLE provide a more
suitable approach in studying transport and mixing in unsteady ows.
8
University of Ghana  http://ugspace.ug.edu.gh
Chapter 3
Dynamical Systems
3.1 Introduction
This chapter introduces the basic understanding of the three types of dynamical
systems. Their properties would be illustrated by solving simple examples as a prelude
to the more complicated concepts of LCSs and FTLE to be discussed later.
Introduction to Dynamical Systems
Dynamical systems are systems in which a function describes the time dependence
of a point in a geometrical space. It can also be dened as a model that explains
the temporal evolution of a system. These systems are governed by a set of variables
known as the state variables. If the values for these variables are known for any given
time, then the state of the system can be determined. Therefore, the state variables
provide a detailed description of the system at any given point in time. The general
form of the dynamicalsystem is expressed as ′x (t; t0,x0) = v (x (t; t0,x0) , t) (3.1)x (t0; t0,x0) = x0
From equation (3.1), t, t0, x0 and x represent the initial time, nal time, starting and
nal positions respectively. In addition, x(t; t0, x0) and v (x (t; t0,x0) , t) represent
9
University of Ghana  http://ugspace.ug.edu.gh
the state variables and the velocity vectors of the dynamical system respectively.
Applications of dynamical systems are ubiquitous in our daily activities, some of
these include the swinging of a pendulum, the growth of bacteria and the movement of
uids. Based on their dependence on time, these systems are grouped into three main
categories, namely Time-independent, Time-periodic and Time-dependent systems.
3.2 Linear and Non-linear systems
The equations of all dynamical systems are classied under these two main categories.
They can either be linear or non-linear. Linear systems are systems that follow the
principle of superposition. The principle of superposition is dened by the following
laws:
i. Law of additivity: Let f(u1) = y1 and f(u2) = y2. Therefore, f(u1+u2) = y1+y2.
ii. Law of homogeneity: let f(u) = y and A be a scalar, then f(Au) = Ay
For a system to be considered linear, it must satisfy both laws. However, most real life
applications of dynamical systems such as the predator-prey model and the Naiver-
Stokes equation are dened by non-linear equations. Non-linear systems are systems
which do not follow the principle of superposition. With these systems, it is realised
that either only one of the laws of superposition is satised or none at all. How-
ever, there are two main methods employed in solving for non-linear systems namely:
the linearisation method and Lyapunov function. The linearisation method involves
linearising a non-linear system. The rst step includes determining the equilibrium,
xed or critical points of the non-linear system. Secondly, the Jacobian matrix of each
xed point is computed. Finally the eigenvalues and eigenvectors are determined for
the equilibrium points. This helps to determine the stability of the system based on
the eigenvalues. However, the linearisation method does not work on all non-linear
systems. For example, if complex eigenvalues with no real parts are obtained after
linearisation, it would be dicult to determine the nature of the non-linear system.
Therefore, in such cases the Lyapunov function can be used to determine the stability
of the equilibrium point. Lyapunov functions are scalar functions that are used to
determine the stability of non-linear dierential equations [13]. These functions are
based on the following denitions
10
University of Ghana  http://ugspace.ug.edu.gh
Consider a non-linear system of the form:
x′ = f(x)
where the function f = (f1, ..., fn) and the partial derivatives are assumed to be
continuous in an open set U ∈ Rn consisting of the origin. Let V (x) be a vector
continuous function dened in U and V (0) = 0. Then, the following denitions hold
I. V (x) is positive denite in U if and only if V (x) > 0 for all x 6= 0, x ∈ U.
II. For V (x) > 0 in U for all x ∈ U, V (x) is said to be positive semi-denite.
III. For V (x) < 0 in U for all x ∈ U, V (x) is said to be negative denite.
IV. For V (x) 6 0 in U for all x ∈ U, V (x) is said to be negative semi-denite.
These denitions are essential in checking for the stability, asymptotic stability and
instability of the trivial solutions or the xed points of non-linear autonomous systems.
In subsequent sections we will consider some examples of both linear and non-linear
dynamical systems.
3.3 Phase Portrait of dynamical systems
A phase portrait is dened as a geometrical representation of the solutions of a dynam-
ical system in the phase plane. Where the phase plane simply refers to the Cartesian
plane. Also, the phase portrait gives a detailed insight about the general properties
of the system by visualizing how the solutions of a system of dierential equations
would behave with time. The phase plots for linear systems can be determined by
the eigenvector method. With this method, the eigenvalues determine the direction
of the trajectories that separate the portrait into distinct regions based on the be-
haviour of the system. These trajectories are known as the separatrices. Also, the
critical/equilibrium points of a dynamical system are the points where the system has
a constant solution. In sketching the plots of linear systems, the sign of the eigenvalue
is critical. Positive eigenvalues generate plots with trajectories that move away from
the origin with time. They are known as unstable manifolds. Subsequently negative
eigenvalues produce plots with trajectories moving towards the origin with time and
they are referred to as stable manifolds.
11
University of Ghana  http://ugspace.ug.edu.gh
In the case of non-linear systems, we may have to apply the linearisation theorem and
then consider the system as a linear function. Also the Lyapunov function is applied
to non-linear systems that cannot be linearised. All computations for most parts of
the project are in R2.
3.3.1 Real and distinct eigenvalues
(i) When the eigenvalues are both positive, that is (0 ≤ λ1 ≤ λ2), the phase portrait
is an unstable node.
(ii) When the eigenvalues are both negative, that is (λ1 ≤ λ2 ≤ 0), the plot is a
stable node.
(iii) When the eigenvalues are positive and negative λ1 ≤ 0 ≤ λ2, the plot is a
saddle.
3.3.2 Repeated eigenvalue
(i) When the eigenvalues are positive (λ1, λ2 ≥ 0) then it is unstable.
(ii) When the eigenvalues are negative (λ1, λ2 ≤ 0) then it is stable.
3.3.3 Complex eigenvalues
(i) When the real part of eigenvalues are both positive eg: 2± 5i. The phase plot
is an unstable spiral.
(ii) When the real part of eigenvalues are both negative eg: −3± 4i. The plot is a
stable spiral.
(iii) When the real part is zero eg: 0± 3i. The eigenvalues are purely imaginary.
The next section explores the three main types of dynamical systems and looks at
their respective phase plots.
12
University of Ghana  http://ugspace.ug.edu.gh
3.4 Time independent systems
They are systems of ordinary dierential equations which do not depend on the inde-
pendent variable which is time. They are also referred to as autonomous systems. In
this section, we would consider a system in Rn before looking at cases in R2. These
collection of equations maybe written in the form

′
x = V1 (x1, x2, · · · xn)

1
′
x2 = V2 (x1, x2, · · · xn)
 ′x3 = V3 (x1, x2, · · · xn) (3.2) .. .′xn = Vn (x1, x2, · · · xn) .
The autonomous equations (3.2) can also be represented in the formx1  
 
V1 
x2  V2 x = and v =  ..   .. . . 
xn Vn.
Therefore generalizing equation (3.2) gives
′
x = V (x). (3.3)
The next section introduces some examples of autonomous systems and their general
solutions.
Example 3.4.1. Consider a dierential equation of the form: ax′′+bx′+c(x) = f(t).
The equation is a linear non-homogeneous second order ODE. The equation is said
to be homogeneous if f(t) = 0. The equation above can be represented as a linear
system of equations:
x′ = y (3.4)y′ = − cx− by.
a a
13
University of Ghana  http://ugspace.ug.edu.gh
Where equation 3.4 can also(be)expre(ssed as a)m(at)rix of the form:
x′ 0 1 x
′ = − − (3.5)y c b y
a a
For example, we want to consider the general solution of the equations: x′′+2x′−3x =
0. The equation can be written as a linear system of equations just as in (3.4) where
a = 1, b = 2 and c = −3: x′ = y (3.6)y′ = 3x− 2y
Using the matrix method, equation 3.6 is given as:
( ) ( )( )
x′ 0 1 x
′ = (3.7)y 3 −2 y
The eigenvalues of the matrix w∣ould be determ∣ ined by the determinant of the matrix.∣∣∣∣−λ 1 ∣∣∣∣ = 0 (3.8)3 −2− λ
Therefore we have the characteristic equation given as
λ2 + 2λ− 3 = 0 (3.9)
The eigenvalues obtained by solving the equation 3.9 are given as λ1 = −3 and λ2 = 1.
Since the eigenvalues are real and distinct, the eigenvectors are linearly independent.
Therefore, we determine the eigenvectors associated to each eigenvalue by substituting
the eigenvalues into equation 3(.8. Fo)r(λ1 )= −3( )
3 1 x 0
= (3.10)
3 1 y 0
From equation 3.10 we have the ex(press)ion 3x+ y = 0 in both cases. Let the value of
1
x = 1, the eigenvector is given as . Using the same approach, the eigenvector
−3
14
University of Ghana  http://ugspace.ug.edu.gh
( )
1
for λ = 1 is . The general solution is of the form C1 exp(λ1t)x0 +C2 exp(λ2t)x0.
1
Substituting the values of the eigenvecto(rs an)d eigenvalues(res)ults in
1 1
X(t) = C1 exp(−3t) + C2 exp(t) (3.11)−3 1
Where C1 and C2 represent the arbitrary constants and they can be determined for
given initial conditions. The phase plot is given below.
Figure 3.1: A Saddle
From gure (3.1), the black lines denote the separatrices and the red curves represent
the trajectories of the solutions. The blue arrows on the separatrices indicate the
direction of the separatrix. The separatrix with the blue arrows converging to the
centre is the stable manifold and it is obtained from λ = −3. This is the point where
trajectories moving towards the origin experience stability.Also λ = 1 results in the
separatrix with the blue arrows diverging from the centre. It is referred to as the
unstable manifold and trajectories moving towards these points experience instability.
Futhrmore, the stable manifold(stability) and unstable manifold(instability) occur as
t → −∞ and t → ∞ respectively. From (3.1) it is observed that the manifolds
determine the behaviour of the solutions represented as red curves as they move
towards the centre.
Example 3.4.2. Consider another system of dierential equations written in matrix
15
University of Ghana  http://ugspace.ug.edu.gh
form as: ( )
′ −3 −2
x = x (3.12)
5 −1
Let A represent the matrix in 3.12 and determine the eigenvalues of the Matrix A.
The determinant is given as
∣∣∣∣∣−3− λ −2 ∣∣∣
∣
= 0 (3.13)
5 −1− λ∣
The characteristic equation of 3.13 is given as
λ2 + 4λ+ 13 = 0.
The complex eigenvalues are given as λ = −2 ± 3i. Since the complex roots are
conjugates of each other, the eigenvectors are the same for both eigenvalues. The
eigenvector associated with the eigenvalue is dened by the linear system Ax = λx
where A is a square matrix, λ is a complex eigenvalue of A and x represents a vector.
Therefore we select λ =(−2 + 3i. )( ) ( )
−1− 3i −2 x 0
= (3.14)
(5 1)− 3i y 0
2
The eigenvector is given as and
−1− 3i
the solution of the equation is given as ( )
′ 2x = exp [(−2 + 3i) t] (3.15)
−1− 3i
The Euler formula is applied to equation 3.15 to separate it into real and imaginary
parts. [ ( )]
x′
2
= exp(2t) (cos 3t+ i sin 3t)
−1− 3i
16
University of Ghana  http://ugspace.ug.edu.gh
[( )] [( )]
′ 2 cos(3t) 2 sin(3t)x = exp(2t) + i exp(2t) (3.16)
− cos(3t) + 3 sin(3t) − sin(3t)− 3 cos(3t)
Let V1(t) represent the real part[(given as )]
2 cos(3t)
V1(t) = (3.17)− cos(3t) + 3 sin(3t)
and V2(t) represent the imagina[ry(part given as )]
2 sin(3t)
V2(t) = (3.18)− sin(3t− 3 cos(3t))
The general solution is a linear combination of the real and imaginary parts and it is
given as
XG = C1 exp(−2t)V1(t) + C2 exp(−2t)V2(t).
The phase portrait is given below.
Figure 3.2: A stable spiral
In gure (3.2), the red continuous curves represents the spiral and the black arrows
on the red curves depicts the direction of the behaviour of solutions of the system.
Since the real parts of both eigenvalue are negative, the trajectories generated will
converge towards the origin as t→∞. Therefore, the phase plot is a stable spiral.
17
University of Ghana  http://ugspace.ug.edu.gh
Example 3.4.3. Considera nonlinear autonomous system given asf(x, y) = x (3− x− 2y) (3.19)g(x, y) = y (2− x− y)
We want to nd the equilibrium points and determine the nature solutions of the
system, and plot the possible phase portrait of the system.
Since it is a non-linear system, we compute for the xed points of the system. The
xed points are obtained by setting the right hand sides of equation 3.19 to zero.
This condition is used to determine all the xed points. The rst xed point when
both values of x and y are zero is (0, 0).When x = 0, the xed point is given as (0, 2).
Also when y = 0, the xed point is (3, 0). Finally when both x and y are non-zero
the xed point is (1, 1).
The next step involves determining the linear stability of the xed point by computing
the Jacobian of the non-∣linear s∣yste∣m. The Jacobian is given∣ as:∣∣∣∣ df df ∣dx dy ∣∣∣ ∣∣∣∣3− 2x− 2y −2x ∣∣J = = ∣∣ (3.20)dg dg −y 2− x− 2y
dx dy
The eigenvalues of the matrix in 3.20 is determined by substituting the values of the
xed points into the equation and solving the characteristic equation for each point.
When (f (x, y), g(x, y)) = (0, 0) the ma(trix in)3.20 results in
3 0
J = (3.21)
0 2
The characteristic equation is given as λ2−5λ(+6)= 0. T(he)eigenvalues are λ = 3 and
0 1
λ = 2 and the corresponding eigenvector are and . Since both eigenvalues
1 0
are real positive and distinct, (0, 0) is an unstable node for the linearised system and
hence, by the linearisation theorem, also for non-linear systems.
For (f(x, y), g(x, y)) = (0, 2) the mat(rix 3.20 b)ecomes
−1 0
J = (3.22)
−2 −2
The characteristic equation is λ2 + 3λ + 2. The eigenvalues are λ = −1 and λ = −2
18
University of Ghana  http://ugspace.ug.edu.gh
( ) ( )
0.4472 0
and the eigenvectors are and . Both eigenvalues are real negative
−0.8944 1
and distinct, therefore (0, 2) is a stable node for the linearised system and likewise
for the non-linear.
When (f(x, y), g(x, y)) = (3, 0) the m(atrix 3.20)gives
−3 −6
J = (3.23)
0 −1
The characteristic equation is λ2 + 4λ + 3 = 0. T(he)eigenva(lues are λ)= −3 and
1 −0.9487
λ = −1 and their corresponding eigenvectors are and The xed
0 0.3162
point (3, 0) is also a stable node.
Finally for (1, 1), 3.20 results in ( )
−1 −2
J = (3.24)
−1 −1
The characteristic equation is given as λ2 + 2λ− 1 = 0.(The comp)lex eig(envalues a)re
√ 0.8165 0.8165
given as λ = −1± 2i and the eigenvectors are given as and .
−0.5774 −0.5774
The point (1, 1) is a stable spiral for the linearised system and hence, by linearisation
theorem, also for the non-linear systems.
Figure 3.3 represents the phase portrait of equation 3.19. From the gure, the points
(0, 2) and (3, 0) are the stable nodes and they are represented by the shaded circles.
Also, the points (0, 0) and (1, 1) are denoted by circular rings and they are the unstable
nodes. It is observed that the trajectories of the unstable nodes converge towards the
stable points. Finally, the point (1, 1) acts like the saddle in gure (3.1).
3.5 Time-Periodic systems
Time periodic systems are dierential equations that depend explicitly on time and
have their solutions recurring after a period of time (T ). Subsequent solutions are
given at 2T , 3T and etc. The phase portrait of these systems displays vividly the
repetitive behaviour of their solutions over time. Furthermore, their phase portrait
are mostly closed orbits in phase plane and continuous curves. They are generally
19
University of Ghana  http://ugspace.ug.edu.gh
Figure 3.3: A Phase portrait of the non-linear time-independent dynamical system.
expressed as
′
x = v(t, x) (3.25)
Time-periodic systems are obtained from physical, chemical and biological systems
that exhibit self-sustained oscillations. Some examples are the predator-prey model,
the harmonic equation and the simple pendulum. The last example would be expa-
tiated on because it would be benecial in understanding the concept of boundaries
in time-dependent systems.
Example 3.5.1. Consider the second order equation given below
x′′ + x3 = 0.
The equation is expressed as a rst-order equation.x′ = y (3.26)y′ = −x3
Employing the method of separation of variables, equation (3.26) is solved as shown
20
University of Ghana  http://ugspace.ug.edu.gh
below.
dy −x3
=
dx y
ydy = −x3dx
Integrating both sides results in∫ ∫
ydy = − x3dx
4
y2 = −x + 2C
2
Hence the general solution is given as
√ ( ) 14 2
y = ± x2 C − (3.27)
4
Figure 3.4: An Ellipse
In gure (3.4), the blue elliptical curves represent the solutions of the system. The
black arrows on (3.4) denote the direction of the solutions of the periodic system.
The phase plot is an ellipse.
21
University of Ghana  http://ugspace.ug.edu.gh
Example 3.5.2. Consider the simple pendulum with the equation given as
x′′ + sinx = 0.
The second order dierential equation is expressed as a rst-order of the formx′ = y (3.28)y′ = − sinx
The general solution is determined using the method of separation.
dy − sinx
=
dx y
ydy = − sinxdx
Both sides of the equation are integrated and the general solution of equation (3.28)
is given as
√ 1
y = ± 2 (cosx+ c) 2 (3.29)
The arbitrary constant c represents the parameter of the phase paths and it determines
the nature of the phase plots. For values of c that lie between −1 ≤ c < 1, the phase
plot generates closed circular orbits. Also, when c = 1, the phase plot produces curves
which are xed at the equilibrium points. When c > 1 the phase plots are open orbits
with a whirl-like motion.
22
University of Ghana  http://ugspace.ug.edu.gh
Figure 3.5: Phase plot of a simple pendulum
Figure 3.5 displays the phase portrait of a simple pendulum where the green circles
represent closed paths, the separatrices are also represented as the blue curves and
the red curves depict the whirling motion. The blue curves, just as the separating
lines in the saddle plot, act as the unstable and stable manifold. They are given at the
points (0, 0) and (±π, 0) . Therefore in the separatrix(blue curves), the pendulum can
only swing back and forth. As a result, the phase portrait of the solutions are circles.
Also, outside the boundary, the pendulum moves in one direction. Correspondingly,
the phase plots of solution are in a whirl-like form. Altogether the simple pendulum
example provides adequate information in understanding hyperbolicity in dynamical
systems.
3.6 Time-dependent systems
They are dynamical systems that depend explicitly on time, however unlike the time
periodic systems, their solutions are not recurring. Most realistic dynamical systems
such as the logistic growth equation and the movement of uids are classied under
this category. They are also referred to as non-autonomous systems and they are
23
University of Ghana  http://ugspace.ug.edu.gh
given as
x′ = v (t, x) (3.30)
The logistic growth equation would be solved below and the phase plot would be
generated.
Example 3.6.1. Consider the dierentia(l equat)ion given as
dP − P= kP 1 (3.31)
dt k
The above equation is often used to model population growth in living organisms.
Using the metho∫d of separation of variables, the equation is given as(dP ) ∫= kdt where P 6= 0 and P 6= K
P 1− P
k
Where the initial part of the equation is expressed as a partial fraction as shown
below
1 1 1
∫ ∫ ( ) =P 1− P ∫ +P (k − P )k
dP dP
+ = kdt
P (k − P )
ln(P )−(ln(k − P)) = kt+ c
k − P
ln = −kt− c
P
k − P
= exp (−kt− c)
P
= A exp (−kt) where A = exp(−c)
The solution is given as
k
P (t) = (3.32)
1 + A exp(−kt)
From the general solution, we compute for the initial population which is given as P0.
k k
P (0) = = (3.33)
1 + A exp(−k(0)) 1 + A
24
University of Ghana  http://ugspace.ug.edu.gh
Substituting equation 3.33 into equation 3.32, the general solution for the logistic
growth equation in terms of P0 is given as
k
P (t) = ( ) (3.34)
k − P0
1 + exp(−kt)
P0
The equilibrium points of the logistic growth equation is obtained by equating (3.31)
to zero. These points are P = 0 and P = k.
P
k
t
Figure 3.6: A phase plot of the logistic growth equation
25
University of Ghana  http://ugspace.ug.edu.gh
From gure (3.6), the red sigmoid shaped curve represent the logistic growth curve.
It can be observed that there is no growth when the time is zero. The dashed-lines
represents the carrying capacity which serves as a boundary. The carrying capacity
is dened as the maximum population size of organisms that the environment can
sustain. It can be observed that the population grows gradually until it gets to the
carrying capacity (k) where it maintains the same value after a period of time till it
decreases. Therefore, the population cannot exceed the carrying capacity. This is an
example of a dependent system where the population growth or the solution can be
determined at any point in time.
However, to determine the behaviour of most time dependent systems is cumbersome
because the solution of the dierential equation varies with time. In some of these
systems, there are regions of dynamically distinct behaviour which are presumed to
be partitioned by separatrices. An example is the movement of uids. Nonetheless,
the distinct regions change with time likewise the separatrices. For uid, even when
the velocity eld is of a simple form, the behaviour of the corresponding trajectories
of the uid can be complicated and seemingly unpredictable, an event called chaotic
transport [24].
3.7 Lagrangian Coherent Structures(LCSs)
LCSs are dened as distinguished surfaces of trajectories in a dynamical system that
exert a major inuence on nearby trajectories over a time interval of interest[10].
Also, [25] dened LCSs as mobile separatrices that divide the ow into regions of
qualitatively dierent dynamics. In steady ows such as time independent and time-
periodic ows, the stable and unstable manifolds provide concrete information about
the ow based on the critical points and the separatrices. However, for unsteady ows
such as time dependent ows, the critical points and separatrices changes even within
a short time interval making it dicult to determine the ow of tracer particles.
Therefore, the LCSs acts as boundaries within the ow that prevents cross-mixing
of uid particles that exhibit dierent dynamical behaviours. Furthermore, they are
embedded in the ow as well dened curves and they shape the ow of tracer particles
by two major approaches namely attraction and repulsion. The LCSs can be referred
to as attracting structures when time (T ) is less than zero and they are important in
studying transport of tracer particles in time-dependent systems. Also when T < 0,
we have more attracting structures in the simulated ow. When T > 0, the LCSs
26
University of Ghana  http://ugspace.ug.edu.gh
is referred to as repelling LCSs and particles on opposite sides of these structures
would separate from each other in the ow [25]. They are very useful in studying
mixing in ows and therefore identifying these structures is key in understanding
how ows separate in time-dependent systems. The attracting and repelling LCSs
are collectively referred to as hyperbolic structures. Although coherent structures are
ubiquitous in time-dependent ows, there is no well-dened mathematical theory for
computing them in non-autonomous systems.
3.7.1 LCSs as material lines and surfaces
LCSs are hyperbolic material surfaces/lines in 2D ows and they can be visualised
as continuous and smooth curves that shape the movement of uid particles in a
simulated ow[11]. They are time-evolving material surfaces along which passive
tracers (e.g. dye,temperature, or any material which does not interact with the ow)
develop distinct patterns exceeding typical time scales in the ow [8]. [11] proposed a
mathematical denition for material lines and surfaces in the study of unsteady ows.
The denition is given below.
Denition 3.7.2. Consider a 2D velocity eld given as
dx
= u (x, t) with nite time interval [tint, tfin] (3.35)
dt
where tint and tfin represent the backward integration time and forward integration
time respectively. The trajectories of the particles are represented as x(t) and they are
obtained from the 2D velocity eld equation. In uid mechanics, particle trajectories
are referred to as path-lines and they can be observed as smooth curves on the x-axis.
Therefore, a curve with initial conditions given as τ0 when advected by the velocity
eld for a period of time t which represents all values of time, the conditions of the
curve transforms into τt. The current conditions τt of the curve represent the material
lines. The advection of particle trajectories from an initial position (x0) at (tint) to
a later position in the ow at time (tfin) is called a ow map. The ow map is
mathematically given as
t
Θ fint : x0 → X (tfin; tint, x0) where x0 represent the initial positionint
When a deforming material line stretches over a 2D surface M in the (x,t) space, a
material surface is formed [11].
27
University of Ghana  http://ugspace.ug.edu.gh
3.7.3 Classifying LCSs as attracting or repelling material lines
This denition is based on LCSs being classied as material lines and it was proposed
by [7].
Consider a material lineM(t) with time t, when it is advected by the ow map into a
time-evolving material line M(t), Then the ow map for M(T ) = ΘTt M(t) with time
interval [t, T ].
Let x represent a point on the material line such that x ∈ M(t). Also let η denote
the unit normal at x on M(t) with the linearised ow map given as ∇ΘTt (x). To
be able to determine the behaviour of the material line M(t) as it is advected in
the ow, we determine the normal rate of repulsion. The normal rate of repulsion is
dened as the projection of the linearised ow map (∇ΘTt (x)) onto the normal (η) [15].
Therefore, we will denote it as βTt (x, η) and it is perpendicular to the linearised ow
map. Therefore, if βTt (x, η) < 1, it implies that M(T ) is overall attracting between
[t, T ] along the trajectory starting from x. Subsequently, if βTt (x, η) > 1, then M(t)
is overall repelling between [t, T ] along the trajectory starting from x [7].
The gure below provides an illustration of how LCSs are identied in unsteady ows.
Figure 3.7: The trajectories of particles in an unsteady ow
28
University of Ghana  http://ugspace.ug.edu.gh
From gure 3.7, the red and blue arrows denote particle trajectories of dynamically
distinct properties and the thin black curve separating the particle trajectories rep-
resents the LCSs. We see that the movement of the trajectories in the ow is chaotic
and could lead to mixing in uids of the same properties. Here, trajectories of simi-
lar properties ow only within their specied domain and they don't mix with other
particle trajectories of dierent properties. Furthermore, trajectories moving towards
the boundary tend to bend inwards thus moving away from the LCSs. Therefore, the
LCSs prevents cross-mixing of particle trajectories of dierent properties in the ow.
Thereby acting as the invariant manifolds and critical points of steady ows.
Figure 3.8: The trajectories of two particles with dierent dynamics in the ow
Figure 3.8, represents the advection of two particles at points A and B with initial
time t in an unsteady ow. After some time interval, it is observed that the two points
diverge from each other when it is integrated forward in time to new positions denoted
′ ′
as A and B . A similar phenomena occurs when it is integrated backwards in time
′ ′
from the nal position. In addition, the line A B denotes the distance between the
two points as they diverge from each other. The distance, separation or magnitude
between the two points can be computed using the idea of a nite time Lyapunov
exponent. The next section covers this in details.
Remark 3.7.4. We emphasize that these gures are only a schematic illustration to
help explain the concept of LCSs and FTLE in uid ows.
29
University of Ghana  http://ugspace.ug.edu.gh
3.8 Finite Time Lyapunov Exponent(FTLE)
The nite time Lyapunov exponent is a scalar value which determines the amount
of stretching/separation of the trajectories of particles when advected from an initial
position X ∈ D, where D is the domain, to a subsequent position Y ∈ D within a
nite time interval [t, t+ T ]. The FTLE is a scalar value that varies as a function
of space and time. For any given position in the ow, the maximum stretch at a
nite time can be determined using the FTLE. The next section considers how the
maximal stretch or the maximum FTLE is derived. In this section, vector values are
represented in bold faces.
3.8.1 Maximum Stretch
Consider an arbitrary point X ∈ D at initial time [t0] and another point Y ∈ D very
close to X. The position of y is given as Y = X + δX(t0) where δX(t0) represents
a perturbation of y from x. When these points are advected by the ow, the ow
map of X and Y are represented as ΘTt (X) and Θ
T
t (Y) respectively. The two points
diverge from each other with time. Therefore, to determine the maximum stretch of
both points at nal time (T ), the ow map is given as
δX(T ) = ΘTt (Y)−ΘTt (X)
where we can substitute the expression Y = X + δX(t) resulting in
δX(T ) = ΘTt (X + δX(t0))−ΘTt (X) . (3.36)
The right hand side of the equation above can be simplied by performing a Taylor
series expansion as shown below.
ΘT
dΘT ( ) d2ΘTX
(X + δX(t )) = ΘT (X) + t δX(t ) + t
(X) 2
t 0 t dX 0 d 2
δX(t0) + . . .
X
Substituting this expression into equation 3.36 and applying the Landau notation to
stand for the subsequent values of the expansion results in
dΘTt (X)δX(T ) = δX(t 20) +O (‖ δX(t0) ‖ ) (3.37)dX
30
University of Ghana  http://ugspace.ug.edu.gh
However, the subsequent values of the expansion which are represented by the Landau
notation are very small and therefore they are negligible. Therefore, equation 3.36
becomes
dΘT (X)
δX(T ) = t δX(t ) (3.38)
dX 0
Now because we are interested in determining the average stretching of the points,
we compute for the magnitude. √
T
‖ dΘδX(T ) ‖ = t (X)δX(t ) (3.39)
dX 0
Applying the idea of inner produc√ts, the above expression can be expressed as〈
‖ ‖ dΘ
T
t (x) dΘ
T (x) 〉
δX(T ) = δX(t0), t δX(t )
√ dX dX
0
〈 dΘT (X)∗ dΘT (X) 〉
= δX(t ), t t0 δX(t )dX dX 0
dΘT (X)
From [24], t is a matrix which represents the derivatives of the ow map with
dX
respect to their variations from the initial position and δX(t0) is a vector. The matrix
is given as  ∂ΘT T t (X1) ∂Θt (X1) ∂ΘT (X1) . . . t ∂X1 ∂X2 ∂Xn∂ΘT ( ) ∂ΘT TT  t X2 t (X2) ∂Θ (X )dΘ  . . .
t 2 
t = ∂X1 ∂X2 ∂Xn 
dX .. . . .. 
∂ΘT ( ) ∂ΘT ( ) ∂ΘTt Xn t Xn t (X ). . . n
∂X1 ∂X2 ∂Xn
The deformation gradient is computed for the matrix and it is denoted as
T dΘT
?
t (X) dΘ
T (X)
∆Θt (X) =
t (3.40)
dX dX
where ? represents the transpose and also the deformation gradient will give a sym-
metric matrix[24]. Now substitute e√quation 3.40 into the norm equation gives:
‖ δX(T ) ‖ = 〈δX(t0),∆ΘTt (X)δX(t0)〉 (3.41)
31
University of Ghana  http://ugspace.ug.edu.gh
Since we want to determine the maximum stretch in nite time, the equation above
would be rewritten as √
maxδX(t0) ‖ δX(T ) ‖ = 〈δX(t0),∆maxΘTt (X)δX(t0)〉 (3.42)
To obtain the maximum stretching, we use ideas from linear algebra and choose
the initial position δX(t0) to be aligned with the eigenvector corresponding to the
maximum eigenvalue of ∆ TmaxΘt (X) such that:
∆ TmaxΘt (X)δX(t0) = λmax∆δX(t0)
Where λmax∆ is the maximum eigenvalue that corresponds with ∆maxΘ
T
t (X) and it is
also associated with the vector δ̄X(t0). Substituting the eigenvalue into the equation
3.42 results in: √
maxδX(t0) ‖ δX(T ) ‖ = 〈δ̄X(t0), λmax∆δ̄X(t0)〉 (3.43)
Employing the concept of inner products where
〈X, AX〉 = A〈X,X〉 = A ‖ X ‖2
Equation 3.42 can be given as √
maxδX(t0) ‖ δX(T ) ‖ = λmax∆ ‖ δ̄X(t 20) ‖
√
maxδX(t0) ‖ δX(T ) ‖ = λmax∆ ‖ δ̄X(t0) ‖ (3.44)
As dened in [24], the FTLE is given as
1 √
σTt (X) = log λ| | max
(∆) (3.45)
T
where T represent the nal integration time. The integration time can be negative for
backward integration and also positive for forward integration. This explains the role
played by the absolute va√lue sign in equation 3.45. Therefore, from equation 3.45,( )
λmax (∆) = exp σ
T
t (X) | T | (3.46)
32
University of Ghana  http://ugspace.ug.edu.gh
Substituting equation 3.46 into the equation 3.44, the maximum stretch between the
two points in nite time is given by
( )
MaxδX(t0) ‖ δX(T ) ‖ = exp σTt (X) | T | ‖ δ̄X(t0) ‖ (3.47)
3.8.2 Dening the LCSs as the second derivative of the FTLE
[24] proposed a denition for detecting LCSs in an FTLE eld. This denition was
based on the numerical computations of the FTLE elds which would be studied in
Chapter 4. [24] dened LCSs as the second derivative of an FTLE eld based on this
explanation.
Denition 3.8.3. A second-derivative ridge of σTt is an injective curve c : (a, b)
satisfying the following conditions for each s ∈ (a, b) :
SR1: The vectors c′(s) and ∇σTt (c(s)) are parallel.
SR2: Σ(n, n)= min‖u‖=1Σ(u, u) < 0 where n is a unit normal vector to the curve c(s)
and Σ is though of as a bilinear form evaluated at the point c(s).
The ridges represent smooth curves in the unsteady ow and also the Hessian which
d2σT (x)
is denoted as Σ = t 2 . Additionally, since σ
T
t varies with time we indicate thedx
time on each curve c(s) at which the σTt is computed. eg: if the FTLE is computed
at time t1, c(s) becomes ct1(s).
33
University of Ghana  http://ugspace.ug.edu.gh
Chapter 4
Numerical Simulations and Analysis
This chapter introduces a numerical approach of detecting LCSs with the FTLE
method. Firstly, we carry out numerical simulations on two types of the double gyre
motions namely: time-independent and time-dependent gyres. The numerical compu-
tations generate the velocity elds and the FTLE plots of both the time-independent
and time-dependent ows. The double gyre provides suitable grounds for studying
complex time dependent systems. The next section involves FTLE computations on
ocean surface velocity elds of the West African Atlantic coast obtained from [1].
The main goal is to identify the LCSs and how they vary in time using realistic ow
elds. Afterwards, a section of the data which covers the Gulf of Guinea is extracted
for further analysis. The algorithm that computes the FTLE's of the ocean current
velocity data is developed by [6] and it was employed on a double gyre example and a
simulated ow in the Argentine basin. Initially, the algorithm loads the ocean current
velocity data into the workspace. Then, it requests for the following parameters: the
longitude, the latitude, the zonal and meridional velocities. Furthermore, the code
requests values for the following: the total integration time, the time step, the zonal
and meridional resolutions of the tracer grid, the method of integration, the direction
of integration and the number of ow maps to generate. Finally, using [2], plots for
the ow maps are developed. These steps are indicated in the sections on the double
gyre example and the numerical computations of the FTLE.
34
University of Ghana  http://ugspace.ug.edu.gh
4.1 Double gyre model
A gyre refers to any whirling mass of uid (the point where the uid ow revolves
around an axis line) formed naturally in uids. However, it is mostly studied in
oceanography, where the wind velocity on the ocean surface causes a large system of
circulating ocean currents to be formed. These systems act as transport barriers by
separating uid particles in a geophysical ow. A double gyre is formed when two
counter-rotating gyres revolve around each other. Also, the double gyre model is a
numerical representation of the double-gyre pattern that occurs mostly in geophysical
ows. This model provides a practical way of understanding the FTLE and LCSs. It
may be dened by the stream function
Ψ(x, y, t) = A sin (πf(x, t)) sin(πy) (4.1)
f (x, t) = a(t)x2 + b(t)x
where a(t) =  sin(ωt) (4.2)b(t) = 1− 2 sin(ωt)
where A represents the magnitude of the velocity vectors, ω is the frequency of os-
cillation and  denotes the amplitude of the motion of the separation point. The
eects of these parameters on the double gyre will be explained in later sections. The
velocity elds of the trajectories of double gyre is obtained by dierentiating equation
4.1 with respect with x and y respectively.U = −∂Ψ = −πA sin (πf(x, t)) cos(πy)∂y
V = ∂Ψ = πA cos (πf(x)) sin(πy)
df (4.3)
∂x dx
df = 2 sin(ωt) [x− 1] + 1
dx
where U represent the zonal velocity(velocities along the coordinates of the x-axis)
and V denotes the meridional velocity (velocities along the coordinates of the y-axis).
The equations above will generate the trajectories of the gyre as well as depict the
LCSs in the ow.
35
University of Ghana  http://ugspace.ug.edu.gh
4.1.1 Computational method for the Double gyre
This section explains the steps taken in the numerical integration of the double gyre
model. Firstly, we carry out numerical computations to determine the zonal and
meridional velocities from equation 4.3. These velocities are computed over a cartesian
grid with the x-axis(longitude) given as [2, 0] and y-axis also given as [0, 1] . The
parameter values given in equation 4.3 have major implications on the double gyre
model. When the magnitude of the velocity vectors (A) is small the distance between
the counter-rotating gyres and the boundary is minimized and vice-versa when it is
large. Smaller values of A are preferred since it produces more rened curves close
to the separatrix. Also, ω determines the frequency of oscillations throughout the
integration time. More oscillations are completed for smaller frequency with larger
integration time and conversely for large frequency in smaller integration time. Lastly,
the value of  controls the movement of the separatrix. A matlab function called quiver
is used to plot the velocity vectors of the double-gyre. Also, to determine the FTLE
elds, the Runge-Kutta algorithm developed by [2] is used to advect uid particles
from an initial position(x0) at time(t) to a subsequent position x(t+ T ; t, x0) at time
t + T. The algorithm advects the uid particles by integrating the velocity dataset.
Afterwards it computes for the spatial gradient by determining the nite dierence
for each point in the ow map. Finally, the Runge-Kutta algorithm determines the
FTLE values (σ) at each point by substituting the gradient values into 3.45. This
process is repeated for subsequent ranges of time. A matlab function called pcolor is
used to plot the FTLE elds.
4.1.2 Time-independent double gyre
A time independent gyre is generated when  is set to zero in equation 4.3. When
 = 0, the separatrix remains xed at the centre for the total integration time. The
equations governing the ow are given by:
U = −πA sin (πf(x)) cos(πy)
V = πA cos (πf(x)) sin(πy) (4.4)df = 1
dx
For the numerical simulations, a total integration time of 20 and a time step of 0.1 is
selected. Also, the selected values of the parameters for the numerical computation
36
University of Ghana  http://ugspace.ug.edu.gh
of the double gyre are given as A = 0.1,  = 0 and ω = 2π/10. These values are akin
to the values used for the double gyre example in [24].
(a) (b)
Figure 4.1: The plot of the trajectories of the velocity vectors and the FTLE plot respec-
tively
37
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.1a, indicates the trajectories of the gyre moving in opposite directions sep-
arating the plot into two equal but dynamically distinct regions. Also, the gyres are
separated at the points [1, 0] and [1, 1] by heteroclinical trajectories. Heteroclinical
trajectories act as the separatrix by the dividing distinct regions in the ow. Since
 = 0, the hetoroclinical trajectories remain xed at the points [1, 0] and [1, 1] at all
times. Figure 4.1b represent the 2D plot of the double gyre indicating the FTLE
ridges in the ow. From the gure, the highest repelling LCSs are represented as red
curves that lie on the separatrix. The magnitude of the FTLE values are shown on
the colour-bar. There is a quick separation of uid particles at this point in the ow.
4.1.3 Time-dependent double gyre
Time-dependent gyre occur when the value of () is not zero in 4.3. For the computa-
tion of the time-dependent gyre, the parameter values used for the time independent
gyre are maintained except  which is given as 0.1. Also, the total integration time
and the domain is maintained for the numerical simulations of the time-dependent
double gyre. Furthermore, the intervals was increased from 0.1 to 5 to provide a clear
view of the velocity vectors plots. Choosing a large interval for the velocity elds does
not result in large errors. It produces the same result as that of a smaller interval.
The only dierence is that larger intervals give a more vivid representation of the
numerical simulation as compared with the smaller intervals. The velocity vectors
and the FTLE plots are shown below.
38
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.2: The plot of the trajectories of the double gyre at t = 1, t = 3, t = 6 and t = 19
respectively.
The gures show the velocity vectors of the double gyre at t = 1, t = 3, t = 6 and
t = 19. Figure 4.2A is similar to what we see for the time independent double gyre
case shown in gure 4.1A. The two gyres are of similar size with ows moving in
opposite directions. In gure 4.2B, the position of the heteroclinical trajectories shift
rightwards and they are xed at the point [1.2, 0] and [1.2, 1]. This causes the gyres
to have distinct sizes with the initial gyre larger than the other gyre. Also, for gure
4.2C, the heteroclinical trajectories move back to the initial position as indicated in
4.2A and the gyres are of the same size. The time-dependent gyres are said to have
completed a period at t = 6. This also occurs for t = 12 and t = 18. Therefore, the
time-dependent double gyre completes 3 periods during the total integration time.
Finally, in gure 4.2D the heteroclinical trajectories move leftwards to the position
[0.8, 0] and [0.8, 1]. Also, the gyres are of dierent sizes with the second gyre larger
than the initial.
39
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.3: The FTLE plot at t = 1, t = 3, t = 6 and t = 19 respectively.
Figures 4.3A, 4.3B, 4.3C and 4.3D represent the corresponding FTLE eld plots of
the vector velocities given above. For all gures, it can be observed that the red curves
represent the magnitude of the High FTLE values . They are given in the interval
[0, 0.2] on the colour-bar. At the position of the red curve, trajectories experience the
quickest divergence in the ow. Therefore, the red curves represent the repelling LCSs
as they divide the ow into two distinct regions. In addition, the magnitude of the
Low FTLE values lie in the range [0, 0.14] on the colour-bar. The rate of divergence
of velocity vectors is reduced as you move down the colour-bar. Since the FTLE plots
are in 2D, the attracting and repelling LCSs can be visualized as material curves.
4.2 Numerical Computations of the FTLE on Satel-
lite data
In this section, the importance of the LCSs acting as transport barriers would be
demonstrated on a satellite data of the Atlantic Ocean. The data is obtained from [1]
and it spans from 1st October 2014 to 13th June 2018. The data les are sorted from
the initial date to the last date of data collection. Afterwards using a Matlab toolbox
40
University of Ghana  http://ugspace.ug.edu.gh
called Nctoolbox which provides read only access to common data model datasets,
the parameters of the dataset are extracted. They consist of the depth, salinity, tem-
perature, meridional velocity, height, time(t), longitude, latitude, zonal velocity and
mixed layer depth. For our analysis, we only extract the velocities at the surfaces.
The values of the meridional velocity (U), zonal velocity (V), latitude (lat), longi-
tude (lon) and time(t) are extracted from the data using a Matlab code developed
by [20]. The longitude and latitudes of the entire dataset is given as [-35◦W,15◦E]
and [-9.8750◦S,30.1250◦N] respectively. The velocities are given as 3D matrices where
the rst, second and third dimension consists of the longitudes, latitudes and time.
Similar to the double-gyre models, the Runge-Kutta algorithm is used to advect the
uid particles by integrating the velocities. Some other values that are important in
the computational procedure are the total integration time, the time step of integra-
tion, the axis resolution and the direction of integration. The total integration time
determines how long the uid particles are advected in the ow. Longer time interval
produces more rened LCSs in the ow. The time step of integration divides the to-
tal integration time to determine the next position of the uid particles in the ocean
ow. Also a ner resolution produces well rened LCSs. Furthermore, the numerical
integration is either performed backward or forwards in time. Backward integration
helps in determining the position from which a particle was advected over a given
time step whiles the forward integration time computes for the position to which a
particle will be advected during a time step. Finally, the Runge-Kutta algorithm is
used to compute for a selected number of ow maps based on the values selected for
the parameters and using a Matlab package produced by [2]. For the purposes of
these project, we are mainly interested in observing the eect of the repelling LCSs
on a simulated ow therefore we employ the backward time integration.
4.2.1 Velocity elds of the Atlantic Ocean
Just like the velocity elds of the double gyre models, the velocity plots of the satellite
data reveal important details about the geophysical ow. It indicates the movement
of uid ows. Figures 4.4A and 4.4B represent the surface velocities in the Atlantic
Ocean.
41
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.4: The surface velocity elds in A the Atlantic Ocean and B zoomed view of the
velocity elds at Ghana's coastline respectively.
42
University of Ghana  http://ugspace.ug.edu.gh
From gure 4.4A, the trajectories are represented as blue arrows and the white section
of the gure represents the land mass of West African and parts of north Africa.
Furthermore, gure 4.4A also indicates the land mass of countries along the Atlantic
ocean such as Ghana, Cote d'ivoire, Liberia and a lot more. Figure 4.4B represents a
enlarged section of gure 4.4A, specically Ghana's coastline. However, our interest
lies in the movement of uid trajectories in the simulated ow which is clearly shown
in gure 4.4B. In both diagrams, it can be observed that the uid ows are exhibiting
chaotic movement.
Remark 4.2.2. From gure 4.4A, there is a white straight line formed at longitude
0◦; this is a defect in the satellite velocity data. This defect is also evident in all the
FTLE elds plots in subsequent sections.
4.2.3 FTLE plots for Numerical Computations on the Atlantic
Ocean
Figures 4.5a and 4.5b depict the evolution of the FTLE elds in the entire West
Africa Atlantic coast. However only the rst and last ow maps are plotted to give
a clear representation of the evolution of the LCSs in the ocean ow. The resolution
on the tracer grids for the satellite data is 0.25 and also the integration is performed
backwards in time. The total time of integration and the time step used for 4.4 is
maintained.
Figure 4.5: The FTLE elds of the Atlantic Ocean at T = 1 and T = 4 respectively.
43
University of Ghana  http://ugspace.ug.edu.gh
In the gures, the wiggly lines in the plot represent the presence of numerous LCSs,
and those that are whitish indicate the most repelling/attracting structures. From
4.5a there is a concentration of high repelling LCSs between longitude 30◦ and latitude
11◦ N. However, in 4.5b, there are few high repelling LCSs formed at the same point.
We next zoom into the Gulf of Guinea to investigate the nature of the LCSs around
the coast of Ghana.
4.3 LCSs in the Gulf of Guinea
The longitude and latitude lie in the range [-10◦W,10◦E] and [-6.25◦S,6.125◦N]. Six
FTLE plots are generated for each total integration time given as: T = 50, T = 70
and T = 100. Also, the six FTLE plots represent the snapshot of the ow maps for
the following time intervals: T = 1, T = 5, T = 9, T = 13, T = 17 and T = 21. This
helps in giving a vivid view of the progression of the LCSs in the Gulf of Guinea.
The resolution is 0.25 and it is maintained for all the integration times. In addition,
computations were performed for backward time integration.
44
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.6: The rst, fth, ninth, thirteenth, seventeenth and twenty-rst ow maps of
the Gulf of Guinea with total integration time of 50 days and time step of 5 days.
45
University of Ghana  http://ugspace.ug.edu.gh
From gure 4.6, the LCSs are visualized as time dependent curves in the uid ow.
With the repelling LCSs represented by yellowish and whitish curves and the attract-
ing LCSs are also represented by the reddish curves and black patches in the ow.
Also the grey section in the gures represent the coastlines of Cote D'Ivoire, Ghana,
Togo, Nigeria and parts of Cameroon. Figures 4.6A and 4.6B indicate attracting
LCSs in many parts of the ow with the highest attracting point for 4.6C located
between 5◦E,1◦S and 10◦E,1◦S . Also for 4.6B the highest attracting point are lo-
cated between 10◦W,1◦S and 0◦E,1◦S. Furthermore, repelling curves are formed in
most parts of the ow with a concentration of high repelling curves centred at about
0◦E,4◦N for both gures 4.6A and 4.6B. Subsequently, gures 4.6C and 4.6D also
indicates less attracting LCSs formed in the Gulf of Guinea and there are a lot of re-
pelling LCSs stretching along 0◦E,4◦N. Figure 4.6E shows the highest case of repelling
LCSs stretching from the Coast of Ghana to further parts of the Gulf of Guinea. It
is most likely that there will be a rapid divergence of uid particles at this point in
the simulated ow. Finally, in gure 4.6F, the number of repelling curves decreases
as compared with 4.6E. The dominant repelling curve stretching from 10◦E,4◦N to
10◦W,1◦S. Also the highest attracting point is located between 0◦E,1◦S and 5◦E,1◦S.
At these points in the ow, there will be a high rate of mixing of uid particles. There
is a thin straight line at 0◦E that divides the FTLE plots into two equal parts. This
is as a result of a defect in the satellite data.
46
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.7: The rst, fth, ninth, thirteenth, seventeenth and twenty-rst ow maps of
the Gulf of Guinea with total integration time of 70 days and time step of 2 days.
47
University of Ghana  http://ugspace.ug.edu.gh
In gure 4.7A there are more attracting LCSs in the simulated ow with less repelling
LCSs. Figure 4.7B shows an increase in the repelling curves with high FTLE elds
formed at 5◦E,3◦N. Also in gure 4.7C, there is a decline of repelling LCSs curves with
the major repelling line stretching from 0◦E,6◦S. Figure 4.7D indicates an increase in
the repelling curves as compared with 4.7C. Figure 4.7E shows a particular repelling
LCS that bounds 0◦E,4◦N and 10◦W,1◦S. Fluid particles within this region would be
conned for sometime due to the repulsion eect of the repelling LCSs. Finally, gure
4.7F shows only traces of repelling curves in the ow.
48
University of Ghana  http://ugspace.ug.edu.gh
Figure 4.8: The rst, fth, ninth, thirteenth, seventeenth and twenty-rst ow maps of
the Gulf of Guinea with total integration time of 100 and time step of 10 days.
49
University of Ghana  http://ugspace.ug.edu.gh
Figures 4.8A, 4.8B shows traces of repelling LCSs in the Gulf of Guinea. In gure
4.8C, there is one main curve that stretches across the Gulf of Guinea with traces
of high FTLE elds in other parts of the ow. Figures 4.8D and 4.8E show only
traces of repelling LCSs scattered within the ow. Figure 4.8F has the highest case
of repelling LCSs in the simulated ow as compared with all the other gures.
4.4 Discussion of Results
From all the numerical computations on the Gulf of Guinea, it was observed that
there was a lot of repelling LCSs in the ow particularly centred along the coastal
belt of Ghana at 50 days. The coastal belt of Ghana is our region of interest. When
the total integration time was increased to 70 days, there were still regions of high
FTLE along the coast of Ghana, however, the concentration of repelling curves was
less as compared with the initial integration time. Also when the total integration
time was increased from 70 to 100 days, there was a further decline in the number
of repelling curves in the ow. Therefore, it was realized that more repelling LCSs
are formed for shorter integration times and for longer integration times there are
less repelling LCSs in the simulated ow. This may be due to the fact that at longer
integration time, the repelling curves have exited the Gulf of Guinea. The main goal
of this numerical simulation was to determine the nature of the LCSs located along
Ghana's coastline so as to be able to predict what will happen in the event of an
oil spill. In the case of an oil spill along the coast of Ghana at T = 50, due to the
presence of high repelling curves close to the coastline, the oil spill will most likely
be pushed onshore. At T = 70, it is most likely that the oil spill will be conned to
the coast of Ghana because the repelling LCSs would serve as transport barriers by
preventing the oil spill from crossing over to a dierent part of the Gulf of Guinea. At
this point, the oil spill can be contained, isolated and cleaned to prevent further harm
to aquatic organisms and also humans. Finally, for T = 100, due to the decrease in
the repelling LCSs and also there being a lot of attracting LCSs, the oil spill would
likely be advected along the main repelling LCSs away from the coast of Ghana.
However, oil spills have a smoothing eect on water and therefore it can cause a
change in the dynamics of the simulated ow. Since this has not been looked at, it
can be considered in future works so as to provide a better understanding of the ow
dynamics in the event of a spill.
50
University of Ghana  http://ugspace.ug.edu.gh
4.5 Implications for Oil spills on the Gulf of Guinea
Crude oil spill is the release of liquid hydrocarbons into the environment from oshore
platforms, drilling rigs and wells. Due to it's immiscible property and calming eect
on water, it has the propensity to spread for several nautical miles forming an oil
slick. There has been an upsurge of crude oil spills since the 1990's and in most third
world countries there has not been any proper monitoring system in place to check for
these spills. This is because most of the existing models such as the remote sensors
are very expensive to use and also due to the complexity of some other models. Oil
spills have a hazardous eect on aquatic plants and animals. For the former it cuts
o the supply of sunlight and air over the area it covers and with the latter due to
it's toxic components it kills organism that are exposed to it.
In Ghana, majority of the populace depend on the marine waters for shing. Ghana
consumes approximately 950, 000 metric tons of sh annually and also import 60
percent of its sh. Most of these products come from the Gulf of Guinea. In light
of the fact that Ghana extracts crude oil on a daily basis from the Gulf, it is very
alarming that it also obtains other products from the same source. Crude oil has
been shown to contain poisonous chemical constituents which can aect most marine
organisms when it is ingested or by external exposure. These organisms can either die
as a result or they are harvested by sherfolk who sell them for domestic consumption.
Subsequently, consuming these contaminated marine organisms can lead to some short
term diseases such as: memory loss, headache, nausea and vomiting and a lot more.
Also, it can cause some long term diseases such as: Lung cancer, skin cancer and
leukaemia. Therefore, the implications of oil spills on the gulf of Guinea are to be
taken seriously.
Due to the fact that there have not been any reported spills or there have been
smaller spills which have not been termed signicant enough to report. It is therefore
important to implement a model which can be used to track crude oil if there should
be a spill. The FTLE model can be used to detect trajectories of the particles in a
unsteady ow. Using the Cauchy green strain tensor, it computes for the eigenvectors
of the system. Based on the backward and forward time integration, the FTLE eld
is able to expose hyperbolic LCSs in the ow by numerical computations using the
velocity elds. Both integration steps are vital in determining repelling structures in
the ow. For the backward time integration, it computes the position from which a
particle was advected over a time step [20]. Also, forward time integration computes
51
University of Ghana  http://ugspace.ug.edu.gh
the position to which a particle will be advected during a time step[20]. It is hoped
that computations of the FTLE eld could be used to forecast the movement of
substances in the ocean. The rst stage of forecasting is called now-casting and it is
dened as the prediction of the current state of a time dependent system using the
available information [21]. This helps to detect transport barriers such as vortices
and ocean eddies in the ow.
52
University of Ghana  http://ugspace.ug.edu.gh
Chapter 5
Conclusion
Transport and mixing are important concepts in tracking of uid particles in time de-
pendent ows as they have been shown to detect coherent structures that determine
the shape of the ow. This project, employed the Finite Time Lyapunov Expo-
nent (FTLE) method to help identify Lagrangian Coherent Structures (LCSs)in an
unsteady ow, with a long term plan of tracking oil spills along the coast of Ghana.
Chapter 1 introduced the problem of pollution,the motivation of this research and the
concepts of LCSs and FTLE. The second chapter reviewed some research papers on
LCSs and the FTLE. Chapter 3 starts with a comprehensive study of time-dependent
systems. The solutions and phase plots of these systems presented in detail. For
instance, the saddle plots of autonomous system and the periodic plots of the sim-
ple pendulum. It was realized that they provided a standard for understanding the
notion of stable and unstable manifolds which separate dynamically distinct regions.
Two denitions of the LCSs were presented. First approach considers it as material
line for which a denition was given. In the second approach the LCSs is dened as
a second derivative ridge of the FTLE eld.
Chapter 4 mainly comprises of the double gyre model and the detection of LCSs in the
Gulf of Guinea. The rst part looked at the time-independent and time dependent
double gyre example. For the former, it was realized that the boundary remained at
the same position after the total integration time. This was due to the amplitude of
motion being set to zero. Also, for time dependent gyre, the value of amplitude of
motion was varied and this causes the boundary to move. The velocity vector plots
displayed the trajectories of the gyres rotating in opposite directions. Also the velocity
vectors clearly showed the movement of the boundary within the time-dependent gyre.
The FTLE plots indicated the highest repelling ridges as red curves and it is likely that
53
University of Ghana  http://ugspace.ug.edu.gh
uid particles close to these points will experience a rapid divergence. Secondly in this
chapter, the FTLE model was applied on a real satellite data on the ocean velocity
elds of the coast of West Africa with our region of interest being the coastline of
Ghana. However, for a better analysis to be made, we considered the Gulf of Guinea
to enable us investigate the progression of the LCSs in the simulated ow. The surface
velocities, the longitude and latitude was extracted for the numerical simulations on
the Gulf of Guinea. From the analysis it was realized that more repelling LCSs are
formed along the coast of Ghana for a shorter time of integration as compared with
a longer time. This nding was important in investigating the occurrence of an oil
spill along Ghana's coastline. The main result suggests that it is most likely oil spills
will be conned at the coast of Ghana for shorter periods however as time progresses
the spill would most likely be advected away from the coastline of Ghana. One of the
limitations of this project is that the velocity elds that we obtained did not include
active tracer like oil particles.
An extension of this project would involve advecting a tracer (pollutant) in the sim-
ulated ow and computing for the velocity elds to determine the advection path of
the tracer as time progresses. This would serve as a good grounds to identify other
transport barriers such as eddies, vortices and fronts that may be embedded in the
ow along the coast of Ghana. As a result, this will then enable us to be able to give
an accurate prediction of the advection path of an oil spillage in case it occurs.
54
University of Ghana  http://ugspace.ug.edu.gh
References
[1] (Accessed in June 2018). PRODUCT USER MANUAL FOR Global Ocean Multi
Observation Products. COPERNICUS MARINE ENVIRONMENT MONITOR-
ING SERVICE (April,2018) CMEMS-MOB-PUM-015-001.
[2] (Accessed 19 July 2018). A Mapping package for Matlab. Mmap(May 20th, 2018)
http://www.eos.ubc.ca/~rich/map.html.
[3] Britannica, E. et al. (2009). Encyclopædia britannica. Chicago: Common Law.
[4] Campbell, C. J. and Laherrère, J. H. (1998). The end of cheap oil. Scientic
American, 278(3):7883.
[5] Devaney, R. L., Siegel, P. B., Mallinckrodt, A. J., and McKay, S. (1993). A
rst course in chaotic dynamical systems: theory and experiment. Computers in
Physics, 7(4):416417.
[6] D.L.Vokov (April, 2008). Manual on how to use OceanFTLE1.0 Matlab package.
[7] Farazmand, M. and Haller, G. (2012). Computing Lagrangian coherent structures
from their variational theory. Chaos: An Interdisciplinary Journal of Nonlinear
Science, 22(1):013128.
[8] Hadjighasem, A. (2016). Extraction of Dynamical Coherent Structures from Time-
Varying Data Sets. PhD thesis, ETH Zurich.
[9] Haller, G. (2011). A variational theory of hyperbolic Lagrangian coherent struc-
tures. Physica D: Nonlinear Phenomena, 240(7):574598.
[10] Haller, G. (2015). Lagrangian coherent structures. Annual Review of Fluid
Mechanics, 47:137162.
[11] Haller, G. and Yuan, G. (2000). Lagrangian coherent structures and mixing in
two-dimensional turbulence. Physica D: Nonlinear Phenomena, 147(3-4):352370.
55
University of Ghana  http://ugspace.ug.edu.gh
[12] Huhn, F., Kameke, A., Pérez-Muñuzuri, V., Olascoaga, M., and Beron-Vera, F.
(2012). The impact of advective transport by the South Indian Ocean Countercur-
rent on the Madagascar plankton bloom. Geophysical Research Letters, 39(6).
[13] La Salle, J., Lefschetz, S., and Alverson, R. (1962). Stability by Liapunov's direct
method with applications. Physics Today, 15:59.
[14] Michel, A., Wang, K., and Hu, B. (2001). Qualitative theory of dynamical sys-
tems, volume 239. CRC Press.
[15] Nave Jr, G. K. and Ross, S. D. (2017). Trajectory-free calculation of attracting
and repelling manifolds. arXiv preprint arXiv:1705.07949.
[16] Nencioli, F., d'Ovidio, F., Doglioli, A., and Petrenko, A. (2011). Surface coastal
circulation patterns by in-situ detection of Lagrangian coherent structures. Geo-
physical Research Letters, 38(17).
[17] Olascoaga, M. J. and Haller, G. (2012). Forecasting sudden changes in environ-
mental pollution patterns. Proceedings of the National Academy of Sciences.
[18] O'neill, B. (2006). Elementary dierential geometry. Academic press.
[19] Onu, K., Huhn, F., and Haller, G. (2015). LCS Tool: A computational platform
for Lagrangian coherent structures. Journal of Computational Science, 7:2636.
[20] Ourmieres, Y., Mansui, J., Molcard, A., Galgani, F., and Poitou, I. (2018). The
boundary current role on the transport and stranding of oating marine litter: The
French Riviera case. Continental Shelf Research, 155:1120.
[21] Peacock, T. and Dabiri, J. (2010). Introduction to focus issue: Lagrangian co-
herent structures.
[22] Peacock, T. and Haller, G. (2013). Lagrangian coherent structures: The hidden
skeleton of uid ows. Physics today, 66(2):4147.
[23] Peng, J. and Peterson, R. (2012). Attracting structures in volcanic ash transport.
Atmospheric environment, 48:230239.
[24] Shadden, S. C., Lekien, F., and Marsden, J. E. (2005). Denition and prop-
erties of Lagrangian coherent structures from nite-time Lyapunov exponents in
two-dimensional aperiodic ows. Physica D: Nonlinear Phenomena, 212(3-4):271
304.
56
University of Ghana  http://ugspace.ug.edu.gh
[25] Shadden, S. C., Lekien, F., Paduan, J. D., Chavez, F. P., and Marsden, J. E.
(2009). The correlation between surface drifters and coherent structures based on
high-frequency radar data in Monterey Bay. Deep Sea Research Part II: Topical
Studies in Oceanography, 56(3-5):161172.
[26] Tallapragada, P., Ross, S. D., and Schmale III, D. (2011). Lagrangian coherent
structures are associated with uctuations in airborne microbial populations. Chaos:
An Interdisciplinary Journal of Nonlinear Science, 21(3):033122.
[27] Voth, G. A., Haller, G., and Gollub, J. P. (2002). Experimental measurements
of stretching elds in uid mixing. Physical review letters, 88(25):254501.
[28] White, R. (2008). Toxic cities: Globalizing the problem of waste. Social Justice,
35(3 (113):107119.
57