J. Ocean Univ. China (Oceanic and Coastal Sea Research)https://doi.org/10.1007/s11802-021-4493-x
ISSN 1672-5182, 2021 20 (4): 790-802
http://www.ouc.edu.cn/xbywb/
E-mail:xbywb@ouc.edu.cn

Adjoint Method-Based Algorithm for Calculating the Relative
Dispersion Ratio in a Hydrodynamic System
JI Fei1), JIANG Wensheng2), *, and GUO Xinyu1), 3)
1) State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, Ministry of Natural
Resources, Hangzhou 310012, China
2) Key Laboratory of Marine Environment and Ecology (Ministry of Education of China), Ocean University of China,
Qingdao 266100, China
3) Center for Marine Environmental Studies, Ehime University, Matsuyama 790-8577, Japan(Received February 19, 2020; revised December 8, 2020; accepted April 9, 2021)© Ocean University of China, Science Press and Springer-Verlag GmbH Germany 2021
Abstract Relative dispersion ratio (RDR) can be used to quantify the deviation behavior of a water parcel’s trajectory caused by adisturbance in a hydrodynamic system. It can be calculated by using a standard method for determining relative dispersion (RD),which accounts for the growth of the deviation of a cluster of particles from a specific initial time. However, the standard method forcomputing RD is time consuming. It involves numerous computations on tracing many water parcels. In this study, a new method
based on the adjoint method is proposed to acquire a series of RDR fields in one round of tracing. Through this method, the continuous variation in the RDR corresponding to a time series of the disturbance time t can be obtained. The consistency and efficiency ofthe new method are compared with those of the standard method by applying it to a double-gyre flow and an unsteady ArnoldBeltrami-Childress flow field. Results show that the two methods have good consistency in a finite time span. The new method has anotable speedup for evaluating the RDR at multiple t.
Key words

relative dispersion; particle tracking; adjoint method; computational efficiency

1 Introduction
Particle tracking is widely used in oceanographic studies involving either in situ observation or numerical simulation. Its applications include rescue, marine oil spill detection, and marine environmental research (Hackett et al.,2009; García-Garrido et al., 2015; Kim et al., 2015). A
good example of applying water parcel tracking is theautonomous Lagrangian current explorer float, which also
plays an important role in ocean observations, such as theglobal ARGO program (Davis et al., 1992). Particle tracking is also used to find the source of certain substances
distributed in seawater through backward tracking in various areas, including coasts and oceans (Durgadoo et al.,2017; Gelderloos et al., 2017).
However, with particle tracking, a small disturbance onparticles (either temporal or spatial) may cause a large
path deviation. Particles released almost simultaneously atnearly the same location will exhibit dispersion behavior(Batchelor, 1952; Kuznetsov et al., 2002; LaCasce, 2008).
Particles are normally assumed to be moving with theflow; as such, dispersion reflects the character of the flow
* Corresponding author. Tel: 0086-532-66782977
E-mail: wsjang@ouc.edu.cn

field. Thus, the trajectory of a single particle is insufficient to represent the complete information about masstransport in seas. As a result, particle dispersion, which is
caused by some mechanisms such as mixing in a hydrodynamic system, is also defined to characterize materialtransport and advection velocity (LaCasce and Ohlmann,
2003; LaCasce, 2008).
To quantify the dispersion effect during particle tracking, Batchelor (1952) proposed the concept of relative dispersion (RD). LaCasce (2008) provided a simpler form thatis described by releasing a cluster of particles in the vicinity of x0 at the initial time t0. The RD of point t0 at time T
with disturbance at t0 is defined in Eq. (1) as the varianceof the displacement of that cluster of particles

D T , x0 , t0  =1
1ii,j j N x T , xi 0 , t0   x T , x j 0 , t0
N ( N  1)

 , (1)2where N is the total number of particles in a release experiment, t0 is the releasing time, xi0 (1 ≤ i ≤ N ) is the releasing position of the ith particle at the initial releasingtime t0 in the vicinity of x0. The operator || || is the norm in
an n-dimensional Euclidean space and denotes the separation of particles from the center of mass (Batchelor,
1952). x(T, xi0, t0) is the position of the particle at time T

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

released at xi0 at the initial time t0. If xi0 and t0 are kept
constant and only T varies, x(T, xi0, t0) represents the trajectory of the particle being released at position xi0 at time
t0. T, xi0 and t0 are commonly separated with either ‘|’ or‘;’ to reflect the Lagrangian characteristics and emphasize
the varying feature of the variable T and the fixing featureof the particle position xi0 at time t0 in defining the trajectory. T. In this study T, xi0 and t0 vary together, so neither ‘|’ nor ‘;’ is used to avoid confusion.
RD is an important index used to describe the characteristics of particle separation. First, the value of RD represents the repelling (or attracting) feature of particles in aflow field. Second, the change in RD over the time elapsed
after a disturbance may reflect different dispersion mechanisms (Kraichnan, 1966; Bennett, 1984). For example, if
dispersion is dominated by eddies with the same scale as
the particle’s separation distance, the growth of RD is inaccordance with the power law. If dispersion is dominatedby large-scale circulation, RD increases exponentially intime. If the separation between particles is larger than thatin a circulation system, the motion of each particle shows
an independent behavior that linearly increases RD withtime. If the dispersion is caused by a random walk, RD isconstant. Lastly, RD shows the sensitivity of a particle’s
trajectory to its initial state, i.e., it is an index of the firstkind of predictability (Lacorata et al., 2014).
D(T, x0, t0) depends on the initial distances among particles in a cluster. To eliminate the influence of the initial
separation, we define the normalized RD as the relative
dispersion ratio (RDR) as ,expressed as D(t, x0, t0)/D(t0,x0, t0), which is used to represent the dispersion feature of
a flow field. The physical properties of RD (such as representing the repelling or attracting the features of particles and reflecting the mechanism of dispersion) can be
inherited by the RDR. According to the form of the RDR,
it is also an important intermediate variable, i.e., the squareof the amplification factor of the distance between the particle pairs in the calculation of the finite time Lyapunov
exponent (FTLE; Pierrehumbert and Yang, 1993).
The RDR or the derived variable FTLE has been applied to many oceanography problems. Sanderson (2014)
verified the interface between a river and an ocean in an
estuary in terms of the RDR calculated from a 3D flow
field. Cucco et al. (2016) used RDR to indicate the influence of wind data on the predictability of sea surface
transport, which is crucial to maritime search and rescue
operations. The spatial structure of the RDR is also utilized to explain water exchange in seas. For example,
Fiorentino et al. (2012) applied this method and explained
the mechanism of sewage accumulation in coast water.
Using high-frequency radar data and high-resolution sea
surface drifter data, Shadden et al. (2009) studied the exchange between water masses and calculated FTLE to verify the barrier effect of a high-FTLE ridge.
In observations and numerical experiments, the algorithm for calculating RD or RDR consists of two steps.
First, a cluster of particles around one specific location islabeled at the initial time and each particle’s exact position is retained. Second, the RD in Eq. (1) is calculated by

