# Closely spaced co-rotating helical vortices: Long-wave instability

A. Castillo-Castellanos1<EMAIL_ADDRESS>S. Le
Dizès1  1Aix Marseille Université, CNRS, Centrale Marseille, IRPHE, Marseille,
France

###### Abstract

We consider as base flow the stationary vortex filament solution obtained by
Castillo-Castellanos et al. (2021) in the far-wake of a rotor with tip-
splitting blades. The cases of a single blade and of two blades with a hub
vortex are studied. In these solutions, each blade generates two closely
spaced co-rotating tip vortices that form a braided helical pattern in the
far-wake. The long-wave stability of these solutions is analysed using the
same vortex filament framework. Both the linear spectrum and the linear
impulse response are considered. We demonstrate the existence of two different
types of instability modes, that are characterized by their symmetry. The
first type corresponds to the local pairing of consecutive turns of the
helical pattern. It is shown to be well described by the instability of an
uniform helical vortex with a core size given by the mean separation distance
of the vortices in the pair. The second type of unstable modes modifies the
separation distance between the vortices in each pair. It seems to be
associated with local pairing but excites significantly higher frequencies.
These unstable modes are also found to spread spatially with a weaker rate.

##  1 Introduction

Rotating blades, such as those of a helicopter rotor or a horizontal-axis wind
turbine, generate concentrated vortices at their tips, transported downstream,
creating a distinctive helical pattern. These helical vortices are associated
with several practical issues actively investigated. One of these issues is
the wake stability with respect to external disturbances, which has practical
relevance since instabilities may accelerate the vortex break-up and the
transition towards turbulent wakes. The interaction between tip vortices and a
solid surface, like a trailing rotor blade, or a wind turbine located
downstream, may cause significant noise, vibration, and fatigue problems. Wake
stability is also fundamental to optimize wind farms design, since wakes may
extend up to 50 km downstream under stable atmospheric conditions (Cañadillas
et al., 2020). Strategies to maximize power generation include dynamically
varying the yaw, pitch, and tip-speed ratio to excite the natural instability
of helical wakes (see, for instance Huang et al. (2019); Frederik et al.
(2020)).

One alternative to accelerate the vortex breakdown is to modify the wake
structure altogether. Brocklehurst & Pike (1994) introduced a modified vane
tip to split the tip vortex into two co-rotating vortices. The associated
spreading of vorticity has been suggested as indicative of reduced noise and
increased aerodynamic efficiency, but there is little information regarding
the change in the wake structure, which is essential to optimize the air-foil
design. The present work aims to fill in some of these gaps by focusing on the
dynamics (and instabilities) of closely-spaced helical vortices. A recent
experimental investigation of two closely spaced helical vortices generated by
single-bladed rotor, Schröder et al. (2020, 2021) displays a rapid increase of
the core radii and subsequent merging related to the development of a
centrifugal instability (Bayly, 1988). In this case, the centrifugal
instability is triggered by patches of opposite signed vorticity formed by a
protruded fin during the roll-up process. A theoretical analysis by Castillo-
Castellanos et al. (2021) presented the wake geometry produced by a tip-
splitting rotor for all wind turbine and helicopter flight regimes using a
filament approach. The resulting wake deviates from a helical pattern in a
favour of an epicycloidal pattern produced by two interlaced helical vortices
inscribed on top of a larger helical curve. From afar, the vortex structure is
reminiscent of a helical vortex, while simultaneously resembling a vortex pair
aligned with the locally tangent flow. Given these similarities, we expect the
solutions to display features from both systems. In particular, similar
instability mechanisms.

Studies on the stability of helical vortices started a century ago (see for
instance Levy & Forsdyke, 1928) although experimental evidences of short-wave
and long wave instabilities are quite recent (Leweke et al., 2014).
Theoretical predictions for the long-wave instability have been obtained by
Widnall (1972) for a uniform helical vortex. Her work was extended by Gupta &
Loewy (1974) for several helices and by Fukumoto & Miyazaki (1991) to account
the effect of an axial flow within the vortices. As demonstrated by Quaranta
et al. (2015, 2019), the long-wave instability is a local pairing instability
(Lamb, 1945). The instability modes are associated with displacements towards
each other of neighbouring loops at specific locations. Quaranta et al. (2019)
have also shown that their growth rate is well-predicted by considering 3D
perturbations on straight vortices (Robinson & Saffman, 1982). For several
helices, additional theoretical results have been obtained for the global
pairing mode which preserves the helical symmetry of the flow (Okulov, 2004;
Okulov & Sørensen, 2007) using Hardin’s expressions for the induced velocity
of a helical vortex (Kawada, 1936; Hardin, 1982). This mode has been further
analysed in the non-linear regime by Selçuk et al. (2017). Finally, recent
numerical works have also evidenced the presence of the long-wave instability
in the context of rotor-generated vortices (Bhagwat & Leishman, 2000; Walther
et al., 2007; Ivanell et al., 2010; Durán Venegas et al., 2021). Helical
vortices are also unstable to short-wavelength instabilities, due to the
modification of the core structure by curvature, torsion and strain (Kerswell,
2002; Blanco-Rodríguez & Le Dizès, 2016, 2017). This requires detailed
information of the inner core structure, which can be obtained through matched
asymptotic techniques (Blanco-Rodríguez et al., 2015), or through direct
numerical simulations (Selçuk et al., 2017; Brynjell-Rahkola & Henningson,
2020). Pairs display a variety of complex behaviors (Leweke et al., 2016). For
instance, consider a vortex pair of circulations $\Gamma_{1}$ and $\Gamma_{2}$
and core radii $a$, separated by a distance $d$. If $d$ is large enough, the
system is expected to translate with constant speed for
$\Gamma_{1}=-\Gamma_{2}$, or rotate around the vorticity barycentre with
constant rotation rate for all other cases, while $a$ is expected to grow
following a viscous law. Pairs are expected to merge as a critical core size
is eventually reached (Meunier et al., 2002; Josserand & Rossi, 2007). This
behaviour is modified by the development of short-wave instabilities (Meunier
& Leweke, 2005) and long-wave instabilities (Crow, 1970; Jimenez, 1975). This
picture becomes increasingly complex as we consider the interactions between
multiple vortex pairs, like the four vortex system considered by (Crouch,
1997; Fabre & Jacquin, 2000; Fabre et al., 2002).

For this work, we focus on the long-wavelength stability of two closely-spaced
helical vortices using a filament approach. This approach has the advantage of
filtering the short-wavelength perturbations. Additionally, we are interested
in the regime observed in the far-field, since perturbations are expected to
be quickly advected away from the rotor (Durán Venegas et al., 2021). We shall
use as base flow the steady solutions that have been obtained in Castillo-
Castellanos et al. (2021). We shall see that closely-spaced helical vortices
display two kinds of local pairing modes, reminiscent to those observed in
helical vortices and vortex pairs.

The paper is organized as follows. In section 2 we present the framework of
the vortex method applied to the pair of helical vortices. We describe the
numerical procedures to obtain the base flow and analyse its stability. In
section 3, we apply our numerical approach to uniform helices for validation.
We then consider two different configurations. The stability properties for a
pair of helical vortices without a central vortex hub are presented in section
4, while the stability properties for two pairs of helical vortices with a
central vortex hub are presented in section 5. In both cases, the structure of
the unstable modes is presented in detail and we explore the influence of the
main geometric parameters. Additional information regarding the spatio-
temporal development of the instability is provided in the appendix A, where
the linear impulse response for the first configuration is studied using the
approach developed by Durán Venegas et al. (2021). Finally, we present our
main conclusions in section 6.

##  2 Methodology

###  2.1 Finding stationary solutions using a filament approach

####  2.1.1 Filament framework

We use the same vortex filament approach as in Durán Venegas & Le Dizès
(2019). Vorticity is considered to be concentrated along thin vortex filaments
moving as material lines in the fluid according to

$\frac{d\boldsymbol{X}_{j}}{dt}=\boldsymbol{U}(\boldsymbol{X}_{j})-{U}_{z}^{\infty}\boldsymbol{\hat{e}_{z}}$
(1)

where $\boldsymbol{X}_{j}=(r_{j},\theta_{j},z_{j})$ is the position of the
$j$-th vortex filament, ${U}_{z}^{\infty}$ an external velocity field, and
$\boldsymbol{U}(\boldsymbol{X}_{j})=(U_{r},r_{j}\Omega,U_{z})$ is the velocity
induced by the vortices at $\boldsymbol{X}_{j}$. In a frame with rotation rate
$\Omega_{R}$, it is convenient to parametrize the position of vortex filaments
in terms of two angular coordinates, $\zeta=\Omega_{R}(t-t_{0})$ and
$\psi=\Omega_{R}t$, and transform (1) as

$\frac{\partial\boldsymbol{X}_{j}}{\partial\psi}+\frac{\partial\boldsymbol{X}_{j}}{\partial\zeta}=\frac{1}{\Omega_{R}}\left[\boldsymbol{U}(\boldsymbol{X}_{j})+\boldsymbol{U}^{\infty}(\boldsymbol{X}_{j})\right]$
(2)

where
$\boldsymbol{U}^{\infty}(\boldsymbol{X}_{j})=(0,-r_{j}\Omega_{R},-{U}_{z}^{\infty})$
contains the contributions from the rotating frame and the axial velocity. In
this context, $\zeta$ is often referred to as the wake age and corresponds to
a spatial coordinate, while $\psi$ is a proxy of time (Leishman et al., 2002).
Filaments are discretized in straight segments
$[\boldsymbol{X}^{n}_{j},\boldsymbol{X}^{n+1}_{j}]$ in order to compute the
velocity field using the Biot-Savart law. The equation is desingularized using
a cut-off approach with a Gaussian vorticity profile. Local contributions to
the velocity field at $\boldsymbol{X}_{j}$ are obtained by replacing the
straight segments by an arc of circle passing through
$[\boldsymbol{X}^{n-1}_{j},\boldsymbol{X}^{n}_{j},\boldsymbol{X}^{n+1}_{j}]$
and use the cut-off formula (see, Durán Venegas & Le Dizès (2019)).

For this work, we consider basic flow solutions that satisfy,

$\frac{dr_{j}}{d\zeta}=\frac{U_{r}(\boldsymbol{X}_{j})}{\Omega_{R}},\quad\frac{d\theta_{j}}{d\zeta}=\frac{\Omega(\boldsymbol{X}_{j})-\Omega_{R}}{\Omega_{R}},\quad\frac{dz_{j}}{d\zeta}=\frac{U_{z}(\boldsymbol{X}_{j})-U^{\infty}_{z}}{\Omega_{R}}.$
(3)

These solutions are $\psi$-independent solutions to (2). They are stationary
in the rotating frame in the sense that flow remains tangent to the vortex
structure at all times:

$(\boldsymbol{U}(\boldsymbol{X}_{j})+\boldsymbol{U}^{\infty}(\boldsymbol{X}_{j}))\times\boldsymbol{T}_{j}=0$
(4)

where $\boldsymbol{T}_{j}$ is the tangent unit vector. As shown in Castillo-
Castellanos et al. (2021), the spatial evolution of the wake as we move away
from the rotor plane can be obtained by numerically solving (3) with boundary
conditions on the rotor, at $\zeta=0$, where the position of each vortex is
prescribed and far-field boundary conditions at $\zeta\to\infty$. In the
following, we shall focus on the wake geometry in the far-field, where
spatially periodic solutions regime can be obtained.

####  2.1.2 Periodic solutions in the far-field

The solution in the far-field has been obtained in Castillo-Castellanos et al.
(2021). It was shown that it is no longer uniform as for a single helix but it
nevertheless exhibits a certain spatial periodicity. In addition to the
azimuthal symmetry $\theta\to\theta+2\pi/N$ for $N$ vortex pairs, solutions
are invariant by the double operation $z\to z+h$ and
$\theta\to\theta+2\pi/\beta$. The parameters $h$ and $\beta$ are chosen such
that there is a single location in an axial period $h$ where the vortices of a
given pair are at the same azimuth. This azimuth is taken to define the mean
radius $R$ and the separation distance $d$ of each vortex pair. Each vortex
pair tends to form a helical braid on a larger helical structure of radius $R$
and pitch $H$ as illustrated on figure 1.

Figure 1:  (a) Geometric parameters of the periodic solutions: separation
distance $d$, radius $R$, axial pitches $H$ and $h$, core size $a$ and twist
parameter $\beta=H/h$. Three dimensional view for (b) $N=1$ and (c) $N=2$ with
hub vortex for ($R^{*}=7$, $H^{*}/N=5.25$, $\beta=3/4$).

