A Recurrence Analysis of Multiple African Easterly Waves during Summer 2006

Abstract

Accurate detection of large-scale atmospheric tropical waves, such as African easterly waves (AEWs), may help extend lead times for predicting tropical cyclone (TC) genesis. Since observed AEWs have comparable but slightly different periods showing spatial and temporal variations, local analysis of frequencies and amplitudes of AEWs is crucial for revealing the role of AEWs in the modulation of TC genesis. To achieve this goal, we investigate the recurrence plot (RP) method. A recurrence is defined when the trajectory of a state returns to the neighborhood of a previously visited state. To verify implementation of the RP method in Python and its capability for revealing a transition between different types of solutions, we apply the RP to analyze several idealized solutions, including periodic, quasiperiodic, chaotic and limit cycle solutions, and various types of solutions within the three- and five-dimensional Lorenz models. We then extend the RP analysis to two datasets from the European Centre for Medium-Range Weather Forecasts global reanalysis and global mesoscale model data in order to reveal the recurrence of multiple AEWs during summer 2006. Our results indicate that the RP analysis effectively displays the major features of time-varying oscillations and the growing or decaying amplitudes of multiple AEWs.

1. Introduction

Recent studies suggest that accurate detection of recurrent, multiple, large-scale tropical waves, such as African easterly waves (AEWs), has the potential to help extend prediction lead times for tropical cyclone (TC) genesis. For example, using the NASA global mesoscale model (GMM, e.g., [1]), the formation, subsequent intensification, and movement of hurricane Helene (2006) were simulated to a degree of satisfaction from Day 22 to 30 within a 30-day model integration when multiple AEWs were realistically simulated (e.g., [2]). To date, the scalable parallel ensemble empirical mode decomposition method (PEEMD; e.g., [3, 4]) has revealed the role of downscaling processes associated with environmental flows in determining the timing and location of hurricane Helene’s formation (e.g., [5]), supporting the model’s practical predictability over extended-range time scales. Prior to Helene’s formation, observed AEWs had comparable but slightly different periods with both spatial and temporal variations. Thus, a local analysis of the frequencies and amplitudes for multiple AEWs is desired in order to monitor the evolution of AEWs that may influence the timing and location of TC genesis.

To effectively reveal space-varying features and/or temporal oscillations in simulations and reanalysis data, recurrence plots (RPs) are employed (e.g., [6]). Recurrence is defined when the trajectory of a state returns back to the neighborhood of a previously visited state. Thus, recurrence may be viewed as a generalization of periodicity to brace quasiperiodicity [7] and chaos [8]. Measuring the patterns associated with recurrence provides valuable information regarding the oscillatory nature of a system. To verify implementation in Python, we test the RP method using several types of idealized solutions and three types of solutions within the three- and five-dimensional Lorenz model (3DLM and 5DLM). We then apply the method in order to analyze two global datasets.

The paper is organized as follows. Section 2.1 introduces the European Centre for Medium-Range Weather Forecasts (ECMWF) global reanalysis data and NASA GMM data (e.g., [1, 2, 4, 9]). We focus on a 30-day period in 2006 over which multiple AEWs were observed during the NASA African Monsoon Multidisciplinary Analyses (NAMMA) field campaign (e.g., [2, 10]). Section 2.2 introduces the 3DLM and 5DLM (e.g., [11, 12, 13]). Section 2.3 introduces RP methods for performing the local analysis. In Section 3.1, the method is applied in order to analyze four basic types of solutions, including (1) periodic, (2) quasiperiodic, (3) chaotic solutions, and (4) Gaussian white noise. While Section 3.2 presents an analysis of the idealized spiral sink, Section 3.3 extends the RP analysis for steady-state solutions of the 3DLM. Section 3.4 discusses limit cycle solutions of the 3DLM and 5DLM and their RP analysis. Section 3.5 presents an analysis of global reanalysis and GMM data in order to reveal the time-varying recurrent behavior of AEWs. Appendix A discusses the mathematical approach for the unique version of the 3DLM and 5DLM, referred to as the Version 2 (V2) system. V2 systems are used for either linear or nonlinear simulations. Additionally, the linearized V2 system is applied to perform an eigenvalue analysis to obtain the growth rates and oscillatory frequencies of solutions.

2. Global data, the Lorenz model, and the recurrence analysis method

2.1 Global reanalysis and model data

For this study, we select the latest European Centre for Medium-Range Weather Forecasts (ECMWF) global reanalysis (ERA-interim) dataset. The data spans the time period from January 1979 to the present day. The dataset has a sampling rate of 6 h and a horizontal grid spacing of $0.75°$ , yielding a spatial resolution of approximately 78 km. The following gridded products are available: 3-h surface fields, daily vertical integrations, 6-h upper air atmospheric fields for pressure, potential temperature, potential vorticity levels, and vertical coverage from the lower troposphere to the stratosphere on 60 model layers. Based on previous studies (e.g., [2, 5]), here, we focus on the 30-day period from August 22 to September 20, 2006. During this period, the NASA African Monsoon Multidisciplinary Analyses (NAMMA) field campaign documented multiple AEWs and provided a great opportunity to characterize the frequency of AEWs as well as their evolution and impact over continental western Africa. Figure 1a and b display multiple AEWs that were detected by a change in sign of meridional winds from the ERA-interim dataset and the NAMMA observations, respectively.