791

recording the displacement of the particles in the cluster
as time proceeds. This ‘forward problem’ framework cancapture a series of subsequent effects of a single initial
disturbance, but it needs a considerable computation effort. To find the RDR at one point, at least two particles
are tracked through integration. Normally, RDR spatial
coverage is required to analyze the dispersion behavior of
a flow field, and a considerable number of particles areneeded to be deployed in the study area. Consequently, thisprocess requires numerous computation tasks. Furthermore, many time instances or time spans, namely t0 and T,
are required for the study, thereby enlarging the computation tasks.
In this study, a new method for calculating RDR is proposed by using the adjoint method to improve the computation efficiency of RDR. The adjoint method is a general approach to solve an inverse problem, which passes
the derivative of one variable over state variables or someparameters backward to the whole process of a dynamic
system (Le Dimet and Talagrand, 1986). In ocean studies,
the adjoint method is widely utilized in sensitivity analysis, data assimilation, parameter estimation, optimal observation network design, and stability analysis (Errico,
1997; Daescu and Langland, 2013).
In this study, the RDR derived via the new method iscalled Adj-RDR to be distinguished from the RDR basedon the traditional method. In one round of calculation, aseries of Adj-RDR fields can be obtained on the basis ofparticle dispersions at the final time T but with different
initial time t. Thus, the potential disturbance at each time
instance t to the final position of the particle at time T isdetermined. The new method is applied to two ideal flow
fields to verify its consistency and efficiency.
The remaining parts of the paper are organized as follows. The existing RDR acquisition methods based on thedefinition and flow mapping are introduced in Section 2.
The derivation of the new algorithm based on the adjoint
method is also described in Section 2. Two ideal flow field
experiments for the verification of our new method aredemonstrated in Section 3. The results are discussed in
Section 4, and the conclusions are presented in Section 5.

2 Methods
In Section 2.1, the standard RDR calculation method ispresented to explain the new method. An algorithm basedon flow mapping (Haller, 2002) is also proposed to extend the RDR as a continuous spatial function of the initial state. In Section 2.2, the new algorithm based on theadjoint method for calculating RDR is introduced.

2.1 RDR Calculation Based on a Standard Algorithm
The trajectory of a particle in a flow field is controlled
by the following equation:x u( x(t, x , t ), t ) ,
t

(2)

where u( x (t, x , t ), t ) denotes the velocity of the particle

792

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802( x , t ) at time t . Its solution can be formally written asxT  x T , x , t   x   u  x  t, x , t  , t  dt .
T

(3)

t

Eq. (3) can be regarded as the mapping of x , the initial
position of particle ( x , t ) at time t into its final positionxT at time T, which can be written as
Tt : x  x T .

(4)

In Eq. (1), RD is calculated on the basis of the relative
displacement among N particles. It involves the distance
between N(N−1)/2 pairs of particles. To simplify the derivation without losing generality, we define the distance ofan arbitrary pair of particles as

d i , j T , x0 , t0   x T , xi 0 , t0   x T , x j 0 , t0 ,2

so that D(T, x0, t0) is the arithmetic mean of di,j(T, x0, t0)
and defined as
D  T , x0 , t 0  = T , x0 , t0 

 i  j d i , j  t0 , x0 , t0 Td T , xt , t    xt  Tt  xt   Tt  xt  xt .
T

 xt   1

 1  0 T  n       1   n  , 0  n (10)

If x is chosen near the point on the undisturbed trajectory xt = x(t, x0, t0) when x0 and t0 are fixed, then the twoparticles ( x , t ) and (xt, t) can form a particle pair, and theirinitial displacement is xt  x  xt . At the final time T,
their displacement isx T , xt , t   x T , x , t   x T , xt , t  x T , x , t   x T , x0 , t0  ,which can be expanded at xt asx T , xt , t   Tt  xt  xt  r  xt  ,xt  0(9)

The eigenvalue decomposition of CtT0 ( xt ) can be expressed as

1i , j  N

where r(∆xt) is the residual term with lim(8)d(T, xt, t) changes with the norm of disturbance ||∆xt||
and the direction of disturbance ∆xt/||∆xt||. Therefore, R(T,x0, t0) is a statistical variable, although it is normalizedwith the arithmetic average of di,j(t, x0, t0) = ||∆xt||2.
Haller (2002) provided a definition with a higher certainty by using the Cauchy-Green strain tensor at xt

CtT0

.

(7)

where [ ]T denotes the transposition of a vector or a matrix.
When the higher-order infinitesimal is omitted, the RDof a particle pair with the initial distance ∆xt can be expressed asT