The pitch of this larger helical structure is directly related to $h$ and
$\beta$ by

$H=h\beta.$ (5)

The parameter $\beta$ characterizes the axial periodicity of the solution. It
becomes axially periodic only if $\beta$ is a rational number $p/q$. In that
case, it means that the double helix makes $q$ turns on itself as the large
helix does $p$ turns. The axial period is thus $pH=qh$. In the following, we
only consider rational values of $\beta$.

As soon as the vortex core size $a$ is fixed, the far-field is then defined by
5 non-dimensional parameters which are

$R^{*}\equiv\frac{R}{d},\quad
H^{*}\equiv\frac{H}{d},\quad\varepsilon^{*}\equiv\frac{a}{d},\quad\beta,\quad
N.$ (6)

The solution is obtained by solving the system (3) with the prescribed
symmetry. The frame velocities $(\Omega_{R},U^{\infty}_{z})$ are unknown
quantities. However, these quantities are proportional to the vortex
circulation $\Gamma$ as it is the unique quantity of the vortex system
involving time. For this reason, we can fix $\Gamma$ to 1.

As shown in Castillo-Castellanos et al. (2021), the problem can be treated as
a non-linear minimization problem using an iterative procedure. The
convergence is rapid if we start for each pair from an initial guess given by
an undeformed double helix on a larger helix and estimates for
$(\Omega_{R},U^{\infty}_{z})$ obtained from uniform helices. The converged
solution is found to exhibit spatial variations but in most cases the initial
guess turns out to be a good approximation of the solution (see Castillo-
Castellanos et al., 2021, for details).

In the present study, we consider two different configurations: one composed
of a single helical pair ($N=1$) without central hub (figure 1b), and another
composed of $N=2$ helical pairs with a central hub (figure 1c). Also, we fix
$\varepsilon^{*}=0.1$ and vary the remaining parameters ($R^{*}$, $H^{*}$ and
$\beta$). In the following, solutions that satisfy (3) are denoted
$\boldsymbol{X}_{j}^{B}(\zeta)$.

###  2.2 Inviscid global instability analysis

The stability of $\boldsymbol{X}_{j}^{B}$ is analysed by considering the
evolution of infinitesimal perturbation displacements

$\boldsymbol{X}_{j}(\zeta,\psi)=\boldsymbol{X}_{j}^{B}(\zeta)+\boldsymbol{X}^{\prime}_{j}(\zeta,\psi)$
(7)

Equation (2) is linearized around $\boldsymbol{X}_{j}^{B}(\zeta)$ to obtain a
linear dynamical system

$\displaystyle\frac{\partial r_{i}^{\prime}}{\partial\psi}$ $\displaystyle=$
$\displaystyle-\frac{\partial
r_{i}^{\prime}}{\partial\zeta}+\frac{1}{\Omega_{R}}\sum_{j=1}^{N}\left[\frac{\partial
U_{r}(\vec{X}_{i}^{B})}{\partial r_{j}}r_{j}^{\prime}+\frac{\partial
U_{r}(\vec{X}_{i}^{B})}{\partial\theta_{j}}\theta^{\prime}_{j}+\frac{\partial
U_{r}(\vec{X}_{i}^{B})}{\partial z_{j}}z^{\prime}_{j}\right]$ (8a)
$\displaystyle\frac{\partial\theta_{i}^{\prime}}{\partial\psi}$
$\displaystyle=$
$\displaystyle-\frac{\partial\theta_{i}^{\prime}}{\partial\zeta}+\frac{1}{\Omega_{R}}\sum_{j=1}^{N}\left[\frac{\partial\Omega(\vec{X}_{i}^{B})}{\partial
r_{j}}r_{j}^{\prime}+\frac{\partial\Omega(\vec{X}_{i}^{B})}{\partial\theta_{j}}\theta^{\prime}_{j}+\frac{\partial\Omega(\vec{X}_{i}^{B})}{\partial
z_{j}}z^{\prime}_{j}\right]$ (8b) $\displaystyle\frac{\partial
z_{i}^{\prime}}{\partial\psi}$ $\displaystyle=$ $\displaystyle-\frac{\partial
z_{i}^{\prime}}{\partial\zeta}+\frac{1}{\Omega_{R}}\sum_{j=1}^{N}\left[\frac{\partial
U_{z}(\vec{X}_{i}^{B})}{\partial r_{j}}r_{j}^{\prime}+\frac{\partial
U_{z}(\vec{X}_{i}^{B})}{\partial\theta_{j}}\theta^{\prime}_{j}+\frac{\partial
U_{z}(\vec{X}_{i}^{B})}{\partial z_{j}}z^{\prime}_{j}\right]$ (8c)

which has the following general form

$\frac{\partial\boldsymbol{q}^{\prime}}{\partial\psi}=\mathsfbi{L}(\zeta)\boldsymbol{q}^{\prime}(\zeta,\psi)$
(9)

where $\boldsymbol{q}^{\prime}=(\boldsymbol{X}_{1},\boldsymbol{X}_{2},\cdots)$
is the total displacement vector, and $\mathsfbi{L}$ is a $N$ by $N$ block-
matrix containing the Jacobian terms. The spatial derivative (on $\zeta$) is
evaluated using the same finite differences scheme as the base solution, while
the velocity gradient is derived from the discretized Biot-Savart equations.
Each sub-matrix of $\mathsfbi{L}$ can be written as

$\mathsfbi{L}_{ij}=-\delta_{ij}\frac{\partial}{\partial\zeta}+\frac{1}{\Omega_{R}}\frac{\partial(\boldsymbol{U}(\boldsymbol{X}_{i}^{B}))}{\partial(\boldsymbol{X}_{j})}$
(10)

and contains $3n_{s}$ by $3n_{s}$ elements where $n_{s}$ is the total number
of discretization point of each vortex. The base flow solution is periodic
with respect to $\zeta$ with a period $\zeta_{B}$ that satisfies
$z_{j}(\zeta_{B})=h$ and $\theta_{j}(\zeta_{B})=2\pi/\beta$. The difference
between $\zeta_{B}$ and $\theta_{j}(\zeta_{B})$ comes from the induced angular
velocity $\Omega(\boldsymbol{X}_{j})$ in (3)b. For small values of $\beta$,
the ratio $C_{\theta}\equiv\zeta_{B}/\theta_{j}(\zeta_{B})$ is close to 1, but
decreases as $\beta$ increases (figures 2a and 2b). Each vortex contains 96
points per period $\zeta_{B}$. For the perturbations, we consider a longer
domain with a period $\zeta_{p}=n_{p}\zeta_{B}$ with $n_{p}=32/\beta$ such
that $z_{j}(\zeta_{p})=n_{p}h=32H$ and
$\theta_{j}(\zeta_{p})=n_{p}2\pi/\beta=64\pi$. This gives $n_{p}=24$ and
$n_{p}=128$ for the extreme values of $\beta$ considered ($\beta=3/4$ and
$\beta=4$). Our objective is to apply periodic boundary conditions to the
perturbations. We therefore discretize the operators $\mathsfbi{L}_{ij}$ such
that they satisfy this property.

(a) (b)

Figure 2:  Evolution of the ratio
$C_{\theta}\equiv\zeta_{B}/\theta_{j}(\zeta_{B})$ for ($N=1$, $R^{*}=7$): (a)
$C_{\theta}$ as function $R/H$ for $\beta=3/4$, and (b) as function of $\beta$
for different $R/H$.

The final problem reduces to a linear system with constant coefficients that
admits the formal solution

$\boldsymbol{q}^{\prime}(\zeta,\psi)=\exp(\mathsfbi{L}\psi)\boldsymbol{q}_{0}^{\prime}(\zeta)$
(11)

for any initial perturbation
$\boldsymbol{q}^{\prime}_{0}(\zeta)=\boldsymbol{q}^{\prime}(\zeta,0)$. In
standard fashion, treating equation (11) as an eigenvalue problem can be used
to describe the asymptotic limit $\psi\to\infty$. Here, because the domain is
long but periodic, it also means that we consider perturbations that span the
whole calculation domain.

We may decompose the linear operator as

$\mathsfbi{L}\mathsfbi{\Phi}=\mathsfbi{\Phi}\mathsfbi{\Lambda}$ (12)

where $\mathsfbi{\Phi}$ is a matrix whose columns
$(\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\cdots)$ are the eigenvectors of
$\mathsfbi{L}$, and $\mathsfbi{\Lambda}$ is a matrix whose diagonal entries
are the corresponding (complex) eigenvalues
$(\lambda_{1},\lambda_{2},\cdots)$. Here, $\sigma_{m}\equiv\Real(\lambda_{m})$
characterizes the temporal growth (or decay) of the perturbations, while
$\omega_{m}\equiv\Imag(\lambda_{m})$ characterizes the temporal oscillations.
We may also introduce a dimensionless growth rate $\sigma_{m}^{*}$ and
frequency $\omega_{m}^{*}$, based on the characteristic advection time-scale
of helical pairs $t_{adv}=\Gamma/(H^{2}/N^{2})$. Since $\mathsfbi{L}$ is real-
valued, modes are either real or come in conjugate pairs. Equation (11) can be
written as

$\boldsymbol{q}^{\prime}(\zeta,\psi)=\sum_{m}\boldsymbol{\phi}_{m}(\zeta)\exp{\\{\sigma_{m}\psi+i\omega_{m}\psi\\}}b_{m}+c.c.$
(13)

where the coefficients $b_{m}$ correspond to the coordinates of
$\boldsymbol{q}^{\prime}_{0}(\zeta)=\boldsymbol{q}^{\prime}(\zeta,0)$
expressed in the eigenvector basis, and $c.c.$ indicates the complex
conjugate.

By construction, eigenvectors have the same dimensionality as
$\boldsymbol{q^{\prime}}$,

$\boldsymbol{\phi}_{m}=(\boldsymbol{\phi}_{m}^{(1)},\boldsymbol{\phi}_{m}^{(2)},\cdots)$
(14)

where $\boldsymbol{\phi}_{m}^{(j)}$ is the $m$-th eigenvector of the $j$-th
vortex filament. These eigenvectors represent a displacement vector and can be
expressed in any coordinate system, e.g.,

$\boldsymbol{\phi}_{m}^{(j)}=(\check{r}_{m}^{(j)},\check{\theta}_{m}^{(j)},\check{z}_{m}^{(j)})$
(15)

indicates the components of $\boldsymbol{\phi}_{m}^{(j)}$ in global
cylindrical coordinates. Conversely, displacement perturbations may also be
expressed in terms of the local radial, azimuthal, and axial coordinates
associated with a uniform helix $\mathscr{H}$ of radius $R$ and pitch $H$,

$\boldsymbol{\phi}_{m}^{(j)}=(\check{\rho}_{m}^{(j)},\check{\phi}_{m}^{(j)},\check{s}_{m}^{(j)}).$
(16)

Due to the spatial periodicity, the eigenvectors can be expanded on a discrete
Fourier basis:

$\boldsymbol{\phi}_{m}=\sum_{n=0}^{n_{s}-1}\boldsymbol{\hat{q}}_{mn}\exp{\left\\{\frac{\mathrm{i}2\pi
n\zeta}{\zeta_{p}}\right\\}}$ (17)

where the azimuthal wavenumber $(k=n\Delta k,\Delta k=\beta/n_{p})$ is
normalized to ensure that $k=1$ corresponds to one helix turn, and as
mentioned above $n_{s}$ is the number of points used to discretize each
vortex.

##  3 Validation of the stability analysis: uniform helices

We validate our approach against existing results on uniform helical vortices
of pitch $H$, radius $R$ and core size $a$. A notable difference with respect
to Widnall (1972) is that we require finding a moving reference frame
compatible with (3). Since the motion of helical vortices is described by
constant rotation rate $\Omega$ and axial velocity $U_{z}$, these vortices are
unperturbed by an additional rotation of angular velocity $\Omega^{a}$ and
translation of axial velocity $U_{z}^{a}$ provided that

$\frac{U_{z}^{a}}{\Omega^{a}}=\pm\frac{H}{2\pi}$ (18)

is satisfied, where the sign is positive for right-handed helices and vice
versa. For this comparison, we consider a rotating frame of reference

$\Omega_{R}=\Omega\mp U_{z}\frac{2\pi}{H},\quad\quad U_{z}^{\infty}=0$ (19)

which is commonly used in numerical simulations.

(a) (b) (c)