Figure 1.

A time-altitude cross section of meridional winds from the ERA-interim dataset (a) and global mesoscale model simulations (c). The color bar to the right of each figure represents the wind velocity in meters per second. Panel (b) displays NAMMA observations for a comparison (e.g., [10]).

As a result of higher spatial and temporal resolutions, GMM simulation data [1, 9] are also selected. The GMM that produced the dataset is composed of three major components: finite volume dynamics [14], National Center for Atmospheric Research (NCAR) Community Climate Model physics, and the NCAR Community Land Model. Simulations spanning the 35 day period began on August 22, 2006. However, our analysis only focuses on the first 30 days. The data set has a 15-min integration time step, 17 vertical levels, and a horizontal resolution of approximately 28 km. Figure 1c provides the time-altitude cross section of meridional winds from the GMM dataset.

2.2 The five-dimensional Lorenz model (5DLM)

In the real world, as a result of nonlinearity, system solutions may possess time-varying frequencies and amplitudes. To understand the performance of the selected method in analyzing nonlinear time series data, the 3DLM and 5DLM are used to provide three types of solutions at different ranges of parameters for verifications. The 5DLM, which extends the 3DLM and includes the 3DLM as a subset (e.g., [13]), is defined by Eqs. (1)–(5), as follows:

$$\frac{\mathit{dX}}{\mathit{d\tau }}=-\mathit{\sigma X}+\mathit{\sigma Y},$$ E1

$$\frac{\mathit{dY}}{\mathit{d\tau }}=-\mathit{XZ}+\mathit{rX}-Y,$$ E2

$$\frac{\mathit{dZ}}{\mathit{d\tau }}=\mathit{XY}-{\mathit{XY}}_1-\mathit{bZ},$$ E3

$$\frac{{\mathit{dY}}_1}{\mathit{d\tau }}=\mathit{XZ}-2{\mathit{XZ}}_1-d_o{Y}_1,$$ E4

$$\frac{{\mathit{dZ}}_1}{\mathit{d\tau }}=2{\mathit{XY}}_1-4{\mathit{bZ}}_1.$$ E5

Here, $\tau$ is the dimensionless time. In this study, the following parameters are kept constant: $b=8/3$ , $d=19/3$ , and $\sigma =10$ . By comparison, we vary the value of the Rayleigh parameter, $r$ , that represents a heating or forcing term. The term ${\mathit{XY}}_1$ in Eq. (3) plays a role in providing negative nonlinear feedback associated with small-scale dissipative processes (e.g., $-d_o{Y}_1$ and $-4{\mathit{bZ}}_1$ ) to stabilize the solutions (e.g., [13]). Equations. (1)–(3), without inclusion of the feedback term ( ${\mathit{XY}}_1$ ), are reduced to become the 3DLM. As discussed below, the 5DLM requires higher values of $r$ for the onset of chaos (e.g., [13]). Recently, the 5DLM has been re-derived and analyzed in detail by Moon et al. [15] and Felicio and Rech [16], showing the model’s robustness. Furthermore, the 5DLM has been extended to higher-dimensional LMs (e.g., [17, 18, 19, 20]) and a generalized LM (e.g., [21, 22]).

The 3DLM produces three types of solutions, including steady-state solutions, chaotic solutions, and nonlinear periodic solutions. We present the first two types of solutions below and discuss the third type of solution in Section 3.4. Figure 2a provides $Y$ and $Z$ modes of the solution for the strange attractor using the 3DLM with $r=25$ . Irregular oscillations surrounding nontrivial critical points that are defined in Appendix A are indicative of chaotic solutions. In Figure 2b, the 5DLM with $r=25$ clearly produces a steady-state solution that spirals into the non-trivial critical point. Figure 2c from the 5DLM with $r=43.5$ appears more like Figure 2a and depicts an irregular oscillatory motion. Such behavior within the 5DLM only occurs when $r\ge r_c=42.9$ and when the other parameters are held at the aforementioned values. The corresponding solutions for Y1 and Z1 are shown in Figure 2d.

Figure 2.

Phase space plots for the 3DLM (a) and 5DLM (b–d). All of the plots show Lorenz strange attractors. (a) Y vs. Z modes using the 3DLM with $r=25$ . (b and c) Y vs. Z modes using the 5DLM with $r=25$ and $r=43.5$ , respectively. (d) $Y_1$ vs. $Z_1$ modes with $r=43.5$ (see details in [13]).