According to this definition, the RDR is expressed in
the form of
R  T , x0 , t0  Td T , x0 , t0    x T , x0 , t0    x T , x0 , t0  ,CtT0  xt   Tt  xt   Tt  xt  .

1
 i  j d i , j  T , x0 , t 0  .
N (N  1)
1i , j  N

i  j di, j
1i , j  Nparticles near xt.
Incorporating Eq. (5) into the matrix form of Eq. (1) obtains(5)r  xt xt

0and the Jacobian of the mapping is x1 T , xt , t 
x1 T , xt , t  



,
x
,

x
t
t
xn  t , xt , t   1 t ,Tt  xt    xn T , xt , t xn T , xt , t  xn  t , xt , t   x1  t , xt , t (6)where Tt ( xt ) describes the deformation rate of the infinitesimal perturbation of the initial position xt or thedevelopment rate of the distance between two adjacent

where the eigenvalue λi indicates the extension rate of the
distance between the particles in a pair in the direction
denoted by the eigenvector ξi.
Dispersion is led by the maximum eigenvalue in the
direction of the corresponding eigenvector; as such, themaximum eigenvalue is often used as the index of thesquare of the amplification factor of the distance between
the particles in a pair in the calculation of FTLE (Haller,
2002). In our study, the same notion is applied to calculate RDR:
R(T , xt , t )  max(1 , 2 , , n ) .

(11)

The calculated R(T, xt, t) corresponds to the maximum
RDR among all the initial displacement directions. In thisstudy, the method described in this section is called thestandard method.

2.2 RDR Calculation Based on the Adjoint Method
The above deduction reveals that if Tt ( xt ) is obtained, then the RDR along the trajectory can be acquiredon the basis of dispersion from t to T. However, in thestandard RDR calculation described in Section 2.1, theindependent integral of Eq. (3) is required to determine
Tt ( xt ) for different values of t. This procedure needs
numerous computation tasks to calculate the RDR forcontinuous t. Therefore, a new method based on the adjoint method is proposed to simplify the calculation.

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

First, a set of cost functions Jk, k = 1, 2,···, n is defined

as follows:J k  x, x   sign  xk T , x , t   xk T , x0 , t0      t  T  kT   x  t, x , t   x  t, x0 , t0   dt ,Ttwhere γk = [δk1, δk2,···, δkn]T with δki being the Kronecker
delta,  (t  T ) is the Dirac delta function, and sign( ) is
the function that takes the sign of the variable with theabsolute value of 1. It represents the distance between thetwo particles (x0, t0) and ( x , t ) in the kth dimension attime T which can be in its simple form:J k  x, x   xk T , x , t   xk T , x0 , t0  .

 x Lk ( x, xk , x )   x J k ( x, x ) .

(13)

 x J k ( x, x )  sign(xk (T , x , t )  xk (T , x0 , t0 )) Tt  (t  T )γ k  x x (t, x , t )dt .

T
 x

J k  x, x    [ xk ]T    u  x  t, x , t  , t   dt , (14)
t
 t






T

T

t

t

(15)

If x  xt , then

Lk ( x, xk , x ) 



(12)

 T
] is the vector of the Lawhere xk (t )  [ xk1 , xk 2 , , xkn
grangian multiplier and k = 1, 2,···, n.
x  t, x , t  is the trajectory of particle ( x , t ) , as such, thesecond term of Lk ( x , xk , x ) in Eq. (14) is always 0 because of Eq. (2). Therefore, we have

Accordingly, a set of functionals with the constraints canbe defined as

Lk x , xk , x  J k  x T  xk

793

T

(16)

Eq. (14) is integrated by parts to obtaink

 x  x    [ x]  xt  [ x ]  u  x t, x , t  , t dt
T


k

T



 T
k




T
x 
T
 J k   x (T , x , t )   xk (T )  [ x ]T  xk (t )    x T  k  xkT  u  x  t, x , t  , t  dt .
t

t




(17)

When the Jacobi operator is applied to Eq. (17), and Eq. (16) is inserted, the following equation is derived:





T

x Lk x, xk , x   xk T   x x T, x , t   xk  t   x   TTsign  xk T , x , t   xk T , x0 , t0     k    x k    x u  x  t, x , t  , t     t  T  γ kT   x x  t, x , t  dt , (18)t  t





T

where

Comparing the above equation with Eq. (6), we can express the Jacobian matrix of the flow mapping Tt ( x )as the gradient of the cost function,  x J k :x k  sign  xk T , x , t   xk T , x0 , t0   xk .
If we setTT

 x k  TT    [ x k ]   x u  x  t, x , t  , t     t  T  γ k , t (19)

with xk T   0 , i.e., x k  0 , then Eq. (19) and its initialcondition define a set of the adjoint functions of x k , so
Eq. (18) becomes x Lk x , xk , x   xk  t  .
Eq. (15) is transformed to x J k ( x , xk , x )   xk (t ) .

(20)

From Eq. (13), we can derive x J k  x (T ,x , t ) xk (T ,x , t ) sign  xk (T , x , t )  xk (T , x0 , t0 )  k.
xn (t ,x , t )  x1 (t ,x , t )Ttsign(x1 (T , x , t )  x1 (T , x0 , t0 )) x J1T .
 x   T sign(xn (T , x , t )  xn (T , x0 , t0 )) x J n 
(11)

By inserting Eq. (20) into Eq. (21) and considering thedefinition of x k we haveTt ( x )    x1 (t ), , x n (t )  .

(22)

In Eqs. (20) and (22), x cannot be xt, but we can evaluate Tt ( xt ) with Tt ( x ) when x is approaching xt.
When Eqs. (19) and (22) are examined, no apparent restriction is imposed on x  xt ; thus, Tt ( xt ) can be expressed asTt ( xt )    x1 (t ), , x n (t )  ,(23)

where x k is obtained from Eq. (19) when x is set to xt
in that equation. Then, the RDR can be calculated with
Eq. (11).

794

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

As an adjoint equation, Eq. (19) is numerically integrated from t  T to t  t0 with its initial conditionx k T   0 . The integration time range [t0, T] is split into
M-1 sections with ∆t = (T − t0)/(M − 1) .
Eq. (19) is separately integrated from tm = (m − 1)∆t + t0
to tm+1 = m∆t + t0 separately as follows:
tm1

t

m

  x   TT  k    x k   x u  x  t, x , t  , t  dt    t 



tm1

t

m

  t  T  γ kT dt .

3) Calculate Tt ( xt ) according to Eq. (26) based on
the information in step 2.
4) Calculate CtTm ( xtm ) according to Eq. (9) and obtain
RDR based on Eq. (11).
However, xtm will not always be on the grid points similar to those at the beginning, and some particles may evenleave the study area. The interpolation is also used to get
RDR on each grid point.