Figure 3:  Stability of a uniform helix for $R/H=\pi/5$ and $a_{e}/R=0.1$: (a) $\sigma_{m}^{*}$ as function of $\omega_{m}^{*}$ in the rotating frame; (b) Fourier spectra of the mode indicated as a red mark in (a); and (c) $\sigma_{m}^{*}$ as function of $k_{0}$.  | $a_{e}/R=0.10$ | $a_{e}/R=0.20$ | $a_{e}/R=0.33$
---|---|---|---
|  1 Peak |  2 Peak |  3 Peak |  1 Peak |  2 Peak |  3 Peak |  1 Peak |  2 Peak |  3 Peak
$\sigma_{m}^{*}$ | 1.575 | 1.511 | 1.323 | 1.554 | 1.497 | 1.381 | 1.538 | 1.485 | 1.395
$\omega_{m}^{*}$ | 3.405 | 10.202 | 16.980 | 3.186 | 9.639 | 15.985 | 3.028 | 9.160 | 15.190
$k_{0}$ | -0.500 | -1.500 | -2.500 | -0.500 | -1.516 | -2.516 | -0.500 | -1.516 | -2.516
$c_{0}$ | 1.006 | 1.004 | 1.003 | 1.005 | 1.003 | 1.001 | 1.004 | 1.002 | 1.001
Table 1:  Growth rate, frequency, wavenumber and phase velocity for the most
unstable modes of a uniform helix for $R/H=\pi/5$ and different core sizes.
Here, $\Delta k=0.01625$.

(a) (b)

Figure 4:  Stability of a uniform helix for $R/H=\pi/5$ and different core
sizes: (a) $\sigma_{m}^{*}$ and (b) $\omega_{m}^{*}$ as function of $k_{0}$,
where $\omega_{m}^{*}$ is evaluated in the advection frame.

The frequency spectra for a helix of pitch $H/R=\pi/5$ and $a_{e}/R=0.1$ is
presented in figure 3a. By construction, the spectrum is symmetric with
respect to zero and only the positive frequencies are displayed. For all three
values, the most unstable frequencies are located near $\omega_{m}^{*}=3$ with
additional local maxima near odd multiples of this frequency (see left panel
in table 1). Figure 3b displays the Fourier spectrum of $\check{z}_{m}$ for
the most unstable mode and $a_{e}/R=0.33$, characterized by a single peak at
$k=-0.500$. We can show that each eigenvector has the form of a complex wave
with a single dominant wavenumber $k=k_{0}$, such that a direct correspondence
between $\sigma_{m}$, $\omega_{m}$ and $k$ can be established, see figure 3c.

Perturbations propagate along the structure as the sum of travelling waves

$z^{\prime}(\zeta,\psi)=\sum_{m}\hat{z}_{m}\cos(\omega_{m}\psi+k_{0}\zeta/C_{\theta})\exp{\left\\{\sigma_{m}\psi\right\\}}$
(20)

with phase velocity

$c_{0}=-\omega_{m}C_{\theta}/k_{0}.$ (21)

A similar behaviour is observed for radial and azimuthal perturbations. As
noted by Brynjell-Rahkola & Henningson (2020), the frequencies obtained in a
frame rotating with $\Omega_{R}$ can be mapped into a second reference frame
rotating with $\Omega_{R}^{\prime}$ through

$\omega_{m}^{\prime}=\omega_{m}+(\Omega_{R}^{\prime}-\Omega_{R})k/C_{\theta}.$
(22)

For instance, in our example most of the tangential velocity comes from the
moving frame, i.e., $|\Omega|\ll|\Omega-\Omega_{R}|$, such that perturbations
are advected with $c_{0}$ close to 1, see table 1. If, instead, we consider a
frame moving with the vortex elements, perturbations are advected with $c_{0}$
close to 0 for the same wavenumber (figure 4). These results are in good
agreement with the stability curves presented in figure 5d by Widnall (1972)
and figure 3 by Quaranta et al. (2015) for the same parameters.

Our approach manages to compute the stability properties of uniform helices
without making any particular assumptions regarding the spatial structure of
the base flow and perturbations. This will be useful for studying the helical
pairs since their geometry makes difficult deducing such spatial structures.

##  4 Stability of one vortex pair without central hub vortex

###  4.1 Typical displacement modes for leapfrogging wakes

For each pair, we may introduce the following decomposition

$\boldsymbol{X}_{+}\equiv(\boldsymbol{X}_{1}+\boldsymbol{X}_{2})/2,\quad\quad\boldsymbol{X}_{-}\equiv(\boldsymbol{X}_{1}-\boldsymbol{X}_{2})/2$
(23)

where $\boldsymbol{X}_{+}$ characterises the large-scale pattern traced by the
vorticity barycentre, while $\boldsymbol{X}_{-}$ represents the rotation of
the vortex pair relative to $\boldsymbol{X}_{+}$. Figure 1b depicts a
leapfrogging wake with $\beta=p/q=0.75$, where $\boldsymbol{X}_{+}$ is
represented as a tube enclosing the vortex pair. The vortex pair completes a
rotation every $p=3$ turns, while the structure repeats itself every $q=4$
turns. Note that $\boldsymbol{X}_{+}$ and $\mathscr{H}$ are close but not
exactly equal due to the effect of self-induction. In a similar vein, we may
introduce the following decomposition

$\boldsymbol{\phi}_{m}^{+}\equiv(\boldsymbol{\phi}_{m}^{(1)}+\boldsymbol{\phi}_{m}^{(2)})/2,\quad\quad\boldsymbol{\phi}_{m}^{-}\equiv(\boldsymbol{\phi}_{m}^{(1)}-\boldsymbol{\phi}_{m}^{(2)})/2$
(24)

where $\boldsymbol{\phi}_{m}^{+}$ and $\boldsymbol{\phi}_{m}^{-}$ indicate the
displacement modes of $\boldsymbol{X}_{+}$ and $\boldsymbol{X}_{-}$,
respectively.

Figure 5a presents the frequency spectrum, which is characterized by a set of
modes distributed over three contiguous lobes at low frequencies, and a second
set of modes over an additional lobe at higher frequencies. The maximum growth
rate is observed near $\omega_{m}=4.8$, with additional local maxima (in
descending order) near $\omega_{m}=51.9$, $\omega_{m}=14.5$, and
$\omega_{m}=24.0$. Dimensionless growth rates are larger than the equivalent
helical vortex. This can be explained by the change in the effective distance
separating neighbouring loops due to rotation of the pair. Low frequency modes
are clearly reminiscent of the unstable modes for the equivalent uniform
helices (figure 5b). Figure 6 shows the eigenvectors $\check{s}_{m}^{+}$ and
$\check{s}_{m}^{-}$ corresponding to the most unstable modes along the local
axial direction. Based on their structure, displacements modes can be
classified as being symmetric ($S_{1}$ to $S_{3}$) or anti-symmetric ($A_{1}$)
with respect to the large helical pattern. For symmetric modes
$\check{s}_{m}^{+}$ has a nearly constant envelope and $\check{s}_{m}^{-}$ has
a sinusoidal envelope, whereas anti-symmetric modes display the opposite
behaviour.

(a) (b)

Figure 5:  Frequency spectra for (a) a leapfrogging wake ($R^{*}=7$,
$H^{*}=5.25$, $\beta=3/4$); and (b) a uniform helix ($H/R=0.75$,
$a_{e}/R=0.1$) in the rotating frame.

(a) (b)

(c) (d)

Figure 6:  Complex eigenvectors for the case in figure 5a. Perturbations in
the local axial direction (a,c) $\check{s}_{m}^{+}$, and (b,d)
$\check{s}_{m}^{-}$ shown for modes (a,b) $S_{1}$, and (c,d) $A_{1}$.

(a) (b) (c)

Figure 7:  Fourier spectrum of modes: (a) $S_{1}$, (b) $S_{2}$, and (c) $A_{1}$, for the case in figure 5a.  Mode | Frequency | Wavenumbers | Phase velocity
---|---|---|---
| $\omega_{m}^{*}$ | $k_{0}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $c_{0}$
$S_{1}$ | $4.824$ | $-0.500$ | $(-1.250,+0.250)$ | $(-2.000,+1.000)$ | $(-2.750,+1.750)$ | $1.002$
$S_{2}$ | $14.451$ | $-1.500$ | $(-2.250,-0.750)$ | $(-3.000,+0.000)$ | $(-3.750,+0.750)$ | $1.000$
$A_{1}$ | $51.849$ | $-5.406$ | $(-6.156,-4.656)$ | $(-6.906,-3.906)$ | $(-7.656,-3.156)$ | $0.996$
Table 2: Leading wavenumbers of the modes shown in figure 7. Here, $\Delta
k=0.03125$.

Each mode seems to display a dominant wavenumber $k_{0}$ with additional peaks
(in descending order of magnitude) at wavenumbers $k_{n}\approx k_{0}\pm
n\beta$ for integer values of $n$ (see Fourier spectra in figures 7a-c and
table 2). This pattern is also observed in the spectra of $\check{z}_{m}$ and
$\check{r}_{m}$. From these observations, we infer that perturbations
propagate along the structure as

$s^{\prime}_{j}(\zeta,\psi)\approx\sum_{m}\left[\bigg{[}\hat{s}_{m0}+\sum_{n=1}\hat{s}_{mn}\cos{(n\beta\zeta/C_{\theta}+\hat{\varphi}_{mn})}\bigg{]}\cos(\omega_{m}\psi+k_{0}\zeta/C_{\theta})\exp{\left\\{\sigma_{m}\psi\right\\}}\right]$
(25)

where $\hat{s}_{mn}$ and $\hat{\varphi}_{mn}$ are the amplitudes and phase
differences measured with respect to $k_{0}$. If we consider the leading terms
in (25), for symmetric modes $s_{+}^{\prime}$ roughly corresponds to a
travelling wave, whereas $s_{-}^{\prime}$ approximates a wave that is
modulated in amplitude by a sine function of period $\zeta_{B}$, i.e., the
periodicity of the base flow (figures 7a-b). Anti-symmetric modes display the
opposite behaviour: $s_{-}^{\prime}$ approximates a travelling wave, whereas
$s_{+}^{\prime}$ is modulated in amplitude by a function of period
$\zeta_{B}$. In both cases, the contributions from higher order terms also
correspond to waves modulated in amplitude by multiples of $\zeta_{B}$.

Figure 8:  Base state from figure 1b perturbed by modes (a,b) $S_{1}$ and
(c,d) $S_{2}$. (a,c) Three dimensional and (b,d) developed plan views
illustrating the local pairing modes.  Figure 9:  Base state from figure 1b
perturbed by mode $A_{1}$. (a) Three dimensional and (b) developed plan view.
Arrows indicate the displacements at the center plane in fig. (a).

(a) (b)

Figure 10:  Function $(z_{i}^{\prime}\star z_{j}^{\prime})$ for the case in
figure 1b and modes: (a) $S_{2}$, and (b) $A_{1}$.

Figure 8a presents the deformation due to the symmetric mode $S_{1}$ by
plotting the perturbed geometry
$\boldsymbol{X}_{j}^{B}+\boldsymbol{X}_{j}^{\prime}$. A developed plan view
illustrates $p=1$ localized pairing events every $q=2$ neighbouring turns of
the large-scale pattern (seen as dashed lines in figures 8a-b). An additional
example corresponding to mode $S_{2}$ is shown in figures 8c-d, where $p=3$
localized pairing events are observed every $q=2$ neighbouring turns. This
behaviour was expected since displacements of the pair are predominantly
aligned, the large-scale pattern behaves much like a uniform helix where
perturbations with wavenumber $k_{0}=p/q$ repeat after $p$ cycles and display
local pairing events at $q$ azimuthal locations (Widnall (1972); Quaranta et
al. (2015)).

Anti-symmetric modes behave in a different manner. Here, the two vortices move
towards (or away from) one another such that $\boldsymbol{X}_{+}$ deforms much
less and only at specific positions (figure 9a). In other words, the pair
predominantly displays an anti-symmetric motion with respect to the helical
structure, hence the name. Displacements are localized and not necessarily
aligned with the rotation of the pair. For instance, at the azimuths where
displacements are perpendicular to the line connecting the two vortices the
structure is twisted back and forth, whereas when displacements are parallel,
the separation distance expands and contracts which could potentially trigger
the merging of the vortex pair. This localization can be deduced from the
envelopes of the corresponding eigenvectors, illustrated as dashed lines
enclosing each vortex in figure 9b. Despite the fact that displacements mainly
concern the small-scale structure, anti-symmetric modes still correspond to a
form of large-scale pairing. If we consider a longitudinal cut, the
displacements of a given vortex are in phase opposition with respect to one of
the vortices in the neighbouring turn but not with its companion (see arrows
at the bottom part of figure 9a). A more quantitative way to illustrate this
pairing is through the spatial cross-correlation