2.3 Recurrence plot analysis

The recurrent nature of a system can offer important insight into its dynamics, such as periodicities or other oscillatory behavior (e.g., growing or decaying oscillations). In the case of deterministic systems, the recurrence of states (i.e., a state returning to an arbitrarily close area at a point from a previous time within the phase space) is an important feature (e.g., [6]). Such recurrences can be visualized using a tool known as recurrence plots (RPs). RPs allow for the representation of multidimensional systems in a two-dimensional visualization; thus, they easily provide insight into high-dimensional systems. The RP plots a black point for each time $\left(i,j\right)$ in which there is a recurrence, more precisely given by the following equation:

$$R_{i,j}=\mathrm{\Theta }\left(\epsilon -\parallel \overrightarrow{x_i}-\overrightarrow{x_j}\parallel \right)i,j=1,\dots ,N,$$ E6

where $N$ is the number of states ( $\overrightarrow{x_i}$ ) considered, $\epsilon$ is the threshold distance for recurrence, $\parallel \cdot \parallel$ is the Euclidean norm, and $\mathrm{\Theta }\left(\cdot \right)$ is the Heaviside function (e.g., [23]). The individual point $\left(i,j\right)$ does not offer any information regarding the states of the system at times $i$ and $j$ . However, the phase space trajectory can be reconstructed from the collection of all of the recurrence points present within the RP (e.g., [6]). As shown in Figure 3a, a typical RP for a simple periodic solution displays uniform diagonal lines. Here, it should be noted that the RP has a main black diagonal line with an angle of $\pi /4$ , referred to as the line of identity (LOI).The threshold distance for recurrence ( $\epsilon$ ) is an important parameter to consider. If the value is too small, there may not be enough recurrences plotted, resulting in insufficient information required for conducting an accurate analysis. For a value that is too large, the appearance of many recurrent points may mask important structures that would otherwise be present within the RP. While there is much discussion regarding the proper selection of $\epsilon$ , a general rule is to find the smallest possible value that is capable of providing a sufficient number of recurrences for identifying recurrent structures within a dataset (e.g., [6]). Figure 2.10 of Reyes [24] demonstrates the impact with different values of $\epsilon$ on the RP for the 3DLM.

Figure 3.

Recurrence plots of (a) periodic solutions produced using $y=\mathit{sin}\left(2\mathit{\pi t}/4\right)$ , (b) quasiperiodic solutions produced using $y=\mathit{sin}\left(16\mathit{\pi t}\right)+\mathit{sin}\left(8\sqrt{\left(2\right)}\mathit{\pi t}\right)$ , (c) chaotic solutions produced by the 3DLM, and (d) Gaussian white noise.

3. Results

In this section, we first discuss the RP analysis for idealized solutions, various types of solutions from the 3DLM and 5DLM. We then apply the RP to analyze the global reanalysis and the GMM 30-day simulation data in order to detect recurrent, multiple AEWs that play an important role in the modulation of TC genesis (e.g., Hurricane Helene (2006)).

3.1 Analysis of four basic types of solutions

Figure 3 provides RPs for four types of solutions, including (a) periodic solutions that are characterized by uniform diagonal lines extending from edge to edge of the plot and spaced with an equal distance from one another; (b) quasiperiodic solutions with a similar structure to the periodic RP, with the exception that the distance between the diagonal lines varies; (c) chaotic solutions from the 3DLM with a Rayleigh parameter of $r=28$ whose plots have many short diagonal lines spaced various distances from one another; and (d) Gaussian white noise, characterized by a homogenous RP consisting of many individual points and little to no organized recurrent structures (e.g., [6]).

Other features of the RP analysis, as discussed below, include (1) vertical or horizontal lines or clusters, which indicate small or no change in solution amplitudes over a certain period of time, and (2) white bands, which suggest transitions from one type of dynamics to another. More than one distinct line or dot structure may also appear within a single RP. Line structures and the length and number of lines, as well as their positioning, describe the dynamics for each system (e.g., [25]). When applying this method to real data, the overall structures should still be visible, although they may or may not be as neatly depicted.

Previously, performance of the RP was demonstrated using idealized data. In practice, most real data has some amount of random noise. Therefore, understanding the capabilities of the method in analyzing noisy data, before applying them to real or non-idealized datasets, is important. Here, we present an analysis of composite data consisting of the raw data and additional noise. Figure 4a displays the RP for idealized periodic solutions generated from $y=\mathit{sin}\left(2\mathit{\pi t}/16\right)$ , $0\le t\le 128$ , and $\mathrm{\Delta }t=0.2$ . Figure 4b provides the RP for the same dataset superimposed with Gaussian white noise with an amplitude of 0.1. While the overall diagonal line structure is still present, some visual differences are apparent. The most prominent differences include thicker diagonal lines and a fluctuation in thickness. The differences are due to so-called false recurrence resulting from the noise itself and the slightly higher value of $\epsilon$ required in order to obtain a reasonable plot.