3 Application of the Proposed Method
in Two Ideal Flow Fields

In particle tracking studies, two common questions arecommonly encountered. One is, ‘what is the final particle
distribution that can occur due to the disturbance at time
(named as disturbance time t0)?’. The other one is, ‘whatkind of disturbance at t0 can cause the final distribution ofparticles?’. The first one is called the forward tracking
problem, which can be used to describe the spreading of
T
m
m 1
m
m
x k   x k    x u    x k   γ k  mM ,pollutants at a certain time (named as result time T) and

being discharged at different t0, e.g., different phases of am  1, 2, , M , (24)tidal period. The second one is called the backward tracking problem, which can be utilized in instances such aswhere x k( m )  x k ((m  1)t  t0 ) is the value of x k delocating the oil spill source at t0 from the oil coverage in
fined at discrete temporal points, tm = (m − 1)∆t + t0, u(m) =water at T. In both cases the final time T is normally reu(x(tm, xt, t), tm) , and δmM is the Kronecker delta. xt is thegardedas a fixed value, and t0 is assumed as a varying
position of the particle (x0 , t0) at time t , so u(m) can alsovariable. The RDR field at a continuously varying t0 is
be expressed as u(m) = u(x(tm, x0, t0), tm).
required to solve these problems.
If Eq. (24) is written in a matrix form under the initial
The RDR in the two problems can be obtained by usingM
condition of x k   0 , we havethe traditional method described in Section 2.1. They canalso be more efficiently solved by applying the new methT
 x  m 1 , , x  m 1    I   u m   t   x  m  , , x  m   od presented in Section 2.2. In this section, two widely
x
1
1
n
n



 




used analytical flow fields, namely, two-dimensional doum < M, (25)ble-gyre flow and three-dimensional unsteady ABC flow,
are adopted to verify the validity and efficiency of the
M 1
M 1
and  x1  , , x n     I .
new method proposed in Section 2.2. These flow fields
are benchmark models for studying the Lagrangian chaosThen, inserting Eq. (25) into Eq. (23), we can have theand mixing either 2D (Shadden et al., 2005; Coulliette
Jacobian of the mapping defined in Eq. (4) at different
et al., 2007) or 3D dynamic systems (Dombre et al., 1986;
time steps tm along the trajectory, i.e.,
Haller, 2001). In the present study, two-dimensional double-gyreflow and three-dimensional unsteady ABC flowm 
i   T 
T
are
used
in forward and backward tracking problems, re tm xtm    I   x u
t  .
(26)
spectively. The RDRs of these two flow fields at differenti  M 1 
initial times are calculated by using the standard method
and the newly proposed adjoint method. The results of
The RDR expressed in Eq. (11) can be obtained based
both methods are compared to verify the effectiveness ofon the Jacobian of flow mapping in Eq. (26). In this way,
the new method.
the influence of the disturbance at a series of time inFor the convenience of reference, the RDRs calculated
stance tm to the final location of the particle can be detervia
the new method and the standard algorithm are called
mined efficiently by evaluating the gradient of the flow
Adj-RDRand standard RDR, respectively.
field along the trajectory.
The left-hand side is discretized with an implicit timeforward scheme. The right-hand side can be integrated analytically by considering a Dirac function in the integrand.
Thus, the numerical scheme of Eq. (19) is expressed asfollows:

 

The calculation of RDR via the adjoint method can bedivided into four steps:
1) Place a numerical particle on every grid point in the
study area at a certain initial time.
2) Track each particle in the flow field from the initialtime to the final time T, and keep the position of the particle at time tm, xti and xu(i) at the same point at thesame time, where i = 1, 2,···, M.

3.1 Double-Gyre Flow CaseIn this section, the forward tracking problem is solved
for a double-gyre flow field. This flow field is designed
as two adjacent closed gyres flowing at a constant speed.
They are separated by a straight line named as the flow
axis, which oscillates periodically with a small amplitude
so that the two gyres can expand and contract alternately.

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

The flow field is given byu1   πAsin(πf ( x1 , t ))cos(πx2 )f ( x1 , t ) ,u2  πAcos(πf ( x1 , t ))sin(πx2 ) x

1

where
f ( x1 , t )  a(t ) x12  b(t ) x1 ,a (t )   sin(t ) , b(t )  1  2 sin(t ) .
In this study, the parameters are set as ω = 2π/10, ε = 0.1,
A = 0.1, i.e., the period of the oscillation of the flow axis
is 10 with an amplitude of 0.1 in the x1 direction, and the
speed on the flow axis is 0.1, all variables in the ideal
flow field are dimensionless. The calculation area is set as[0, 2] in the x1-axis direction and [0, 1] in the x2-axis direction. Resolution is set as ∆x1 = ∆x2 = 0.002 so that the795