$(z_{i}^{\prime}\star
z_{j}^{\prime})(\Delta\theta)\equiv\int_{\infty}^{\infty}z_{i}^{\prime}(\zeta/C_{\theta})z_{j}^{\prime}(\zeta/C_{\theta}+\Delta\theta)d\zeta$
(26)

where $\Delta\theta$ is the delay in angular position. For instance, for
$S_{2}$ the displacements between neighbouring turns, i.e.
$\Delta\theta=2\pi$, are well anti-correlated since perturbations are in
opposition of phase (figure 10a). Conversely, for mode $A_{1}$ the auto-
correlation between consecutive turns is negative, while the cross-correlation
is positive (figure 10b), meaning that vortices move in opposite directions,
but one of them is paired with one of the vortices in the neighbouring loop.
However, due to the geometry, the choice of the characteristic wavenumber
$k_{0}$ is not obvious. For instance, increasing $H^{*}$ from $5.25$ to $10.5$
shifts $k_{0}=5.406$ to $k_{0}=5.938$ for mode $A_{1}$, while increasing
$\beta$ from $3/4$ to $7/2$ shifts the wavenumber to $k_{0}=5.552$. However,
perturbations between consecutive pairs remain anti-correlated, suggesting the
dominant wavenumber is selected by the geometry as the one that amplifies the
local pairing.

In the appendix A, we also analyse the long-wave instability of our solutions
using the linear impulse response approach developed in Durán Venegas et al.
(2021). The idea is to solve the linear perturbation equation (11) using a
Dirac as initial condition and analyse the spatio-temporal growth of the
resulting wave packet. A sufficiently long domain is considered so that the
wavepacket does not reach the boundaries during the length of the simulations.
As shown in the appendix, the temporal modes can be recovered but their study
is more difficult to perform with the linear impulse response as all the
instability modes are simultaneously excited. However, the linear impulse
response provides additional information by telling us how the instability
spreads in space. In particular, we are able to show that although the most
unstable anti-symmetric perturbations propagate at a similar speed than the
symmetric perturbations, their spreading in space is much less important. The
transition from convective to absolute instability is therefore expected to be
associated with the symmetric perturbations. Moreover, as for the linear
spectrum, we also demonstrate that the spatio-temporal evolution of these
symmetric perturbations can be well-described by the spatio-temporal evolution
of the linear impulse response on a uniform helical vortex of large core size.

###  4.2 Geometric dependency of the most unstable modes

(a) (b) (c)

(d) (e) (f)

Figure 11:  Stability curves, growth rates and dominant wavenumbers for
$\beta=3/4$: (a,b,c) as function of $H^{*}/R^{*}$ for $R^{*}=7$ ; (d,e,f) as
function of $R^{*}$ for $H^{*}/R^{*}=0.75$.

(a) (b) (c)

Figure 12:  Growth rates, dominant wavenumbers, and phase velocity as a
function of $\beta$ for $(R^{*}=7,H^{*}=5.25)$.

Given the pairing occurs between neighbouring turns, the growth rates are
expected to display a strong dependency on the relative pitch $H^{*}/R^{*}$.
In general, the maximum growth rate is larger than the maximum growth rate
obtained for uniform helices with the total circulation of the vortex pair
(figure 11a), which is also larger than $\sigma=\pi/2(\Gamma/H^{2})$ obtained
for an array of point vortices of circulation $\Gamma$ and separated by $H$.
In this limit, the symmetric and anti-symmetric modes have comparable growth
rates (figure 11a-b). For symmetric modes, the dimensionless growth rate
initially increases for mode $S_{1}$ before vanishing, while the modes
$S_{2}$, $S_{3}$, and so on, gradually vanish as we increase $H^{*}/R^{*}$.
For anti-symmetric modes, growth rates initially increase at a faster rate,
but also vanish more quickly, while the dominant wavenumber is approximatively
given by $k_{0}=0.45H^{*}/R^{*}$ for $H^{*}/R^{*}>2$, see figure 11c.

It is also interesting to fix the relative pitch $H^{*}/R^{*}$ and change the
separation distance through $R^{*}$ (figure 11d). On the one hand, for modes
$S_{1}$ and $S_{2}$ the stability curves are essentially unchanged, while the
growth rate of $S_{3}$ decreases as the separation distance becomes smaller
(figure 11c). This is also reminiscent of uniform helical vortices, where the
maximum growth rates also decrease as the effective core size becomes smaller
(Widnall, 1972; Quaranta et al., 2015). Symmetric and anti-symmetric growth
rate curves overlap for $R^{*}<6$. For mode $A_{1}$, the growth rate remains
constant (figure 11e), while $k_{0}$ seems to evolve linearly with $k_{0}\sim
0.89R^{*}$ (figure 11f). This suggests that mode $A_{1}$ tends to select a
particular wavelength along the large scale structure instead of an azimuthal
wavenumber $k$. We shall revisit this matter in the following section.

Now consider the influence of the parameter $\beta$. Figure 12a presents the
stability curves as we move from leapfrogging to sparsely braided wakes. As
seen in figure 12a, growth rates fluctuate around a constant value, which can
be explained since the effective separation distance between neighbouring
vortices is modified only by a factor $\cos(\phi)$ due to relative orientation
of the vortex pairs, where $\phi=2\pi/q$ for $\beta=p/q$. As we increase
$\beta$, anti-symmetric modes are observed over an wider range of wavenumbers,
but $k_{0}$ varies little for the most unstable mode (figure 12b). From this
point of view, the most noticeable change concerns the phase velocity $c_{0}$
which decreases for large $\beta$ (figure 11c). For leapfrogging wakes $c_{0}$
is close to 1 but decreases as the structure becomes braided, indicating that
vortex elements move with an angular velocity smaller than $\Omega_{R}$.

##  5 Stability of two vortex pair with central hub vortex

###  5.1 Typical displacement modes for leapfrogging wakes

Figure 13:  Frequency spectra for (a) two vortex pairs and a hub vortex
($R^{*}=7$, $H^{*}=10.5$, $\beta=3/4$); and (b) two interlaces helices
($H/R=1.5$, $a_{e}/R=0.1$) in the rotating frame. In figs. (a) and (b), colour
indicates the phase difference between vortex pairs and helices, respectively.

(a) (b)

Figure 14:  (a) Three dimensional view of the base state from figure 1c perturbed by mode $S_{0}$, and (b) developed plan view illustrating the global pairing.  Mode | Frequency | Wavenumbers | Phase velocity
---|---|---|---
| $\omega_{m}^{*}$ | $k_{0}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $c_{0}$
$S_{0}$ | $0.015$ | $0.000$ | $(-0.781,+0.719)$ | $(-1.531,+1.469)$ | $(-2.281,+2.219)$ | $0.952$
$S_{1}$ | $9.863$ | $-1.031$ | $(-1.781,+0.281)$ | $(-2.531,+0.469)$ | $(-3.281,+1.219)$ | $0.987$
$S_{2}$ | $19.086$ | $-2.000$ | $(-2.750,-1.250)$ | $(-3.500,-0.500)$ | $(-4.250,+0.250)$ | $0.984$
$A_{1}$ | $51.846$ | $-5.468$ | $(-6.219,-4.719)$ | $(-6.969,-3.969)$ | $(-7.719,-1.719)$ | $0.978$
Table 3: Leading wavenumbers of the modes shown in figure 13. Here, $\Delta
k=0.03125$.

A similar stability analysis can be performed for the case of two interlaced
vortex pairs with a central hub vortex, illustrated in figure 1c. We use
similar geometric parameters, but double the pitch as to preserve the mean
axial distance between neighbouring turns. The results are summarized in
figures 13, 14 and 15. The frequency spectra presented in figure 13a is
characterized by a first set of modes distributed over three overlapping lobes
at low frequencies, and a second set distributed over two additional lobes at
higher frequencies. As in the previous case, low frequency modes are symmetric
modes, while the higher frequency modes are anti-symmetric. The structure of
the eigenvectors is the same as before with dominant wavenumbers at $k_{0}$
and additional harmonic terms at $k_{n}=k_{0}\pm n\beta$, see table 3.

As expected, symmetric modes are reminiscent of the unstable modes obtained
for two equivalent helical vortices (figure 13b). Displacements between
neighbouring turns are out-of-phase for modes in the branch containing $S_{0}$
and $S_{2}$ (in blue), and phase-aligned for modes in the branches containing
$S_{1}$ and $S_{3}$ (in orange). These displacements result in a local pairing
at $2k$ azimuthal positions per turn of the large-scale helix. Mode $S_{0}$
corresponds to a special case. Deformations between neighbouring turns are
out-of-phase, where one vortex pair expands in the radial direction, while the
other one contracts, resulting in a uniform pairing of the large-scale pattern
along the azimuthal direction (figures 14a-b). This is analogous to the global
pairing mode $k=0$ of the interlaced helices (see, for instance Okulov &
Sørensen (2009); Quaranta et al. (2019)).

Finally, we note that anti-symmetric modes are observed over a similar range
of values as in the case of one vortex pair. One notable difference is that
two modes are now obtained for the same $k_{0}$ (each using a different colour
in figure 13a), where each pair has a similar structure but shifted in phase
(not shown).

###  5.2 Geometric dependency of the most unstable modes

(a) (b) (c)

(d) (e) (f)

Figure 15:  Stability curves, growth rates and dominant wavenumbers as
function of $R^{*}$ for $(H^{*}/R^{*}=1.5,\beta=3/4)$ (a,b,c), and as function
of $H^{*}/R^{*}$ for $(R^{*}=7,\beta=3/4)$. Fig. (e) uses dimensionless growth
rate $\sigma^{\dagger}\equiv\sigma\Gamma/(4\pi R^{2})$ in the limit of large
$H^{*}/R^{*}$, see text.

As expected, symmetric and anti-symmetric modes display a different behaviour
as we vary the geometric parameters. Varying the separation distance has a
small influence on the symmetric modes. Similarly to the case of one vortex
pair, the change in growth rates is reminiscent to that of varying the
effective core size in uniform helices: small for $S_{0}$ and $S_{1}$, but
more important for the higher wavenumbers (figure 15a). For anti-symmetric
modes, the change in growth rates is also small (figure 15b), while the
dominant wavenumber $k_{0}$ increases linearly with $R^{*}$ (figure 15c). For
$H^{*}/R^{*}=1.5$, $k_{0}$ was found to be roughly the same as in the case of
one vortex pair with equal effective pitch, suggesting the modes are selected
by the same pairing mechanism described in §4.

Varying the relative pitch shows the transition between two different regimes.
For small $H^{*}/R^{*}$, stability curves have a similar structure as before:
symmetric modes distributed over two or more branches with small $k_{0}$ and
anti-symmetric modes distributed over two overlapping branches with larger
$k_{0}$. Growth rates are larger than predicted rates for the equivalent
uniform helices, but still close to the point vortex prediction
$\sigma\sim\Gamma N^{2}/H^{2}$. For large $H^{*}/R^{*}$, the stability curves
are characterized by two branches at small $k_{0}$ containing modes $S_{0}$
and $S_{1}$, and single branch containing $A_{1}$ at larger $k_{0}$. Modes
$S_{1}$ and $A_{1}$ are shifted towards larger wavenumbers. Instead of
vanishing, the dimensionless growth rate $\sigma^{*}$ continues to increase,
pointing to a different scaling law with $\sigma\sim\Gamma/R^{2}$ for large
$H^{*}/R^{*}$ (figures 15d-e). A similar behaviour is observed for the most
unstable wavenumber $k_{0}$. For small $H^{*}/R^{*}$, it seems to scale with
the effective pitch $H^{*}/NR^{*}$, whereas for large $H^{*}/R^{*}$ it follows
a different scaling law in $H^{*}/R^{*}$ (figure 15f). This change of regime
can be understood as follows. For the case a single pair ($N=1$), the limit of
large $H^{*}/R^{*}$ leads to a pair of parallel co-rotating vortices which are
known to be stable with respect to long wavelength perturbations (Jimenez,
1975). For the case of $N=2$ vortex pairs with a central hub, the same limit
provides a system composed of two co-rotating pairs of vortices of circulation
$\Gamma$ and one counter-rotating vortex of circulation $-4\Gamma$ at the
centre. For this configuration, the instability is necessarily controlled by
the distance between the vorticity centroids ($b=2R$ in our current notation),
and the separation distance $d$. This configuration is similar to four-vortex
systems (Crouch, 1997; Fabre et al., 2002), which are known to be unstable,
with a maximum growth rate scaling with $\Gamma/b^{2}$, and a most unstable
wavenumber varying with $b$ and $d$.