Figure 4.

Recurrence plots of raw data with a periodic solution given by $y=\mathit{sin}\left(2\mathit{\pi t}/16\right)$ (a) and the raw data superimposed with Gaussian white noise that has a finite amplitude of 0.1 (b).

3.2 Analysis of a spiral sink and source

In addition to the four fundamental solutions in Figure 3, real-world model data may include growing or decaying oscillatory solutions. Mathematically, such a solution can be represented as $y=e^{\mathit{\alpha t}}sin\left(\mathit{\beta t}+\theta \right)$ , where $\theta$ is a phase constant. This type of solution can be obtained from a linearized system with a complex eigenvalue, $\lambda =\alpha +\mathit{i\beta }$ , where $\alpha$ and $\beta$ are real numbers. Next, we discuss the conditions under which a recurrence plot can be defined for the general oscillatory solution. To facilitate discussions, we further assume $\beta >0$ and consider three cases with (i) $\alpha =0$ for a simple oscillation, (ii) $\alpha >0$ for a growing oscillation (i.e., a spiral source), and (iii) $\alpha <0$ for a decaying oscillation (i.e., a spiral sink), respectively.

RPs for spiral sinks or sources display a different structure as compared to RPs for periodic and quasiperiodic solutions, as discussed below. Figure 5a displays an oscillatory solution, $y=e^{-0.1t}\mathit{sin}\left(t\right)$ (i.e., $\left(\alpha ,\beta ,\theta \right)=\left(-\mathrm{0.1,1,0}\right)$ ). The solution decays with time over the time interval from 0 to 128. When the solutions produce time-varying amplitudes, the RPs display vertical and horizontal lines. As the solution decays toward zero, the lines become denser (the upper right-hand corner of Figure 5c). Namely, the density of the lines in these RPs increases when the solution is closer to the mean value (i.e., zero for this idealized solution). The appearance of horizontal and vertical lines is explained below.

Figure 5.

A time evolution of the solution $y=e^{-01.t}\mathit{sin}\left(t\right)$ (a), its power spectrum (b), and RP (c) for $\tau \in \left[\mathit{0},\mathit{128}\right]$ .

Here, we provide an additional analysis for the pattern of the RP in Figure 5. As shown in Figure 6, a distance ( $D$ ) between two points $y\left(t_1\right)$ and $y\left(t_2\right)$ with a time lag of $T=\frac{2\pi }{\beta }$ is given by $D=\mid y\left(t_2\right)-y\left(t_1\right)\mid$ , where $t_2=t_1+T$ . We can obtain the following equation:

Figure 6.

An oscillatory mode and two points at $t_1$ and $t_2$ , where $t_2=t_1+T$ .

$$D=\mid e^{\alpha t_2}sin\left(\beta t_2+\theta \right)-e^{\alpha t_1}sin\left(\beta t_1+\theta \right)\mid$$ E7

$$\begin{array}{l}=|e^{\alpha \left(t_1+\frac{2\pi }{\beta }\right)}sin\left(\beta \left(t_1+\frac{2\pi }{\beta }\right)+\theta \right)-e^{\alpha t_1}sin\left(\beta t_1+\theta \right)|\\ =|e^{\frac{2\mathit{\alpha \pi }}{\beta }}{e}^{\alpha t_1}sin\left(\beta t_1+\theta \right)-e^{\alpha t_1}sin\left(\beta t_1+\theta \right)|\\ =|e^{\frac{2\mathit{\alpha \pi }}{\beta }}{e}^{\alpha t_1}sin\left(\beta t_1+\theta \right)-e^{\alpha t_1}sin\left(\beta t_1+\theta \right)|\\ =|e^{\alpha t_1}sin\left(\beta t_1+\theta \right)\left(e^{\frac{2\mathit{\alpha \pi }}{\beta }}-1\right)|.\end{array}$$ E8

Recurrence appears when the distance, $D$ , is small (i.e., $D<\epsilon$ , where $\epsilon$ represents a threshold), requiring

$$|e^{\alpha t_1}sin\left(\beta t_1+\theta \right)\left(e^{\frac{2\mathit{\alpha \pi }}{\beta }}-1\right)|<\epsilon .$$ E9

Therefore, (i) for a simple oscillation, $\alpha =0$ , Eq. (8) is always valid, suggesting that two points with a time lag of $T$ define a recurrent point within the recurrence plot. For a specific $t_1$ where $sin\left(\beta t_1+\theta \right)=1/M\ne 0$ , Eq. (8) leads to:

$$\mid e^{\alpha t_1}\mid <\left|\frac{M}{\left(e^{\frac{2\mathit{\alpha \pi }}{\beta }}-1\right)}\right|\epsilon .$$ E10