grid number is 1000×500. The time step for obtaining the
trajectory by integrating Eq. (3) is 0.002. In this study, a
4th-order predictor and a 5th-order corrector scheme areused to determine the integrals.
For calculating the RDR with the standard algorithm,
particle pairs are released on the grid points with the initial displacement of 10−5. Six groups of experiments are
set at t0 = 0, 4, 8, 12, 16, 20 and T = 20. The results are
displayed in column A of Fig.1, with the most remarkablefeature of the double-gyre pattern. The high RDR occurs
near the flow axis and forms several bands that stretch
like a circle. The influence of the different values of t canbe found in column A in Fig.1 from rows 1 to 6. The highvalue band between the two gyres moves left and rightbecause of the oscillation of the flow axis. The integral
time length also alters the pattern to some extent. The
number of the high-RDR bands and the RDR in the band
increase as the integral time is extended.

Fig.1 RDR of the double-gyre flow calculated by two methods. Column A, RDR calculated by the standard algorithm displayed as log10RDR; Column B, RDR calculated by the adjoint method, displayed as log10RDR; Column C, The absolutedifference between RBR calculated by two methods. Rows 1 to 6 denote the results with initial disturbance time at 0, 4, 8,
12, 16 and 20.

796

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

Adj-RDR is calculated in accordance with the procedure described in Section 2.2 at a time step of 0.02. The
reason for taking this value of ∆t is explained in the discussion. Six figures at t0 = 0, 4, 8, 12, 16, 20 are selected
and presented in column B of Fig.1, which shows theRDR distribution with the same meaning as those in column A. However, the RDRs of different t0 are calculatedby tracking the particles in one run. In Fig.1, the spatial
pattern of Adj-RDR is consistent with the result of thestandard algorithm. Their difference is less than 15% of
the absolute value in column C of Fig.1, but Adj-RDR ishigher than the result obtained via the standard algorithm.
The overestimation of Adj-RDR occurs in the low-value
region of RDR, e.g., the center of the gyres. The difference between the two algorithms in a high-value area,
where we are usually most interested in, is smaller thanthat in a low-value area. Therefore, the result of Adj-RDR
is acceptable in the double-gyre test.
The area-averaged RDRs obtained with both methods
are also compared in Fig.2. Adj-RDR is calculated at atime step interval of 0.02, so it is represented as a continuous line. The RDR based on the standard algorithm isshown as 6 points in the figure. When t0 > 8, i.e., the track
span is shorter than 12, Fig.2 shows a good correspondence between the two RDRs based on different algorithms. The value produced via the new methodis higher
than that obtained via the standard method. If t0 < 8, i.e.,
the track span is longer than 12, then Adj-RDR exponentially grows with the track span, but the standard RDR
appears to reach the maximum value at 107.
The computation time consumed by using the two methods is compared. A desktop computer with a 3.6 GHz AMD
3700x processor is used to calculate the RDR with MATLAB 2017b. ‘Tic’ and ‘toc’ commands are utilized to keep
time. The computation time is the average of the results ofthree rounds of calculations. The standard method takes
138.1 ± 0.5 s to obtain the distribution of the RDR in six
cases (t0 = 0, 4, 8, 12, 16, 20). By comparison, the adjoint
method needs 136.0 ± 0.5 s to determine the distribution of
the RDR in 100 cases at t0 of 0 to 20 with a time step of
0.2.

Fig.2 Area averaged RDRs calculated by both methods continuous line, Adj-RDR; point, RDR by standard algorism.

3.2 Backward Tracking in an Unsteady ABC Flow
As mentioned at the beginning of this section, the RDR
can be used to determine the source of matter in oceansby solving the backward tracking problem. In the application of backward tracking, a new flow field is formulated
by reversing the time sequence and taking the opposite
sign of the original flow field. Thus, a new flow mappingis defined as
Tt : xT  x .

(27)Then, the RDR can be calculated for this flow mappingby using the standard algorithm or the adjoint method as
the forward tracking problem.
The physical definition of the derived RDR differs fromthat obtained through forward tracking, as in Section 3.1.
The high RDR, in this case, indicates the area where particles gather (Haller, 2015).
In this section, the unsteady ABC flow is used as thetest flow field, which is given by
1u1   A  2 t sin(πt )  sin( x3 )  C cos( x2 )




1


u2   A  t sin(πt )  cos( x3 )  B cos( x1 ) .
2u3  C cos( x2 )  B cos( x1 )In this study, the parameters are set as follows: A = 3 , B   2 , and C = 1. The computational domain is
in the cube {x = (x1, x2, x3)| xi ∈ [0, 2π], i = 1, 2, 3} which
is divided into 100×100×100 grids. The track span is set
as 8, which is four periods of motion.
In the backward tracking problem, t varies from t0 = 8
to T = 0. It can be transformed into a forward tracking problem by assuming that t  t0  T  t , and the new velocity
field is as follows:1u1    A  2 (8  t )sin(π(8  t ))  sin(x3 )  Ccos(x2 )




1


u2    A  (8  t )sin(π(8  t ))  cos(x3 )  Bcos(x1 ) .
2
u3  Ccos(x2 )  Bcos(x1 )The original backward tracking problem here is to findthe disturbance at t that can generate the end state at T. In
the transformed forward problem, the goal is to obtain the
RDR when t is from 0 to 8 and the final state is T  8 .
When the standard algorithm is used t = 0, 0.4, 0.8, 1.2,
1.6, that corresponds to t = 8, 7.6, 7.2, 6.8, 6.4 is chosen.
When the new adjoint method is applied, t is selectedfrom 0 to 8 with a time step of 0.02.
The RDR distribution on the three surfaces of the calculation cube is shown in Fig.3 when t  0 , i.e., t = 8 ischosen. The results from the standard algorithm and the

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