##  6 Discussion

In this article, we have studied the long-wave stability properties of
closely-spaced helical vortex pairs using a cut-off filament approach. The
considered base flow configuration corresponds to the far-wake produced by a
rotor emitting two distinct vortices near the tips of each blade and which was
studied previously in (Castillo-Castellanos et al., 2021). Both the temporal
linear spectrum and the linear impulse response have been analysed, but the
linear spectrum has been found to be much more convenient to identify the
different instability modes. We have classified these modes in two main
groups. Symmetric modes are characterized by a block displacement of the
vortex pair, analogous to the local pairing modes in helical geometries. Anti-
symmetric modes are characterized by a mirrored displacement with respect to
the vorticity centre of the pair, with the most unstable mode corresponding to
a local pairing between one member of a pair with the other member of a pair
in the neighbouring turn. These modes are particularly important since they
could trigger the merging of vortex pairs. We have also considered the
dependency of the stability properties with respect to the relative pitch, the
separation distance, and the twist parameter. We have identified the regions
in the parameter space where each mode is dominant. Our observations also
suggest that the pairing mechanism associated with the anti-symmetric modes
amplifies a specific axial wavelength (in the developed plane) instead of an
azimuthal wavelength, reminiscent of the anti-symmetric modes observed in
four-vortex systems. A similar pairing mechanism has also been observed for
the case of two pairs of helical vortices with one central hub. However, in
this case, the instability does not disappear in the limit of large pitch and
exhibits a maximum growth rate scaling with $R/d$.

Our analysis assumes the wavelength of the displacement perturbations to be
large with respect to the core size. For helical vortices, Quaranta et al.
(2015) estimated the limit of validity in the form of wavenumber $k\leq
k_{l}$. For the values used in figure 9, this upper limit corresponds to
$k_{l}\approx 0.40(a_{e}/R)^{-1}$ and $k_{l}\approx 1.60(a_{e}/d)^{-1}$. Since
we consider relatively thin vortices ($a_{e}/R\sim 0.01$ and $a_{e}/d\sim
0.1$), the limit of validity of the long-wave approximation should not be a
concern. Unlike Quaranta et al. (2015), which uses analytical expressions, our
approach considers the filaments as a sequence of straight segments to compute
the Jacobian matrix using semi-analytical expressions. This approach has been
validated using known results for the long-wave instability of uniform
helices. However, it is more general as it does not require prior knowledge of
the spatial structure of the instability modes. It also provides the complete
spectrum and applies to any stationary vortex solution.

By using a filament approach, we have neglected what is occurring in the
vortex cores. Yet, vortex cores are expected to be distorted by curvature and
straining effects (Blanco-Rodríguez et al., 2015). Moreover, these
deformations are also responsible for the short-wave instabilities developing
in vortex cores. Depending on the geometric parameters, these short-wave
instabilities can become dominant. For instance, the elliptical instability is
expected to grow with $\Gamma/d^{2}$ instead of $\Gamma/H^{2}$ although with
different pre-factors (Roy et al., 2008; Blanco-Rodríguez & Le Dizès, 2016),
while the curvature instability (Blanco-Rodríguez & Le Dizès, 2017) is also
expected to be present and important if the vortex core exhibits an axial jet.
None of these have been considered here. However, Brynjell-Rahkola &
Henningson (2020) have shown that it is indeed possible to analyse the
stability of uniform helical vortices with respect to both short and large
wavelengths using a DNS from an initial condition obtained by the filament
solutions together with a prescribed vortex model in the cores. It would be
interested to implement such an approach to our configurations in order to
analyse the competition between both instabilities.

## Acknowledgments

This work is part of the French-German project TWIN-HELIX, supported by the
Agence Nationale de la Recherche (grant no. ANR-17-CE06-0018) and the Deutsche
Forschungsgemeinschaft (grant no. 391677260). The authors are grateful to
Eduardo Durán Venegas and Thomas Leweke for their valuable contributions.

##  Appendix A Space-time impulse response

The space-time impulse response of $\boldsymbol{X}_{j}^{B}$ is studied as in
Durán Venegas et al. (2021). We introduce a Dirac impulse as to excite all the
wavenumber components with equal amplitude and follow the evolution of the
resulting wavepacket. We consider two types of initial perturbation as to
preferentially excite the symmetric and anti-symmetric modes

$\boldsymbol{q}_{0}^{\prime}(\zeta)=\begin{cases}(0,0,A_{0}\delta(\zeta-\zeta_{p}),0,0,~{}~{}A_{0}\delta(\zeta-\zeta_{0}))&\mbox{Case
A (in-phase)}\\\
(0,0,A_{0}\delta(\zeta-\zeta_{p}),0,0,-A_{0}\delta(\zeta-\zeta_{0}))&\mbox{Case
B (out-of-phase)}\end{cases}$ (27)

where $A_{0}$ and $\zeta_{0}$ indicate the amplitude and wake coordinate of
the initial perturbation and $\delta$ is the Dirac delta function. Then, we
use (11) recursively to obtain
$\boldsymbol{q}^{\prime}(\zeta,\psi=m\Delta\psi)$ for $(m=1,2,\cdots)$ until
the exponential regime is established. In general, the exponential regime is
established quickly and our calculation domain is considered to be long enough
to avoid boundary effects.

Temporal growth rates are estimated from the impulse response as follows. At
each time $\psi$, we apply the spatial Fourier transform along $\zeta$ to the
axial displacement. Then, for each azimuthal wavenumber $k$ we can monitor the
amplitude of each Fourier mode and estimate its growth rate by

$\sigma_{j}(k)\approx\frac{\partial(\log\lVert\hat{z}^{\prime}_{j}(k)\rVert)}{\partial\psi}.$
(28)

It is also interesting to consider the growth rate

$\sigma_{V}(V_{\psi})\approx\frac{\log(|z^{\prime}_{j}(\psi,\zeta_{0}+V_{\psi}\psi)|)}{\psi}$
(29)

in the frame moving along the vortex structure with an angular velocity
$V_{\psi}$. In particular, we are interested in the velocity $V_{\psi}^{\max}$
at which the growth rate is maximum and the upper and lower limits at which
the perturbation grows, $V_{\psi}^{+}$ and $V_{\psi}^{-}$. This provides a
quantitative criterium to identify convectively unstable and absolutely
unstable flows (Huerre & Monkewitz, 1990).

(a) (b)

Figure 16:  Impulse response for ($R^{*}=7$, $H^{*}=5.25$, $\beta=3/4$). (a,b)
$l_{2}$-norm of the perturbations for cases A and B, respectively.

Consider the propagation of the initial perturbation corresponding to case A.
The norm of the total displacement $\boldsymbol{q}^{\prime}$ grows
exponentially with constant rate $\sigma^{*}$ while the axial displacements
$z_{1}^{\prime}$ and $z_{2}^{\prime}$ display variable growth rates with a
period comparable to the characteristic time of the vortex pair (figure 16a).
The symmetric component $z_{+}^{\prime}$ is dominant and displays a constant
growth rate equal to $\sigma^{*}$. A similar behaviour is observed for radial
and angular displacements (not shown). The spatio-temporal evolution
corresponding to case B provides similar growth rates, but $z_{-}^{\prime}$ is
now dominant and the oscillations in $z_{1}^{\prime}$ and $z_{2}^{\prime}$ are
less pronounced (figure 16b). Cases A and B have slightly different growth
rates, illustrating a clear dependency on the initial conditions (figure 16a
and 16b). The observed growth rates correspond to the predicted rates for the
most unstable symmetric and anti-symmetric modes, $S_{1}$ and $A_{1}$,
respectively. As expected, if we observe long enough both cases eventually
arrive to whichever one has the largest growth rate.

(a) (b) (c)

Figure 17:  Fourier amplitudes of (a) $z_{+}^{\prime}$ and (b)
$z_{-}^{\prime}$ for case A, where each line corresponds to a different
instant, and (c) growth rates obtained using (28).

(a) (b) (c)

(d) (e) (f)

Figure 18:  Fourier amplitudes of (a,d) $z_{+}^{\prime}$ and (b,e)
$z_{-}^{\prime}$, where each line corresponds to a different instant, and
(c,d) growth rates obtained using (28). Figs. (a-c) correspond to case A with
a high pass filter, and figs. (d-f) to case B with a low pass filter, see
text.

Figure 17a (resp. 17b) displays the Fourier spectra of $z_{+}^{\prime}$ (resp.
$z_{-}^{\prime}$) taken at regular intervals, where the most unstable
wavenumbers are clearly identified. For each wavenumber $k$, the corresponding
growth rate obtained using (28) is shown in figure 17c. Results are in good
agreement with the eigenvalue problem. For instance, in $z_{+}^{\prime}$ the
most unstable wavenumbers correspond to the dominant wavenumbers of $S_{1}$
and $S_{2}$. Additional peaks observed at $k=1$ and $k=2$ correspond to even
harmonics, while peaks in $z_{-}^{\prime}$ (in orange) observed at $k=0.25$
and $k=1.25$ correspond to odd harmonics of $S_{1}$. Similar results are
observed for mode $A_{1}$ and to a smaller degree for mode $S_{2}$.

Issues in recovering the harmonics of mode $S_{2}$ can be explained by the
overlapping between symmetric and anti-symmetric modes. We may take advantage
of the separation in temporal frequency to reduce this overlapping. For
instance, applying a high pass filter to the wavepackets from case A before
evaluating the Fourier coefficients, i.e., filtering the anti-symmetric modes,
we recover the growth rates of symmetric modes over an extensive range of
wavenumbers (figures 18a-c). Conversely, applying a low pass filter to the
wavepackets from case B, i.e., filtering the symmetric modes, we recover the
growth rates of anti-symmetric modes over a similar range (figures 18d-f).
Superimposing the growth rates obtained from cases A and B, we obtain a more
complete picture of the stability curves in good agreement with results from
section 4.

(a) (b)

Figure 19:  Impulse response for ($R^{*}=7$, $H^{*}=5.25$, $\beta=3/4$) for
case B: (a) contours of $\log(|z^{\prime}_{1}(\zeta,\psi)|)$; and (b)
$\sigma_{V}^{*}$ as function of $V_{\psi}^{*}$. The symmetric and anti-
symmetric parts, $z^{\prime}_{+}$ and $z^{\prime}_{-}$, propagate at different
velocities shown respectively in blue and orange.

Finally, consider the space-time evolution of the initial perturbation
corresponding to case B displayed in figure 19a. The perturbation initially
propagates in a narrow wavepacket with a high frequency content characteristic
of anti-symmetric modes, until a second wavepacket with lower frequencies
characteristic of symmetric modes is observed. Despite the overlapping of
frequencies and wavenumbers, it is still possible to estimate the front
velocities using (29) by considering separately the symmetric and anti-
symmetric parts, $z^{\prime}_{+}$ and $z^{\prime}_{-}$. The maximum growth
rate corresponds to a velocity $V_{\psi}^{*}$ close to the phase velocity
$c_{0}^{*}=c_{0}(2\pi/\Omega_{R})/t_{adv}$ of the most unstable mode (figure
19b). As expected, the value of $\sigma_{V}^{*}$ for the symmetric part is
well approximated by an equivalent helical vortex. As shown by Durán Venegas
et al. (2021), for helical vortices $\sigma_{V}^{*}$ is well approximated by
the predicted growth rate for an periodic array of point vortices

$\sigma_{V}^{*}=(\pi/2)(1-(V_{rel}^{*})^{2})$ (30)

where $V_{rel}^{*}=V_{\psi}^{*}-V_{\psi}^{\max*}$ is the velocity relative to
the vortex elements. For the range of values considered, the wave packet of
the symmetric part spreads more rapidly than that of the anti-symmetric part,
suggesting that the transition from convective to absolute instability is well
approximated by equation (30).