The value of $\alpha t_1$ roughly determines if the above inequality is valid. As a result, Eq. (10) suggests that (ii) for a spiral source, $\alpha >0$ , a recurrence point appears when $t_1$ is small and that (iii) for a spiral sink, $\alpha <0$ , a recurrence point appears when $t_1$ is large. This type of RP appears in the red box in Figure 5c. Additionally, (iv) when $sin\left(\beta t_1+\theta \right)=0$ in Eq. (7a), the point at $t=t_1$ can produce recurrence with all of the points at $t_2$ that have very small amplitudes (i.e., the amplitudes are close to zero). Thus, as shown in Figure 5c, the appearance of continuous horizontal lines indicates a transition from the third type of solution (i.e., a spiral sink) into the 4th type of solution with a small or zero amplitude. Note that the distance of two horizontal lines determined in (iv) yields an estimate of a half of the period (i.e., $T/2$ .) In general, Eq. (9) suggests that the aforementioned recurrence threshold should be selected by taking the rate of growth or decay of the oscillatory solutions into account, which will be the subject of a future study.

3.3 Analysis of a steady-state solution within the 3DLM

While the 3DLM with $r>24.74$ produces chaotic solutions, it simulates steady-state solutions when $r<24.74$ . Since the steady-state solution displays a decaying oscillatory mode, which may resemble a weakening AEW, its RP is analyzed below. Figure 7a displays numerical solutions of the Y mode of the 3DLM with $r=10$ . At approximately $t=7.0$ , only small fluctuations in the solution remain (i.e., the solution has a small amplitude over this interval). Figure 7b displays the corresponding recurrence plot with $\epsilon =0.05$ . Initially, horizontal and vertical lines appear, consistent with the previous analysis using Figures 5 and 6. As the solution begins to converge to a steady state, the plot rapidly becomes denser.

Figure 7.

Solutions for the Y mode of the 3DLM with $r=\mathit{10}$ (a) and the corresponding recurrence plot (b). Power spectra for Y mode solutions of the 3DLM V2, $\mathit{FN}=0$ (c) and $\mathit{FN}=\mathit{1}$ (d).

By comparison, Figure 7c displays the power spectrum for the linear 3DLM V2 ( $\mathit{FN}=0$ ) with the same parameter values as the solution in Figure 7a. Eigenvalues for the 3DLM V2 $\mathit{FN}=0$ are approximately $-0.59550\pm 6.17416i,-12.47567$ , indicating that the frequency of a linear solution near the nontrivial critical point is $6.174161/2\pi =0.9826$ , as confirmed by the power spectrum. However, such a spectral analysis does not reveal a local transition from a decaying oscillation to a constant solution. Figure 7d provides the power spectrum for the corresponding nonlinear simulations with $\mathit{FN}=1$ within the 3DLM V2, which produces a result similar to panel (c). Comparable spectra for linear and nonlinear steady-state solutions are due to the fact that perturbations (which measure departures from the nontrivial critical point) become smaller as time proceeds when solutions move closer to the critical point.

3.4 Analysis of limit cycle solutions

At moderate heating parameters, both the 3DLM and 5DLM models produce chaotic solutions when other parameters are held constant. By comparison, for very large values of $r$ , the 3DLM and 5DLM produce limit cycle solutions. A limit cycle is defined as a closed and isolated orbit to which nearby trajectories converge (e.g., [26, 27]). In Figures 8 and 9, limit cycle solutions are simulated in order to compare the 3DLM V2 and 5DLM V2 with the original 3DLM and 5DLM with $r=800$ . In all three panels for both plots, only the trajectory in red is seen, indicating that the V2 with $\mathit{FN}=1$ and the original versions produce the same solutions over the given time intervals.

Figure 8.

Time evolution of X mode solutions (a) and X vs. Y mode solutions (b), (c) of the 3DLM V2 ( $\mathit{FN}=1$ ) and the original 3DLM with $r=800$ . Panels (a) and (b) display solutions over the time interval [0, 5], while panel (c) shows the solution over the interval [2.5, 5] during which the orbit becomes regularly oscillatory in (a).

Figure 9.

Time evolution of X mode solutions (a) and X vs. Y mode solutions (b), (c) of the 5DLM V2 ( $\mathit{FN}=1$ ) and the original 5DLM with $r=800$ . Panels (a) and (b) show solutions over the time interval [0, 8], while panel (c) occurs over the interval [4, 8] during which the solution becomes a limit cycle.

In Figure 8 for the 3DLM, we plot the time evolution of the X mode in panel (a). As shown, the solution begins by an oscillation with small time-varying amplitudes and gradually increases to become regularly oscillatory at approximately $\tau =2.1$ . In the $X-Y$ phase space in panel (b), the solution spirals away from the corresponding nontrivial critical point and then becomes a closed curve, leading to a limit cycle that encloses trivial and nontrivial critical points. The closed and isolated features of the limit cycle are clearly shown in panel (c) over the time interval $\tau \in \left[\mathrm{2.5,5}\right]$ , displaying only one trajectory.