adjoint method are presented in Figs.3a and 3b, respectively. They display the ‘hollow’ pattern of a high value
around the vortex tube. However, the RDR obtained withthe adjoint method has a ‘higher contrast’ than that determined with the standard algorithm. Therefore, Adj-RDR
is higher than the standard RDR in the high-value zone,
whereas the former is lower than the latter in the lowvalue zone. This pattern can also be found in Fig.3c, whichshows the difference between the RDRs acquired by bothmethods.
The RDR distribution at different initial times on thethree slices in the middle of the cube, i.e., x1 = π, x2 = π,and x3 = π, is plotted to express the results more clearly.
The RDRs calculated by the standard algorithm and the
adjoint method are displayed in Figs.4 and 5, respectively.
Both RDR distributions are similar. However, the RDR in

797

Fig.5 has a ‘higher contrast’ than that in Fig.4. Furthermore, the RDR distribution in Fig.5 becomes blurrier as
the time span between t0 and T is prolonged.
The time series of area-averaged Adj-RDR and standard
RDR are shown in Fig.6. Adj-RDR exponentially grows as
t increases. By comparison, the standard RDR grows exponentially when t < 6; as it reaches the maximum value, it
no longer grows as t increases.
The time consumed by using the two methods in the
ABC flow is compared. It is a 3D flow field, so the computation time is very long. For this reason, a parallel MATLAB m file is written and run in 10 threads. With the standard method, obtaining the RDR field in 1 case with t = 8
takes 1331.3 s. With the adjoint method, determining the
distribution of the RDR field in 400 cases with t from 8 to
0 needs 2354.7 s.

Fig.3 RDR distribution of unsteady ABC flow field on three surfaces of the cube. (a), RDR calculated by the standard algorithm as log10RDR; (b), RDR calculated by the adjoint method as log10RDR; (c), the absolute difference between RBR calculated by two methods time range is for T = 0 and t0 = 8.

Fig.4 RDR distribution of unsteady ABC flow field calculated by the standard algorism on three inner sections of the cubewith different initial time. The first row denotes x1 = π; second row denotes x2 = π; third row denotes x3 = π. The columns
from 1–5 correspond to the initial time t = 8, 7.6, 7.2, 6.8 and 6.4. The color in the figures denotes the value log10RDR.

798

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

Fig.5 Same as Fig.4, but calculated by adjoint method.

In the standard method, particle tracking is the mosttime-consuming procedure. In the adjoint method, particle tracking and the calculation of Eq. (26) are both time
consuming. Either particle tracking or particle tracking
plus Eq. (26) consumes more than 90% of the calculation
time of the corresponding method in the two examples in
Section 3. Therefore, we focus on the comparison of thecomputational complexity of particle tracking and the calculation of Eq. (26).
The calculation cost of tracking one particle in an nd dimensional space in one time step needs at least O(nd2)
multiplications when a simple discrete form of Eq. (2), i.e.,
Eq. (28), is used to track the particle:x ( m 1)  I   x u( m 1) t
Fig.6 The time series of area-averaged Adj-RDR and standard RDR. The red line denotes the continuous time seriesof the sum of Adj-RDR change with different initial time
and the black dots denotes the sum of RDR calculated bystandard algorithm. Each dot is averaged from an independent experiment with the disturbed particles released
at t.

4 Discussion
4.1 Computational Efficiency
The RDR field calculated with the adjoint method is
almost continuous with respect to. By contrast, the RDR
field obtained with the standard method is confined to acertain instance of t. The calculation procedure needs tobe carried out several times to determine the influence ofdisturbance from time t to T. Therefore, the new method
is beneficial. The computational complexity of the two
methods is briefly analyzed in this section.

 x
T

(m)

.

(28)

Every calculation of Eq. (28) is assumed as one unit ofcalculation here. By using the standard method, nd + 1
particles should be tracked in the nd -dimensional space todetermine the RDR by calculating Eq. (6) at one point. If
nt time steps are needed to track a particle from t0 to T,
then the computation cost to obtain R(T, x0, t0) is nt × (nd+ 1). The computation cost is nt × (nd + 1) × nr for the disturbance at a single time instance t0 to determine the RDR
field with nr sampling points. If a series of RDR fields isneeded with respect to nc disturbance time instances, e.g.,
the 6 cases in Section 3.1, nc independent calculations withdifferent tracking time steps are needed. The disturbancetime necessary. t is spread evenly between t0 to T for nc
cases. Then, the time step number in each case decreases
from nt to 0, while t is from t0 to T. Therefore, the total
computation time to complete the whole calculation isnc/2 × nt × (nd + 1) × nr.
In the adjoint method, the computation cost includes

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

the computation of Eq. (26) and the computation of particle tracking. The computation cost of Eq. (26) is nd timesof Eq. (28). Particle tracking should be calculated from t0
to T (nt steps) and the calculation of Eq. (28) from T to t0(nt steps) to acquire the RDR by using the adjoint method
at one point. Therefore, for the continuous time series of
Adj-RDR field at nr points, the computation cost is nt ×
(nd + 1) × nr.
Theoretically, the adjoint method saves the computation cost nc/2 times of that of the standard method in solving the inverse problem. In real calculation, the adjoint
method needs to record and read the data of particle position and velocity, so the calculation time tends to be longer
than the theoretical value. The speedup is less than nc/2,
but it is still very notable. In this study, the speedup ratioof the double-gyre example is about 16.9 times, which isless than 50 times in theory. For the unsteady ABC flow,
the speedup ratio is about 113 times, which is also less than
200 times in theory.
In addition to time cost, the memory cost of the adjoint
method is determined to calculate the RDR. The extra random-access memory (RAM) occupied by the Adj-RDR
method mainly occurs in step 3 of the calculation procedure. In step 3, memory needs to save the spatial divergence tensor of the velocity for each particle of the current time step. The tensor is d times larger than the velocity vector, where d is the spatial dimension; as such, theRAM usage of step 3 is also d times of the particle tracking process. For the double-gyre experiment in Section
3.1, the tensor takes up 15.305 megabytes of RAM.
Apart from the RAM usage, extra disk space is requiredin the Adj-RDR method to record the trajectory information compared with that in the standard method. The data
size of the trajectory can be expressed as follows:799