## References

  * Bayly (1988) Bayly, B. J. 1988 Three-dimensional centrifugal-type instabilities in inviscid two-dimensional flows. Phys. Fluids 31 (1), 56–64.
  * Bhagwat & Leishman (2000) Bhagwat, M. J. & Leishman, J. G. 2000 Stability analysis of helicopter rotor wakes in axial flight. J. Amer. Helic. Soc. 45 (3), 165.
  * Blanco-Rodríguez & Le Dizès (2016) Blanco-Rodríguez, F. J. & Le Dizès, S. 2016 Elliptic instability of a curved batchelor vortex. J. Fluid Mech. 804, 224–247.
  * Blanco-Rodríguez & Le Dizès (2017) Blanco-Rodríguez, F. J. & Le Dizès, S. 2017 Curvature instability of a curved batchelor vortex. J. Fluid Mech. 814, 397–415.
  * Blanco-Rodríguez et al. (2015) Blanco-Rodríguez, F. J., Le Dizès, S., Selçuk, C., Delbende, I. & Rossi, M. 2015 Internal structure of vortex rings and helical vortices. J. Fluid Mech. 785, 219–247.
  * Brocklehurst & Pike (1994) Brocklehurst, A. & Pike, A.C. 1994 Reduction of BVI noise using a vane tip. In AHS Aeromechanics Specialists Conference. American Helicopter Society.
  * Brynjell-Rahkola & Henningson (2020) Brynjell-Rahkola, M. & Henningson, D. S. 2020 Numerical realization of helical vortices: application to vortex instability. Theor. Comp. Fluid Dyn 34 (1), 1–20.
  * Cañadillas et al. (2020) Cañadillas, B., Foreman, R., Barth, V., Siedersleben, S., Lampert, A., Platis, A., Djath, B., Schulz-Stellenfleth, J., Bange, J., Emeis, S. & others 2020 Offshore wind farm wake recovery: Airborne measurements and its representation in engineering models. Wind Energy 23 (5), 1249–1265.
  * Castillo-Castellanos et al. (2021) Castillo-Castellanos, A., Le Dizès, S. & Durán Venegas, E. 2021 Closely spaced corotating helical vortices: General solutions. Phys. Rev. Fluids 6, 114701.
  * Crouch (1997) Crouch, J.D. 1997 Instability and transient growth for two trailing-vortex pairs. J. Fluid Mech. 350, 311–330.
  * Crow (1970) Crow, S. C. 1970 Stability theory for a pair of trailing vortices. AIAA journal 8 (12), 2172–2179.
  * Durán Venegas & Le Dizès (2019) Durán Venegas, E. & Le Dizès, S. 2019 Generalized helical vortex pairs. J. Fluid Mech. 865, 523–545.
  * Durán Venegas et al. (2021) Durán Venegas, E., Rieu, P. & Le Dizès, S. 2021 Structure and stability of Joukowski’s rotor wake model. J. Fluid Mech. 911, A6.
  * Fabre & Jacquin (2000) Fabre, D. & Jacquin, L. 2000 Stability of a four-vortex aircraft wake model. Phys. Fluids 12 (10), 2438–2443.
  * Fabre et al. (2002) Fabre, D., Jacquin, L. & Loof, A. 2002 Optimal perturbations in a four-vortex aircraft wake in counter-rotating configuration. J. Fluid Mech. 451, 319.
  * Frederik et al. (2020) Frederik, J. A., Doekemeijer, B. M., Mulders, S. P. & van Wingerden, J.-W. 2020 The helix approach: Using dynamic individual pitch control to enhance wake mixing in wind farms. Wind Energy 23 (8), 1739–1751.
  * Fukumoto & Miyazaki (1991) Fukumoto, Y. & Miyazaki, T. 1991 Three–dimensional distorsions of a vortex filament zith axial velocity. J. Fluid Mech. 222, 369–416.
  * Gupta & Loewy (1974) Gupta, B.P. & Loewy, R.G. 1974 Theoretical analysis of the aerodynamic stability of multiple, interdigitated helical vortices. AIAA J. 12 (10), 1381–1387.
  * Hardin (1982) Hardin, J. C. 1982 The velocity field induced by a helical vortex filament. Phys. Fluids 25 (11), 1949–1952.
  * Huang et al. (2019) Huang, X., Moghadam, S. M. A., Meysonnat, P.S., Meinke, M. & Schröder, W. 2019 Numerical analysis of the effect of flaps on the tip vortex of a wind turbine blade. Int. J. Heat Fluid Flow 77, 336–351.
  * Huerre & Monkewitz (1990) Huerre, P. & Monkewitz, P. A. 1990 Local and global instabilities in spatially developing flows. Annu. Rev. Fluid Mech. 22 (1), 473–537.
  * Ivanell et al. (2010) Ivanell, S., Mikkelsen, R., Sørensen, J. N. & Henningson, D. 2010 Stability analysis of the tip vortices of a wind turbine. Wind Energy 13 (8), 705–715.
  * Jimenez (1975) Jimenez, J. 1975 Stability of a pair of co-rotating vortices. Phys. Fluids 18 (11), 1580–1581.
  * Josserand & Rossi (2007) Josserand, Ch. & Rossi, M. 2007 The merging of two co-rotating vortices: a numerical study. European Journal of Mechanics-B/Fluids 26 (6), 779–794.
  * Kawada (1936) Kawada, S. 1936 Induced velocity by helical vortices. J. Aeronaut. Sci. 3 (3), 86–87.
  * Kerswell (2002) Kerswell, R. R. 2002 Elliptical instability. Annu. Rev. Fluid Mech. 34 (1), 83–113.
  * Lamb (1945) Lamb, H. 1945 Hydrodynamics. Dover publications.
  * Leishman et al. (2002) Leishman, J. G., Bhagwat, M. J. & Bagai, A. 2002 Free-vortex filament methods for the analysis of helicopter rotor wakes. J Aicr 39 (5), 759–775.
  * Levy & Forsdyke (1928) Levy, H. & Forsdyke, A. G. 1928 The steady motion and stability of a helical vortex. Proc. R. Soc. Lond. A 120, 670–690.
  * Leweke et al. (2016) Leweke, T., Le Dizes, S. & Williamson, C. H.K. 2016 Dynamics and instabilities of vortex pairs. Annu. Rev. Fluid Mech. 48 (1), 507–541.
  * Leweke et al. (2014) Leweke, T., Quaranta, H. U., Bolnot, H., Blanco-Rodriíguez, F. J. & Le Dizès, S. 2014 Long- and short-wave instabilities in helical vortices. J. Phys.: Conf. Ser. 524, 012154.
  * Meunier et al. (2002) Meunier, P., Ehrenstein, U., Leweke, T. & Rossi, M. 2002 A merging criterion for two-dimensional co-rotating vortices. Phys. Fluids 14 (8), 2757–2766.
  * Meunier & Leweke (2005) Meunier, P. & Leweke, T. 2005 Elliptic instability of a co-rotating vortex pair. J. Fluid Mech. 533, 125–159.
  * Okulov (2004) Okulov, V. L. 2004 On the stability of multiple helical vortices. J. Fluid Mech. 521, 319–342.
  * Okulov & Sørensen (2007) Okulov, V. L. & Sørensen, J. N. 2007 Stability of helical tip vortices in a rotor far wake. J. Fluid Mech. 576, 1–25.
  * Okulov & Sørensen (2009) Okulov, V. L. & Sørensen, J. N. 2009 Applications of 2d helical vortex dynamics. Theor. Comp. Fluid Dyn 24 (1-4), 395–401.
  * Quaranta et al. (2015) Quaranta, H. U., Bolnot, H. & Leweke, T. 2015 Long-wave instability of a helical vortex. J. Fluid Mech. 780, 687–716.
  * Quaranta et al. (2019) Quaranta, H. U., Brynjell-Rahkola, M., Leweke, T. & Henningson, D. S. 2019 Local and global pairing instabilities of two interlaced helical vortices. J. Fluid Mech. 863, 927–955.
  * Robinson & Saffman (1982) Robinson, A. C. & Saffman, P. G. 1982 Three-dimensional stability of vortex arrays. J. Fluid Mech. 125, 411.
  * Roy et al. (2008) Roy, C., Schaeffer, N., Le Dizès, S. & Thompson, M. 2008 Stability of a pair of co-rotating vortices with axial flow. Phys. Fluids 20 (9), 094101.
  * Schröder et al. (2020) Schröder, D., Leweke, T., Hörnschemeyer, R. & Stumpf, E. 2020 Experimental investigation of a helical vortex pair. In Deutscher Luft-und Raumfahrtkongress 2020.
  * Schröder et al. (2021) Schröder, D., Leweke, T., Hörnschemeyer, R. & Stumpf, E. 2021 Instability and merging of a helical vortex pair in the wake of a rotor. J. Phys.: Conf. Ser. 1934 (1), 012007\.
  * Selçuk et al. (2017) Selçuk, C., Delbende, I. & Rossi, M. 2017 Helical vortices: Quasiequilibrium states and their time evolution. Phys. Rev. Fluids 2 (8), 084701.
  * Walther et al. (2007) Walther, J. H., Guénot, M., Machefaux, E., Rasmussen, J. T., Chatelain, P., Okulov, V. L., Sørensen, J. N., Bergdorf, M. & Koumoutsakos, P. 2007 A numerical study of the stabilitiy of helical vortices using vortex methods. J. Phys.: Conf. Ser. 75, 012034.
  * Widnall (1972) Widnall, S. E. 1972 The stability of a helical vortex filament. J. Fluid Mech. 54 (4), 641–663.


# Time-scale characteristics of Kasaï river hydrological regime variability
for 1940-1999

Jean-Marie Tshitenge Mbuebue
Albert Mbata Muliwavyo
Vincent Lukanda Mwamba
Edmond Phuku Phuati
Albert Kazadi Mukenga Bantu
Franck Tondozi Keto

###### Abstract

The present study was undertaken with the aim of contributing to the
characterization of the nonstationary variability of the hydrological regime
of the Kasaï River using the wavelet analysis for 1940-1999. The rainfalls and
discharge over Kasaï Basin have marked fluctuations with a perceptible
downward trend and some shift around 1950, 1960, 1970, 1983 and 1994. The
results show that rainfalls over Kasaï basin and the discharge at Ilebo
station patterns exhibit a strong annual oscillation and some intermittent
oscillations in 2-8 years (1950-1975, 1983-1995) and 8-16 years (1970-1999)
time scales. The wavelet coherence analysis reveals a weak possible connection
between hydrological variables (rainfalls, discharge) and climate indices
relative to sea surface temperature and atmospheric circulation over Atlantic
tropical, Indian and Pacific Oceans (coherence less than 0.55).

_KEYWORDS_ : wavelet transform, coherence analysis, rainfall, discharge,
climate indices.

Department of Physics, University of Kinshasa, B.P. 190 Kinshasa XI, DR Congo

##  1 INTRODUCTION

The Kasaï river is the chief southern tributary of the Congo River, into
which, at Kwamouth, Congo (Kinshasa), 200 km above Malebo (Stanley) Pool, it
empties a volume approaching one-fifth that of the main stream. This stretch
is interrupted by several spectacular rapids and waterfalls, the river flowing
in deeply trenched valleys at elevations of about 600–900 m. The Kasaï River
eventually crosses a further series of rapids and falls, broadening and
deepening to make navigation possible.

For effective and sustainable management of water resources in the Kasaï basin
the knowledge of their availability and their variability is needed. It is
also an important indication for modeling and exploration of the evolution of
water resources of this region. In order to improve the understanding of the
influence of various factors on the evolution of the hydrological response of
the Kasaï, and more broadly on the development of water resources, from the
regional to a global scale, Kasaï watershed has been subject of few studies
undertaken by De Baker [1], Devroey [2] who provided a lot of information’s
about Kasaï basin. Ntombi et al. [3] and Kisangala et al. [4] have focused on
statistical analyzes based on existing data to establish an understanding of
hydrology of the basin.

In the context of climate change, the determination of preferential scales
between climate variability and hydrological variability of Kasaï watershed is
needed. The overall variability of short-term climate is generally associated
with phases of oceanic and atmospheric phenomena coupled with El Niño (ENSO)
and the North Atlantic Oscillation (NAO). While El Niño Southern Oscillation
(ENSO) affects weather and climate variability around the world, the North
Atlantic Oscillation (NAO) is the climate dominant mode in the region of the
North Atlantic. These oscillations have been used in several studies to
develop accurate models able to predict climate variability.The determination
of the climate change impact on hydrological systems and their water resources
constitutes a major issue of the 21st century for which scientists must
answer.

The main objectives are to:

  * •

identify and quantify the main modes of variability of the climate and
hydrological response of Kasaï river basin;

  * •

characterize the relationships between discharge at Ilebo’s station, gridded
rainfall over Kasaï basin signals and some climate indices.

##  2 STUDY AREA AND DATA

Kasaï basin covers an area of 904,000 square kilometers, between longitudes
15.75 ° and 24.7 ° east and between 0.75 ° South and 12.25 ° South latitude
(Fig. 1). 72.4% of this basin is in DR Congo and the remaining 27.6%,
representing the Southwest part is located in Angola.

Figure 1: Relief and stream order map of Kasaï watershed