For comparison, a limit cycle solution of the 5DLM is plotted in Figure 9. Just as in plots for the 3DLM, the X mode of the 5DLM grows until it becomes oscillatory with a nearly constant amplitude after $\tau =3.2$ . Figure 9b plots solutions of the X and Y modes over $\tau \in \left[0,8\right]$ . The solution spirals outward from its initial position and gradually forms a limit cycle. Panel (c) displays the solution in the $X-Y$ space over $\tau \in \left[4,8\right]$ , showing a closed curve. The above spiral orbit and limit cycle and their transition are analyzed below using the RP.

Recurrence plots for limit cycle solutions exemplify the method’s capacity for local analysis. Figure 10a and b provides RPs for the 3DLM and 5DLM with $r=800$ . In both of these RPs, four distinct line patterns can be sequentially identified: (1) a dark block or a cluster of points (bottom left), (2) vertical and horizontal lines, (3) a large white band or a lack of points (middle), and (4) uniformly spaced diagonal lines (upper right). As shown in Figures 8 and 9, these patterns correspond to a transition from a spiral solution near a nontrivial critical point to a (nonlinear) periodic solution that encloses trivial and nontrivial critical points. Specifically, the solutions initially oscillate with very small amplitudes, progress into oscillations with larger amplitudes, and finally turn into regular oscillations. Since the patterns of the RPs directly correspond to a change in solution types, the RPs are capable of clearly identifying local dynamics and transition. By comparison, a spectral analysis cannot reveal such a local transition.

Figure 10.

A recurrence plot of a transition from a spiral orbit (the bottom left corner) to a limit cycle (the top right corner) within the 3DLM (top) and 5DLM (bottom) using $r=\mathit{800}$ .

3.5 Analysis of global reanalysis and model data

Here, we apply the RP to analyze the case with multiple AEWs. Figure 11 plots the time series for wind velocity at the 850 hPa level for the ERA-interim data (a) and GMM simulations (b). In panel (a), the amplitude of the signal is nearly constant over the entire duration. In panel (b), the signal’s amplitude generally decays first and then grows with time. In comparison to the observations from the NAMMA field campaign (Figure 1b or Figure 4a of [2]), previous studies (e.g., [2, 28]) have indicated that the GMM simulations more accurately represent actual observations as compared to the reanalysis data.

Figure 11.

Time series of the wind velocity at the 850 hPa level for the ERA-interim dataset (a) and the global mesoscale model (b).

Figure 11 displays the AEWs from the ERA-interim and GMM data, to be analyzed below. We first create two sets of idealized oscillatory solutions that mimic the “observed” and simulated AEWs. Figure 12a displays an idealized solution consisting of a periodic solution with a frequency of 1 and a periodic solution with a frequency of 2, separated by periodic solutions with negligible amplitude, similar to the ERA-interim data. In this simplified case, the oscillatory solutions both have the same amplitude. The corresponding power spectrum indicates peaks at 1 and 2 but does not readily offer information regarding the wave transitions (not shown). The RP of the signal, as seen in Figure 12b, displays uniform diagonal lines in each corner and a solid block of recurrent points in the center. These features indicate transitions between all three types of solution components. The first component with a frequency of 1 is shown by the distance between the diagonal lines. The black section in the center corresponds to no or slow variations of the states in time. For the third component (i.e., the second periodic component), uniform diagonal lines are present and are separated by a length of 0.5, consistent with the frequency of 2.

Figure 12.

Idealized oscillatory data used to mimic multiple AEWs with two different periods (a) and decaying and growing amplitudes (c). The corresponding RPs are seen in (b) and (d), respectively.

Figure 12c presents an idealized solution with the following three components: decaying oscillations, small variations with small amplitudes, and growing oscillations. The idealized data resemble the GMM data. The first and third components correspond to a spiral sink and a spiral source, respectively. Figure 12d displays the solution’s RP with horizontal and vertical lines. Such a RP pattern resembles a combined feature of the RP in Figure 7b for a spiral sink and Figure 10b for a spiral outward solution.

The RPs for the above idealized data, shown in Figure 12b and d, are used for a comparison with the RPs of the ERA-interim and GMM data, as shown in Figure 13a and b, to reveal the features of observed and simulated AEWs in Figure 11. Figure 13a presents the RP for normalized ERA-interm data at the 850 hPa level. A solid black region is present in the middle, indicating slow variations in states across this region. In areas surrounding the center, the familiar diagonal line structure associated with oscillatory behavior, and possessing some degree of periodicity, can be identified. The data results in a diagonal line structure in RPs that is comparable to Figure 12b, consisting of periodic solutions. Based on visual inspection, most of the recurrence points fall into diagonal lines, while a few stand alone. Finding an estimated period for the recurrence of AEWs through the RP is feasible. By calculating the distance between them, we obtain a good measure of the time interval between AEW occurrences. Using this approach, the distance between diagonal lines ranges from 3.5 to 5 days, consistent with the observed AEWs.