in the smallest neighborhood along the particle trajectory.
In the standard algorithm, the initial pair is finally separated into two tracks regardless of the size of the initial
disturbance because of the effect of the dispersion. Whenthe particle pairs are apart far enough, the relative motion
between the particle pairs is basically random. Statistically, the average distance is no longer growing. Ding andLi (2007) explained the saturation phenomena in the study
on the nonlinear FTLE. If the initial disturbance is small,the standard RDR is consistent with Adj-RDR in a longer
time (Fig.7).
The RDR obtained by the adjoint method is sensitive tothe time step in the calculation because of the tangent linear approximation. A set of experiments of the double-gyreflow field (as in Section 3.1) is used to discuss the relationship between the time step ∆t and the Adj-RDR. A
standard algorithm with the initial disturbance δx = 10−5 isused as a benchmark result. Figs.8 and 9 illustrate the influence of different time steps on the difference in the two
methods. The smaller the time step, the smaller the difference in the RDR between the two methods. Moreover,in Fig.9, the difference increases when ∆t = 0.032. Thus,
when ∆t is larger than the upper limit, the difference increases rapidly with ∆t. However, when ∆t is less than thetime step, the decrease in ∆t slightly affects the reductionin the difference. In the above experiment, ∆t = 0.02 is used
for the compromise between efficiency and accuracy.
Tangent linear approximation is introduced to Adj-RDR
in the derivation of Eq. (24), and the second and higherorder term of ∆t is disregarded. Therefore, an upper limitof ∆t should be set to ensure that  x u( x (t, xt , t ), t ) does
not change significantly with ∆t. The limit is also related
to the characteristics of the flow field, i.e., the temporalvariation in  x u( x (t, xt , t ), t ) .

File’s data size = Number of time steps × Number of released particles × Dimension of flow filed (float data).
The adjoint method requires a relatively high resolution
in space and time, so the footprint of the trajectory information is a large amount of data. In the case of the double-gyre experiment in Section 3.1, the total size of the
trajectory files is 6.872 gigabytes. In general, the extra
cost of RAM and disk space is acceptable.

4.2 Computation Accuracy
Adj-RDR has the same spatial structure as that calculated by the standard algorithm in forward and backward
tracking applications, as displayed in Section 3. However,
Adj-RDR is higher than the standard RDR, especiallywhen the calculation span is large. In the time series of
the area-averaged Adj-RDR and standard RDR shown in
Figs.2 and 6, Adj-RDR exponentially grows as the calculation time span increases. For the standard RDR, the dispersion ratio grows exponentially during a limited time;afterward, the area-averaged RDR appears to reach the
maximum value. This difference occurs because the adjoint method calculates the tangent linear dispersion ratio

Fig.7 Same as Fig.2, but the results of the standard RDR
with a smaller initial disturbance scale are supplemented.
The red line represents the RDR at each initial time calculated by the adjoint algorithm. The blue, cyan, and yellow dots correspond to the RDR calculated by the standard algorithm, and the initial distances between the particle pairs are 10−5, 10−6, and 10−7, respectively.

800

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

Fig.8 Absolute difference between Adj-RDR and standard RDR (Fig.1A). Columns A, B, C, D, and E denote the results
obtained via the adjoint method with time steps of 0.002, 0.008, 0.032, 0.128, and 0.512, respectively. Rows 1 to 6 denote
the results with initial disturbance times at 0, 4, 8, 12, 16, and 20.

Fig.9 Sum of the absolute difference between Adj-RDR
and standard RDR, and the change with time step.

4.3 Limitations of the New Method
The RDR calculated by the adjoint method is blurrier
in distribution as the time span between t0 and T is prolonged (Fig.5). Therefore, the accuracy of the new method
is limited when the time span is long.
The reason for the time span limitation is that R(T, xt, t)
is calculated on the particle in the adjoint method, i.e.,

RDR is mapped on xt, which varies with t. For comparison, in the standard method, RDR is mapped on the initial
location where the particle is released. In an aperiodic
flow field, when t is farther away from t0, xt is farther
away from the mesh point where the particles are released(Fig.10). In the area near the boundary, some particles
leave the study area. In other words, the number of sampling points in the study area is reduced, and the distribution of sampling points is not uniform when t − t0 is high.
This effect is more significant for the flow field with anopen boundary because particles eventually flow out of
the study area.
However, if a study aims to calculate the structure of
the RDR in a finite time, such as several tidal periods, theadjoint method is valid for the application.

5 Conclusions
In this study, a new method for calculating RDR is
proposed by using the adjoint method. The validity of themethod is verified via two idealized numerical experiments. The results reveal that the spatial structure of RDR
calculated with the adjoint method is consistent with thatobtained with the standard algorithm. The new method is
successful in forward and backward tracking cases. Its

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

801

Fig.10 Adj-RDR distribution of the unsteady ABC flow field mapped on the location of a particle.

computational efficiency is also exhibited with a speedupof more than 15 times in the two examples.
The new method has some limitations. The spatial resolution of RDR decreases as the tracking time is extended.
Moreover, RDR increases as the integration time is prolonged. In many cases, this limitation can be avoided because the RDR of a finite time is needed. RDR calculated
for a long time loses its physical significance because of
the nature of the flow field itself.
In conclusion, the newly developed method for computing RDR via the adjoint method can provide the time
evolution of the dispersion characteristics of a flow field.
The present study helps elucidate the mass or momentum
distribution in seas. In future studies, the new method will
be applied to real cases in oceans.

Acknowledgements
Fei Ji thanks the China Scholarship Council (CSC) for
supporting his stay in Japan and the Ministry of Education,
Culture, Sports, Science and Technology, Japan (MEXT)
under a Joint Usage/Research Center, Leading Academia
in Marine and Environment Pollution Research (LaMer)
Project for supporting his short stay in Japan.