From orographic point of view (Fig. 1), one can identified a depression,
sprinkled with few isolated massifs, at the North of the basin around the 4th
parallel which is between 300 and 500 meters of altitude. The trays exceeding
the altitude of 1000 meters are located on the western edge of the basin
between Wamba and Loange, to the south, as well as in the Southeast part. The
rapids that cut most of the major tributaries of the basin are located between
the plateau at 500 meters of altitude and the central depression.

The longest river in the southern Congo River basin system, the Kasaï River,
measures 2,153 km from its source on the eastern slope of the Bie Plateau in
Angola by 12 ° South and 19 ° East longitude Greenwich, near the plateau 1500
meters above sea level where the Zambezi River also takes its origins [2] to
Kwamouth. The Kasaï River is the largest tributary of the Congo River. It
swells waters of many tributaries whose some of them measure themselves more
than 1,000 kilometers.

Its headstream runs east for 250 miles (402 km), then turns north for nearly
300 miles (483 km) to form the frontier between Angola and Congo. Traffic is
especially heavy on the 510-mile- (820-kilometre-) long waterway from Kinshasa
to Ilebo. Above Ilebo, navigation is impeded by sandbanks, but shallow-draft
boats can reach Djokupunda. There the Mai- Munene Falls require goods to be
carried 40 miles (64 km) by rail to the next navigable reach, from Mukumbi to
Mai- Munene. In its lower reaches, the river enters the equatorial rainforest.
Below its confluence with the Kwango, it forms Wissmann Pool and then receives
the Fimi-Lukenie, bringing the waters of Lake Mai-Ndombe. Thereafter, it is
known as the Kwa.

The climate is equatorial in the northern part of Kasaï basin. It is
characterized by little difference in seasons, with high temperature and
humidity virtually constant. The rainfall distribution throughout the year is
on average around 1500 millimeters.

Figure 2: Mean annual rainfall map of the Kasaï basin (1940-1999). Data:
SIEREM.v1 at a $0.5{{}^{\circ}}\times 0.5{{}^{\circ}}$ resolution. Isohyets in
mm. Figure 3: Rainfalls minus Evapotranspiration (PRE-ETP) map over Kasaï
watershed the Kasaï basin (1940-1999). Data: Climate Research Unit CRU CL 3.10
climatology at a 0.5° resolution (New et al., 2002). Isohyets in mm.

In the South, there is the tropical climate of the Sudan type with two maxima
of rainfall and temperature; the dry season becomes longer as one progresses
from the north towards the south in the basin. The relief of the southern
region of Kasaï basin causes a slight temperature drop.

Variations in the length of the rainy season depending on the latitude are
clearly demonstrated by the graph in Fig. 6. This basin had 38 weather
stations in the colonial era. The mean maximum and minimum absolute
temperatures are respectively 32 °/27 °C and 21 °/13 °C in the north/south of
Kasaï basin (Fig. 4).

Figure 4: a) Minimum, b) Mean and c) maximum temperature over Kasaï watershed
the Kasaï basin (1940-1999). Data: Climate Research Unit CRU CL 3.10
climatology at a 0.5° resolution (New et al., 2002). Isohyets in mm.

The great equatorial forest covers just about any depression of the central
basin until about the 5th parallel; it pushes its ramifications as wide wooded
galleries, going up through the valleys of Kwilu, Kasaï and Lulua until the
7th parallel.

### Climate data

Gridded rainfall (1940-1999) were extracted from the African gridded data set
from the SIEREM database [5] provided by the HydroSciences Montpellier
Laboratory on the basis of several data sources, including the former ORSTOM
(now IRD – Institut de Recherche pour le Développement) rainfall database [6].
This data set SIEREM.v1 was constructed from SIEREM database by kriging
interpolation. SIEREM is an environmental information system for water
resources. Raw historical climate datasets are not appropriate for statistical
analysis because of inhomogeneity and leakage while gridded datasets have
undergone rigorous quality control and have been homogenized [7, 8].

Climate indices used are listed in table 1.

Table 1:  Data sources used for monthly climate indices  Climate indices | Data | Periods | Source
---|---|---|---
NAO | North Atlantic Oscillation [9] | 1950-2013 | “Earth System Research Laboratory (ESRL)
Nino 3.4 | East Central Tropical Pacific SST Index [10] | 1951-2013 | of the National Oceanic and
TNA | Tropical Northern Atlantic Index [11] | 1948-2013 | Atmospheric Administration (NOAA)”
TSA | Tropical Southern Atlantic Index [11] | 1948-2013 | web page http://www.cpc.ncep.noaa.gov/data/ind ices/
Solar | Solar Flux [12] | 1948-2014 |
DMI | Dipole Mode Index [13] | 1870-2014 | The Frontier Research Center for Global Change
|  |  | (JAMSTEC, Japan)” web page
|  |  | (http://www.jamstec.go.jp/frcgc/resear ch/d
|  |  | 1/iod/HTML/Dipole%20Mode%20In dex.html)

##  3 METHODS

Characterization of the hydrological variability (discharge), the rainfall
fluctuations and their comparison with the climate indices are performed by
using the signal processing methods (multitaper, wavelet analysis,…).

Climate variability refers to the climatic parameter of a region varying from
its long-term mean. Every year in a specific time period, the climate of a
location is different. Some years have below average rainfall, others have
average or above average rainfall.

Figure 5: Diagram of the overall approach of the study

###  3.1 HIERARCHICAL CLUSTERING ON PRINCIPAL COMPONENTS

Combining principal component methods, hierarchical clustering and partitional
clustering, help to better visualize data and highlight the main features of
the data set [14, 15].

###  3.2 THE SEGMENTATION METHOD OF HILBERT

The time series segmentation procedure in this case gives several changes of
average. By using a specific algorithm, one or more shift dates (if any) which
separate adjacent segments whose averages have a significant difference under
the Scheffe test are selected [16].

###  3.3 WAVELET TRANSFORM

Recently the wavelet transform has gained a lot of popularity in the field of
signal processing. This is due to its capability of providing both time and
frequency information simultaneously, hence giving a time-frequency
representation of the signal. The traditional Fourier transform can only
provide spectral information about a signal. Moreover, the Fourier method only
works for stationary signals. In many real world applications, the signals are
non-stationary. One solution for processing non-stationary signals is the
wavelet transform. The wavelet analysis allows the passage of a representation
of a signal to another as Fourier analysis, with a different time-frequency
resolution. The wavelet transform $W_{n}(s)$ of a discrete time series xn is
defined as the convolution of xn with a scaled and translated version of a
function l/JO called « mother wavelet » [17]:

$W_{n}\left(s\right)=\sqrt{\frac{\delta
t}{s}}\sum_{n^{\prime}=0}^{N-1}x_{n^{\prime}}\psi_{0}\left[\frac{\left(n^{\prime}-n\right)\delta
t}{s}\right]$ (1)

where $\delta t$ and $s$ are respectively the time step and the scale.

In this paper, Morlet wavelet is used as « mother wavelet » and consists of a
plane wave modulated by a Gaussian:

$\psi_{0}\left(\eta\right)=\pi^{-\frac{1}{4}}e^{i\omega_{0}\eta}e^{-\frac{\eta^{2}}{2}}$
(2)

and $\omega_{0}$ is taken to be 6 [18].

###  3.4 GLOBAL WAVELET SPECTRUM (GWS) OR TIME AVERAGE OF LOCAL WAVELET
SPECTRA

The Global Wavelet Spectrum (GWS) is used to estimate changes of wavelet power
with frequency [19]. It provides an unbiased and consistent estimation of the
true power spectrum of a time series. :

It is calculated by:

$\overline{W}^{2}\left(s\right)=\frac{1}{N}\sum_{n=0}^{N-1}\left|W_{n}\left(s\right)\right|^{2}$
(3)

###  3.5 SCALE-AVERAGED WAVELET POWER (SAWP)

The Scale-Averaged Wavelet Power (SAWP), which is defined as the weighted sum
of the wavelet power, is used to examine fluctuations in power over a range of
a scales (a band). It is given by:

$\overline{W}_{n}^{2}=\frac{\delta_{j}\delta_{t}}{C_{\delta}}\sum_{j=j_{1}}^{j_{2}}\frac{\left|W_{n}\left(s_{j}\right)\right|^{2}}{s_{j}}$
(4)

where $C_{\delta}=0.776$ for the Morlet wavelet, $\delta_{j}$ is the factor of
scale averaging, $\delta_{t}$ is the sampling period, $j_{l}$ and $j_{2}$ are
the scales over which SAWP is computed.

The wavelet –filtered reconstructions of a determined oscillation band-passed,
$j_{l}$ and $j_{2}$, is given by:

$X_{n}=\frac{\delta_{j}\delta_{t}^{\frac{1}{2}}}{C_{\delta}\psi_{0}\left(0\right)}\sum_{j=j_{1}}^{j_{2}}\frac{R\left[W_{n}\left(s_{j}\right)\right]}{s_{j}^{\frac{1}{2}}}$
(5)

where the factor $\psi_{0}\left(0\right)$ is used for removes the energy
scaling.

###  3.6 WAVELET COHERENCE ANALYSIS

The wavelet coherence, $R_{n}^{2}(s)$, is defined as the square of the cross-
spectrum normalized by the individual power spectra [20].

$R_{n}^{2}\left(s\right)=\frac{\left|s\left(s^{-1}W_{n}^{XY}\left(s\right)\right)\right|^{2}}{s\left(s^{-1}\left|s^{-1}W_{n}^{X}\left(s\right)\right|^{2}\right)s\left(s^{-1}\left|s^{-1}W_{n}^{Y}\left(s\right)\right|^{2}\right)}$
(6)

where $R_{n}^{2}(s)\in\left[0,1\right]$,
$W_{n}^{XY}\left(s\right)=W_{n}^{X}\left(s\right)\times W_{n}^{Y}(s)$ is the
cross-wavelet, $n$ represents the time index, $s$ the scale of $Y$ wavelet
transform, $W_{n}^{X}\left(s\right)$ and $W_{n}^{Y}\left(s\right)$ are the
wavelet transform of $X$ and $Y$ and $S$ denotes the smoothing operator for
$W_{n}^{XY}\left(s\right)$.

The wavelet phase difference or relative phase between two time series $X$ and
$Y$ is given as:

$\phi_{n}=\arctan\nicefrac{{\mathcal{J}\left(s^{-1}W_{n}^{XY}\left(s\right)\right)}}{{\mathcal{R}\left(s^{-1}W_{n}^{XY}\left(s\right)\right)}}$
(7)

where $\mathcal{J}\left\\{.\right\\}$ and $\mathcal{R}\left\\{.\right\\}$ are
the imaginary and real parts of the wavelet spectra. Wavelet coherence lies
between 0 and 1. This reflects the time–scale variability of the linear
relationship between two signals. A near-zero coherence indicates no linear
relationship between the signals. Coherence close to 1 indicates a linear
relationship between both signals.

##  4 RESULTS

Hierarchical Clustering on Principal Components Kasaï basin rainfalls shows
that the Kasaï basin can be divided into three areas as in fig. 7.

Figure 6: Typology of rainfall regimes over the Kasaï basin. Based on a
hierarchical cluster analysis of mean monthly rainfall for 1940-1999
(SIEREM.v1data). Figure 7: Map of the 3 regional divisions of Kasaï basin
after hierarchical clustering methods on principal component analysis of
rainfall Figure 8: Percentage of variance explained

Only PC1 will be discussed as it contributes for 80 % of the rainfall
variability in the Kasaï basin (figure 8).

Figure 9: Phases of homogeneous rainfalls over area 1,2 and 3 of Kasaï basin
(1940-1999)

By Hubert segmentation, which gives us the year of different ruptures observed
on the rainfall time series (Fig. 9), we can see that rains on the Kasaï basin
are marked by fluctuations with a tendency on the one hand to the noticeable
decline between 1970 and 1990, 1993 and 1999 and also an upward trend between
1960 and 1970 in area 1, 1990 and 1992.

The rain in area 1 presents five ruptures and therefore six times homogeneous
rainfall:

wet phases 1940-1959; 1959-1969; 1970-1979; 1988-1990; dry phases 1980-1987;
1991-1999;

The rain in area 2 presents five ruptures and therefore six times homogeneous
rainfall:

wet phases 1940-1969; 1988-1990; 1996-1999; dry phases 1970-1975; 1976-1987;
1991-1995;

The rain in area 3 presents three ruptures and therefore four times
homogeneous rainfall:

wet phases 1940-1959; 1970-1979; dry phases 1959-1969; 1980-1987;