Figure 13.

Recurrence plots for the normalized wind velocity data at the 850 hPa level for the ERA-interim data (a) and global mesoscale model data (b). Panels (a) and (b) produce the RP analyses that are consistent with those in Figure 12b and d, respectively.

Figure 13b provides the RP for normalized data from the 850 hPa level generated by the GMM. The line structure present displays some similarities to Figure 12d, with a darker area in the middle and a distinct line pattern surrounding the center. Compared to Figure 13a, one major difference is the appearance of horizontal and vertical lines instead of diagonal lines. As suggested by the idealized solution in Figure 12c and d, this type of structure indicates decaying and growing oscillations, which is also similar to the combined feature in Figures 7b and 10b. The distance between the vertical and horizontal lines can be used to estimate a period for GMM data of approximately 5 days, consistent with observations as well as the RP analysis of the ERA-interim data.

4. Concluding remarks

Accurate detection of recurrent, multiple AEWs and their evolution (e.g., intensification) may provide a good indicator for determining the timing and location of tropical cyclone formation and initial intensification. Since large-scale AEWs have better predictability, such an indicator may help extend the lead time of TC formation prediction. To achieve this goal, in this study, we first deployed the recurrence plot method in Python and verified our implementation using several types of idealized solutions and three types of solutions obtained in the classical 3DLM and 5DLM. We then applied the RP method to analyze two, real-world datasets, ERA-interim data and global mesoscale model data, in order to reveal the time-varying amplitudes and frequencies of multiple AEWs over a 30-day period between late August and September 2006. Compared to traditional global analysis methods (e.g., spectral analysis), the RP method is effective in showing temporal growing or decaying oscillations and the transition between these two types of solutions. As a result, the RP analysis of global data not only produces a good estimate of the period for AEWs but also displays the weakening or intensifying trend of AEW intensities, as shown in Figure 13.

A summary on the performance of the RP in analyzing various types of solutions is provided below. We first performed the analysis for (1) four basic types of data, including periodic, quasiperiodic, and chaotic solutions as well as Gaussian white noise, (2) idealized decaying/growing oscillations, (3) steady-state solutions using the 3DLM, (4) limit cycle solutions using the 3DLM and 5DLM with a focus on the transition from a spiral solution near a nontrivial critical point to a nonlinear periodic solution that encloses trivial and nontrivial critical points, and (5) idealized solutions that mimic the features of selected AEW data with different periods and different amplitudes. For basic types of solutions, the RP produces a consistent analysis as compared to previous studies. For time-varying oscillations, we discussed how recurrence may lead to horizontal and vertical lines in RPs. The distance between two consecutive horizontal (or vertical) lines gives an estimate of half of a period (i.e., $T/2$ ). For the fourth and fifth types of datasets, the corresponding RP analyses clearly display local transitions between solution components with various periods or amplitudes.

Based on its promising performance in revealing the features of idealized solutions, the selected RP method was applied to analyze ERA-interim data and GMM data, yielding a RP analysis consistent with the analysis of the idealized AEWs from the fifth dataset. While the ERA-interim data resemble the idealized data in Figure 12a with comparable amplitudes but different periods, the GMM data resemble the idealized data in Figure 12c with decaying and growing oscillations. The RP analysis of the ERA-interim data in Figure 13a produces an estimated period of 3.5–5 days, and a spectral analysis of the same data yields a dominant period of 4.7 days. Both estimated periods are consistent with the observed results. However, only the RP analysis can reveal local variations or transitions. By comparison, the RP pattern of the global model data in Figure 13b that displays horizontal and vertical lines is similar to the RP plot of the idealized data in Figure 12d. While the RP produces a comparable period of approximately 5 days, it additionally reveals realistic information regarding the weakening and intensifying trend of AEWs.

Recent studies within simplified and generalized Lorenz models [22, 29, 30, 31, 32] suggest the importance of detecting oscillatory components for extending the lead time of weather prediction. The analysis with RP is encouraging but largely produces qualitative information. Our future work includes an application of recurrence quantification analysis (RQA) methods (e.g., [33, 34, 35]) for quantitatively determining the “recurrence rate” and “determinism” in order to quantitatively measure the recurrence and determinism (or “predictability”) of the recurrent solutions.

Acknowledgments

We thank Dr. Y.-L. Wu, Dr. J. Chern, and Ms. S. Faghih-Naini for their scientific discussions and technical help. Thanks are also extended to Dr. B. Bailey, Dr. F. De Sales, and anonymous reviewers for their comments on the original manuscript. We are grateful for the support from the College of Science of San Diego State University.