References
Batchelor, G. K., 1952. Diffusion in a field of homogeneous
turbulence II. The relative motion of particles. Mathematical
Proceedings of the Cambridge Philosophical Society, 48: 345362.
Bennett, A. F., 1984. Relative dispersion: Local and nonlocal dynamics. Journal of the Atmospheric Sciences, 41 (11): 18811886.
Coulliette, C., Lekien, F., Paduan, J. D., Haller, G., and Marsden,
J. E., 2007. Optimal pollution mitigation in Monterey Bay
based on coastal radar data and nonlinear dynamics. Environmental Science & Technology, 41 (18): 6562-6572.
Cucco, A., Quattrocchi, G., Satta, A., Antognarelli, F., De Biasio,
F., Cadau, E., et al., 2016. Predictability of wind-induced sea
surface transport in coastal areas. Journal of Geophysical Research: Oceans, 121 (8): 5847-5871.
Daescu, D. N., and Langland, R. H., 2013. The adjoint sensitivity guidance to diagnosis and tuning of error covariance parameters. Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications (Vol. II). Springer, Berlin, Heidelberg,
205-232.
Davis, R. E., Regier, L. A., Dufour, J., and Webb, D. C., 1992.
The autonomous Lagrangian circulation explorer (ALACE).
Journal of Atmospheric and Oceanic Technology, 9: 264-285.

802

JI et al. / J. Ocean Univ. China (Oceanic and Coastal Sea Research) 2021 20: 790-802

Ding, R., and Li, J., 2007. Nonlinear finite-time Lyapunov exponent and predictability. Physics Letters A, 364 (5): 396-400.
Dombre, T., Frisch, U., Greene, J. M., Hénon, M., Mehr, A., andSoward, A. M., 1986. Chaotic streamlines in the ABC flows.
Journal of Fluid Mechanics, 167: 353-391.
Durgadoo, J. V., Rühs, S., Biastoch, A., and Böning, C. W., 2017.
Indian Ocean sources of Agulhas leakage. Journal of Geophysical Research: Oceans, 122 (4): 3481-3499.
Errico, R. M., 1997. What is an adjoint model? Bulletin of theAmerican Meteorological Society, 78 (11): 2577-2592.
Fiorentino, L. A., Olascoaga, M. J., Reniers, A., Feng, Z., BeronVera, F. J., and MacMahan, J. H., 2012. Using Lagrangian coherent structures to understand coastal water quality. Continental Shelf Research, 47: 145-149.
García-Garrido, V. J., Mancho, A. M., Wiggins, S., and Mendoza,C., 2015. A dynamical systems approach to the surface searchfor debris associated with the disappearance of flight MH370.
Nonlinear Processes Geophysics, 22: 701-712.
Gelderloos, R., Haine, T. W., Koszalka, I. M., and Magaldi, M.
G., 2017. Seasonal variability in warm-water inflow toward
Kangerdlugssuaq Fjord. Journal of Physical Oceanography,
47 (7): 1685-1699.
Hackett, B., Comerma, E., Daniel, P., and Ichikawa, H., 2009.
Marine oil pollution prediction. Oceanography, 22 (3): 168175.
Haller, G., 2001. Distinguished material surfaces and coherent
structures in three-dimensional fluid flows. Physica D: Nonlinear Phenomena, 149 (4): 248-277.
Haller, G., 2002. Lagrangian coherent structures from approximate velocity data. Physics of Fluids, 14 (6): 1851-1861.
Haller, G., 2015. Lagrangian coherent structures. Annual Review
of Fluid Mechanics, 47: 137-162.
Kim, J. J., Stockhausen, W., Kim, S., Cho, Y. K., Seo, G. H., andLee, J. S., 2015. Understanding interannual variability in the
distribution of, and transport processes affecting, the early life
stages of Todarodes pacificus using behavioral-hydrodynamic

modeling approaches. Progress in Oceanography, 138: 571583.
Kraichnan, R. H., 1966. Dispersion of particle pairs in homogeneous turbulence. Physics of Fluids, 9 (10): 1937-1943.
Kuznetsov, L., Toner, M., Kirwan, A. D., Jones, C. K. R. T., Kantha, L. H., and Choi, J., 2002. The loop current and adjacentrings delineated by Lagrangian analysis of the near- surface
flow. Journal of Marine Research, 60 (3): 405-429.
LaCasce, J. H., and Ohlmann, C., 2003. Relative dispersion at thesurface of the Gulf of Mexico. Journal of Marine Research,
61 (3): 285-312.
LaCasce, J., 2008. Statistics from Lagrangian observations. Progress in Oceanography, 77 (1): 1-29.
Lacorata, G., Palatella, L., and Santoleri, R., 2014. Lagrangianpredictability characteristics of an ocean model. Journal ofGeophysical Research: Oceans, 119 (11): 8029-8038.
Le Dimet, F. X., and Talagrand, O., 1986. Variational algorithms
for analysis and assimilation of meteorological observations:
Theoretical aspects. Tellus A: Dynamic Meteorology and Oceanography, 38 (2): 97-110.
Pierrehumbert, R. T., and Yang, H., 1993. Global chaotic mixing
on isentropic surfaces. Journal of the Atmospheric Sciences,
50 (15): 2462-2480.
Sanderson, A. R., 2014. An alternative formulation of Lyapunov
exponents for computing Lagrangian coherent structures. In:
2014 IEEE Pacific Visualization Symposium. Yokohama, Japan, 277-280.
Shadden, S. C., Lekien, F., and Marsden, J. E., 2005. Definitionand properties of Lagrangian coherent structures from finitetime Lyapunov exponents in two-dimensional aperiodic flows.
Physica D: Nonlinear Phenomena, 212 (3-4): 271-304.
Shadden, S. C., Lekien, F., Paduan, J. D., Chavez, F. P., andMarsden, 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): 161-172.
(Edited by Xie Jun)