###  4.1 IDENTIFICATION OF VARIABILITY MODES

The wavelet analysis is used to study the frequency composition and visualize
unsteady fluctuations of hydro variables so as to identify the main modes of
variability of flows and precipitation in the basin of the Kasaï.

Table 1 shows the contribution of the energy bands in the discharge at Ilebo
station and in rainfall in the Kasaï basin. The wavelet analysis of Kasaï flow
shows a plurality of energy bands (Figure 10 and 11): A period of about 1 year
mainly expressed A band of 2-8 years well localized:

  * •

between 1953 and 1973 for both rainfalls and discharge;

  * •

between 1973 and 1995 for rainfalls

A band of 8-16 years with a structure which begin around 1973 till 1999.

One can find a change point around 1970 (Figure 10 and 11) which is associated
to low rainfalls and discharge.

Figure 10: Continuous wavelet power spectrum of a) discharge of Kasaï River at Ilebo station and b) PC 1 of gridded rainfall over Kasaï basin. The left axis is the Fourier period (in yr) corresponding to the wavelet scale on the right axis. The bottom axis is time (yr). The thick black contour designates the 5% significance level against red noise and the cone of influence (COI) where edge effects might distort the picture is shown as a lighter shade. Figure 11: Continuous global wavelet power spectrum of a) discharge of Kasaï River at Ilebo station and b) PC 1 of gridded rainfall over Kasaï basin. Table 2: Contribution of selected energy bands on total variance of discharges for the Kasaï River at Ilebo station and rainfalls in Kasaï basin | Discharge | Rainfall
---|---|---
| Ilebo | Region 1 | Region 2 | Region 3
0-1 year | 42% | 44% | 45% | 45%
1-2 years | 34% | 35% | 30% | 25%
2-8 years | 10% | 5% | 6% | 8%
8-16 years | 6% | 3% | 4% | 6%

###  4.2 COHERENCE BETWEEN RAINFALLS AND DISCHARGES

The discharge of Kasaï River at Ilebo station seems to be strongly related to
the rainfalls with total coherence around 59.63%. Figure 12 attests that the
coherence between discharge and rainfall in the Kasaï basin are essentially
distributed across one year cycle and in 2-8 band (6.5 years).

Figure 12: Wavelet coherence between the discharge of Kasaï River at Ilebo
station and the PC 1 (80%) of rainfalls over Kasaï basin. The 5% significance
level against red noise is shown as a thick contour. The relative phase
relationship is shown as arrows (with in-phase pointing right, anti-phase
pointing lefting straight down).

###  4.3 RELATIONSHIP TO CLIMATE FORCING

Some characteristic oscillations were next related to several well-known
climate indices. The coherence wavelet analyses between discharge time series
at Ilebo station, rainfalls (Figure 12) in the three area of Kasaï basin and
climate index (NAO, NINO 3.4, TNA, TSA, DMI and Solar flux) show period of
high coherence during some periods listed in table 3 and 4. Figure 13 and 14
also shows the characteristic periods of high wavelet coherence between
discharge at Ilebo station, first principal component of rainfalls in Kasaï
basin and selected climate indices.

Figure 13: Wavelet coherence between the discharge of Kasaï River at Ilebo station and climate indices (NAO, Nino 3.4, TNA, TSA, DMI, Solar flux. The 5% significance level against red noise is shown as a thick contour. Table 3:  Contribution of selected energy bands on total variance of discharges for the Kasaï River at Ilebo station and rainfalls in Kasaï basin Observed scales | NAO | Nino 3.4 | TNA | TSA | DMI | Solar
---|---|---|---|---|---|---
2-4 years | 1951-1955 |  | 1982-1987 | 1982-1987 | 1960-1969 |
| 1965-1971 |  |  |  |  |
| 1979-1990 |  |  |  |  |
4-8 years | 1955-1958 | 1932-1945 |  |  | 1950-1957 | 1948-1962
|  | 1952-1962 |  |  | 1980-1990 |
8-16 years | 1950-1960 | 1948-1972 | 1967-1985 | 1967-1985 |  |

Figure 14: Wavelet coherence between the PC 1 of rainfalls over Kasaï basin and climate indices (NAO, Nino 3.4, TNA, TSA, DMI, Solar flux. The 5% significance level against red noise is shown as a thick contour. The relative phase relationship is shown as arrows (with in-phase pointing right, anti-phase pointing left, and BMI leading AO by 90◦ pointing straight down). Table 4: Contribution of selected energy bands on total variance of discharges for the Kasaï River at Ilebo station and rainfalls in Kasaï basin  Observed scales | NAO | Nino 3.4 | TNA | TSA | DMI | Solar flux
---|---|---|---|---|---|---
2-4 years |  | 1955-1957 | 1957-1960 |  | 1940-1948 | 1956-1961
| 1964-1975 | 1967-1976 |  | 1951-1952 |
|  | 1979-1984 |  | 1970-1978 |
|  | 1982-1992 |  |  |
4-8 years | 1954-1967 | 1950-1967 | 1962-1972 |  |  |
8-16 years | 1975-1995 | 1967-1999 | 1957-1999 | 1972-1999 |  | 1980-1997
Figure 15: Comparison between the wavelet –filtered reconstructions of the
monthly discharge over 2-8 band-passed of Kasaï river at Ilebo station and
Nino 3.4 index. Figure 16: Comparison between the wavelet –filtered
reconstructions of the monthly rainfall over 2-8 band-passed on area 1 of
Kasaï basin and Nino 3.4 index.  Figure 17: Comparison between the wavelet
–filtered reconstructions of the monthly rainfall over 2-8 band-passed on area
2 of Kasaï basin and Nino 3.4 index.  Figure 18: Comparison between the
wavelet –filtered reconstructions of the monthly rainfalls over 2-8 band-
passed on area 3 of Kasaï basin and Nino 3.4 index.

Figure 19 presents the relationships over the 2–8-yr band between the first
principle component of the gridded rainfall over Kasaï basin and the Niño3.4
SST index.

Figure 19: Percentage of scale-averaged wavelet coherence over the 2–8-yr band
for the gridded rainfall over Kasaï basin and the Niño 3.4 SST index.

Mean percentage of wavelet coherence between discharge at Ilebo station or
rainfall and some selected climate indices except Nino 3.4 index shows a weak
nonlinear relationship (figure 20).

Figure 20: Mean percentage of wavelet coherence ( 2-8 years band-passed)
between of Ilebo’s discharge, PC 1 of gridded rainfalls, Monthly mean
rainfalls in Kasaï’s area 1, area 2, area 3 and some climate indices (NAO,
NINO 3.4, TNA, TSA, DMI and Solar flux)

##  5 5 DISCUSSIONS AND CONCLUSION

We propose in this article a time-scale analysis of the river discharge
fluctuations at Ilebo and rainfalls over the Kasaï basin. The two temporal
discontinuities around 1970-1975 and 1990-1995 over figure 9 and 10 were also
found in other studies [8], [20]. The wavelet coherence between discharges at
Ilebo and rainfall over Kasaï River shows a high correlation in annual and the
2-8 year band –passed. Wavelet coherence between discharges at Ilebo and
rainfall over Kasaï River and selected climate index confirm a less (figure
20) and nonlinear relationship between the hydrology of the Kasaï basin and
the Pacific, Atlantic and Indian SST fluctuations. One can see covariability
between rainfalls in the Kasaï basin and ENSO phenomenon, particularly at El
Niño years [10] (Figure 15-19). The influence of teleconnections is low but
their effects seem especially important having regard to the covariability
observed with El Niño example. Further investigations should be conducted to
understand the variability of the Kasaï river basin regime variability. What
is the role of deforestation and of human activity and urbanization on water
resource to availability?

## ACKNOWLEDGMENTS

Thanks for the gridded rainfall data providers IRD and HydroSciences
Montpellier, France. The wavelet software was provided by:

Torrence and Compo at http://atoc.colorado.edu/research/wavelets/

A. Grinsted at http://www.pol.ac.uk/home/research/waveletcoherence/

## References

  * [1] D. Backer, “La rivière kasaï entre kwamouth et port - francqui. l amélioration de ses conditions de navigabilité par balisage, signalisation et travaux dâaménagement,” _bulletin des séances de l’Institut Royal Colonial Belge_ , vol. III, no. 1, pp. 252 – 278, 1932.
  * [2] E. Devroey, Le Kasaï et son bassin hydrographique. Goemaere, 1939.
  * [3] Ntombi and Kisangala, “Impact de la lithologie et de lâhydrométrie sur la navigabilité du kasaï en r. d. congo,” pp. 157–164, 2002.
  * [4] Kisangala and Ntombi, “Estimation des réserves régulatrices et l’établissement de la courbe de tarissement du bassin versant du kasaè durant la période dâétiage à lâéchelle de lumbu,” _Revue Congolaise des Sciences Nucléaires_ , vol. 26, no. 2, 2012.
  * [5] [5] J.-F. Boyer, C. Dieulin, N. Rouche, A. Cres, E. Servat, J.-E. Paturel, and G. Mahe, “Sierem : an environmental information system for water resources,” _IAHS publication_ , vol. 308, p. 19, 2006.
  * [6] [6] N. Rouché, G. Mahé, S. Ardoin-Bardin, B. Brissaud, J.-f. Boyer, A. Crès, C. Dieulin, G. Bardin, G. Commelard, A. Dezetter, et al., “Constitution d’une grille de pluies mensuelles pour l’afrique (période 1900-2000),” _Science et changements planétaires/Sécheresse_ , vol. 21, no. 4, pp. 336– 338, 2010.
  * [7] G. Mahé, “Modulation annuelle et fluctuations interannuelles des précipitationscipitations sur le bassin-versant du congo,” pp. 13–26, 1995.
  * [8] G. Mahé, Y. Lhote, J. C. Olivry, and G. Wotling, “Trends and discontinuities in regional rainfall of west and central africa : 1951 - 1989,” _Hydrological Sciences Journal_ , vol. 46, no. 2, pp. 211 – 226, 2001.
  * [9] J. W. Hurrell, “Decadal trends in the north atlantic oscillation : regional temperatures and precipitation,” _Science_ , vol. 269, no. 5224, pp. 676–679, 1995.
  * [10] K. E. Trenberth, “The definition of el nino,” _Bulletin of the American Meteorological Society_ , vol. 78, no. 12, pp. 2771–2777, 1997
  * [11] D. Enfield and E. Alfaro, “The dependence of caribbean rainfall on the interaction of the tropical atlantic and pacific oceans,” _Journal of Climate_ , vol. 12, p. 2093â2103, 1999.
  * [12] “Latest solar radio flux report from drao, penticton,â natural resources canada, 10 february 2016. [online]. available : http ://www.spaceweather.ca/solarflux/sx-4-en.php.”
  * [13] K. Ashok, Z. Guan, and T. Yamagata, “Impact of the indian ocean dipole on the relationship between the indian monsoon rainfall and enso,” _Geophysical Research Letters_ , vol. 28, no. 23, pp. 4499–4502, 2001.
  * [14] S. Le Dien and J. Pagès, “Hierarchical multiple factor analysis : Application to the comparison of sensory profiles,” _Food quality and preference_ , vol. 14, no. 5, pp. 397–403, 2003.
  * [15] L. Lebart, A. Morineau, and K. M. Warwick, “Multivariate descriptive statistical analysis ; correspondence analysis and related techniques for large matrices,” 1984.
  * [16] P. Hubert, J. P. Carbonnel, and A. Chaouche, “Segmentation des séries hydrométéorologiques : application à des séries de précipitations et de débits de l’afrique de l’ouest,” _Journal of hydrology_ , vol. 110, no. 3, pp. 349–367, 1989\.
  * [17] C. Torrence and G. Compo, “A practical guide to wavelet analysis,” _Bulletin of the American Meteorological Society_ , vol. 79, pp. 61 – 78, 1998\.
  * [18] M. Farge, “Wavelet transforms and their applications to turbulence,” _Annual Review of Fluid Dynamics_ , vol. 24, pp. 395 – 457, 1992.
  * [19] Y. Hao, G. Liu, H. Li, Z. Li, J. Zhao, and T.-C. J Yeh, “Investigation of karstic hydrological processes of niangziguan springs (north china) using wavelet analysis,” _Hydrological Processes_ , vol. 26, no. 20, pp. 3062–3069, 2012\.
  * [20] [20] A. Laraque, G. Mahé, D. Orange, and B. Marieu, “Spatiotemporal variations in hydrological regimes within central africa during the xxth century,” _Journal of Hydrology_ , vol. 245, no. 1, pp. 104–117, 2001\.