A.

The perturbation method was applied to derive the so-called V2 system in order to analyze solutions of the 3DLM, as well as the 5DLM, with initial conditions near the nontrivial critical point. The 5DLM V2 and its non-dissipative V2 were previously documented in [13] and [7], respectively. This appendix simply discusses the 3DLM V2. For the 3DLM, in Eqs. (1)–(3) without the feedback term ${\mathit{XY}}_1$ , the total field is decomposed into the basic state ( $X_c,Y_c,Z_c$ ) and the perturbation ( $X^{\prime },Y^{\prime },Z^{\prime }$ ) with respect to the basic state in the form of $X=X_c+X^{\prime }$ for each state variable, yielding a new set of equations:

$$\frac{{\mathit{dX}}_c}{\mathit{d\tau }}+\frac{d{X}^{\prime }}{\mathit{d\tau }}=-\sigma \left(X_c+X^{\prime }\right)+\sigma \left(Y_c+Y^{\prime }\right),$$ EA1

$$\frac{{\mathit{dY}}_c}{\mathit{d\tau }}+\frac{d{Y}^{\prime }}{\mathit{d\tau }}=-\left(X_c{Z}_c+X_c{Z}^{\prime }+Z_c{X}^{\prime }+X^{\prime }{Z}^{\prime }\right)+r\left(X_c+X^{\prime }\right)-\left(Y_c+Y^{\prime }\right),$$ EA2

$$\frac{{\mathit{dZ}}_c}{\mathit{d\tau }}+\frac{d{Z}^{\prime }}{\mathit{d\tau }}=\left(X_c{Y}_c+X_c{Y}^{\prime }+Y_c{X}^{\prime }+X^{\prime }{Y}^{\prime }\right)-b\left(Z_c+Z^{\prime }\right).$$ EA3

Note that the above basic state solutions are determined using the time-independent, nontrivial critical point solutions, defined as (e.g., [11])

$$Z_c=r-1,$$ EA4

$$X_c=Y_c=\pm \sqrt{{\mathit{bZ}}_c}.$$ EA5

Since the basic state solutions are time independent, Eqs. (A1)–(A3) are reduced to become

$$\frac{d{X}^{\prime }}{\mathit{d\tau }}=-\sigma \left(X^{\prime }-Y^{\prime }\right),$$ EA6

$$\frac{d{Y}^{\prime }}{\mathit{d\tau }}=\left(r-Z_c\right)X^{\prime }-X_c{Z}^{\prime }-Y^{\prime }-\mathbf{FN}\left(X^{\prime }{Z}^{\prime }\right),$$ EA7

$$\frac{d{Z}^{\prime }}{\mathit{d\tau }}=Y_c{X}^{\prime }+X_c{Y}^{\prime }-b{Z}^{\prime }+\mathbf{FN}\left(X^{\prime }{Y}^{\prime }\right),$$ EA8

which describe the time evolution of the perturbations. $\mathit{FN}$ is a flag identifying whether the system is fully nonlinear ( $\mathit{FN}=1$ ) or not ( $\mathit{FN}=0$ ). The above system with $\mathit{FN}=0$ or $\mathit{FN}=1$ is referred to as a Version 2 system, denoted by V2. In the V2 system, there are no nonlinear terms that involve the product of two basic state variables. The system only allows interactions between one basic state variable and one perturbation (i.e., $X_c{Y}^{\prime }$ ) or between two perturbation variables (e.g., $X^{\prime }{Y}^{\prime }$ ). By setting $\mathit{FN}=0$ , the system is linear with respect to the perturbation. When the initial perturbations are small, nonlinear terms with two perturbation variables are smaller. Thus, such a linear system is capable of depicting the initial time evolution of the full 3DLM when the perturbation is small. In comparison to the full 3DLM, the nonlinear V2 model ( $\mathit{FN}=1$ ) produces the same solutions, even at the onset of chaotic irregular oscillations, as shown in Figure A1. The time evolutions of the $Y^{\prime }$ solutions for both $\mathit{FN}=0$ and $\mathit{FN}=1$ are shown in Figure A2. The solutions are initially the same and begin to diverge as time further progresses, indicating the impact of nonlinearity. An eigenvalue analysis of the 3DLM V2 is provided in Section 2.1.3 of Reyes [24].

Figure A1.

Time evolution of the X mode solutions of the original 3DLM and 3DLM V2 with $\mathit{FN}=1$ over the time interval $\mathit{0}\le \mathit{\tau }\le \mathit{50}$ with a time step $\Delta \mathit{\tau }=\mathit{0.001}$ .

Figure A2.

Time evolution of $Y^{\prime }$ linear (blue) and nonlinear (red) solutions using the 3DLM V2 with $\mathit{0}\le \mathit{\tau }\le \mathit{110}$ and $\mathrm{\Delta }\mathit{\tau }=\mathit{0.001}$ .