Searching for Chaos in Tropical Cyclone Intensity: A Machine Learning Approach

2.1. Maximum intensity equilibrium

If the PI equilibrium is the only possible state of TC development whose statistical properties are stationary for intensity climatology, how can we use this equilibrium to examine TC intensity chaos? In this study, we will follow the same approach as in Kieu et al. (2022) and assume that the PI equilibrium can be characterized by a low-dimensional phase space where chaos manifests. Note that Kieu et al.’s approach based on the phase-space reconstruction method to directly search for the dimension of a chaotic attractor contains significant subjectivity and is sensitive to data noise (e.g., Kantz and Schreiber 2003). Here, we propose to use ML to quantify the TIP range, which can also confirm the existence of TC intensity chaos, albeit less directly as compared to the phase-space reconstruction method.

A key part in searching for chaos at the equilibrium is therefore to obtain first a statistically stable PI state so that one can analyze it. Given such a state, one can then extract the time series of key TC scales such as U_max, V_max, W_max, or P_min, which can serve as input features for ML training. Details of a long simulation of TC intensity that ensures such a stable PI state and how to extract the required data for the ML approach will be given in Section 3, which are identical to those used in Kieu et al. (2022).

2.2. Machine learning approach

Broadly speaking, machine learning (ML) can be considered as a framework that can search for rules from data. Given an ML architecture, a measure of accuracy, and input data, the rules can be obtained within a prescribed level of accuracy. The key advantage of ML in practical applications lies in its ability to learn rules from input data without a priori knowledge, provided that the input data is sufficiently good (i.e., the input data can ensure several criteria including i) comprehensiveness, ii) relevancy, ii) consistency, and iv) uniformity.) With an inherently large volume of data, climate and weather prediction provide a great domain for ML applications, which justifies the surge of ML applications in atmospheric science recently.

Specifically for TC intensity, ML offers a unique way to study low-dimensional chaos. To set up a context for applying ML to our TC intensity chaos problem, we will focus hereafter on supervised ML, which requires a set of input data and corresponding targets (labels) for training an ML model. At a basic level, supervised ML models need a surjective mapping between an input training dataset $(\mathcal{T})$ and a target dataset $\left(\mathcal{L}\right)$ (i.e., one $y\in \mathcal{L}$ will have at least one $x\in \mathcal{T}$) so that the training can be carried out. For a typical time-prediction problem (i.e., given a state of a system at one time t = 0, one needs to predict the state of the system at a later time $t=\tau$), this mapping can be considered as a propagator from a given initial condition to the later time $\tau$. Mathematically, such a propagator can be expressed as $x(\tau )=M(\tau )x(0)$, where $M(\tau )$ is the propagator from t = 0 to $\tau$ and x(t) is the model state at time t.

For a full-physics model, $M(\tau )$ is nothing but a numerical model with governing equations integrated from t = 0 to $t=\tau$. For ML, $M(\tau )$ is a however nonlinear operator that is learned from a training dataset. In principle, the more data we have, the better an ML model can search for underlying rules and build $M(\tau )$ without any physical equations. Thus, we can feed an ML model with a large amount of data, and let it figure out the best possible relationship between t = 0 and $t=\tau$. For deep learning that is based on neural networks, an ML model with sufficient layers and depth should capture a nonlinear mapping between two time slices, making it suitable for TC intensity prediction.

From this perspective, it is immediate that chaotic systems will pose a challenge to any ML model, because one input may give totally different outcomes after reaching predictability limit T (i.e., one $x\in \mathcal{T}$ would give two different $y_1,\hspace{0.17em}\hspace{0.17em}{y}_2\in \mathcal{L},\hspace{0.17em}\hspace{0.17em}\forall \tau >T$, where two input x₁, x₂ are practically considered to be the same as x if their difference is within some measurement errors, $|x_i-x|\hspace{0.17em}\hspace{0.17em}\le \hspace{0.17em}\hspace{0.17em}ϵ,i=1,\hspace{0.17em}\hspace{0.17em}2$ for a sufficiently small uncertainty ϵ). So, there exists no longer a good mapping between the training and the label datasets, and ML models cannot learn any rule from data.

The deterioration of ML models after entering the chaotic regime as described above suggests, however, a unique way to study chaotic systems. Specifically, we will search for a lead time T beyond which an ML model can no longer be trained from any input dataset, which gives us a direct estimation of the predictability range for a chaotic system. This approach is natural in the sense that an ML model should generally be able to predict the next state of a system from a given input, if the system remains predictable and sufficient training data is provided. As soon as the system enters a chaotic regime, ML models are no longer trainable. From this perspective, ML models are naturally a great tool for studying chaos.

Our aforementioned use of ML to examine predictability is well suited for TC intensity, as this approach serves two purposes: i) it verifies if a low-dimensional representation is sufficient for TC intensity, and ii) it helps estimate the TIP range in that low-dimensional phase space. Our underlying hypothesis is that an ML model can predict TC intensity in a low-dimensional phase space, whose dimensions correspond to several basic TC scales, up to a certain lead time T. Beyond this lead time T, ML models can no longer be trained to predict TC intensity, thus revealing TC low-dimensional chaos and providing us an estimation for TIP.

The results in Kieu et al. (2022) provide a pathway to verify this hypothesis with ML. Specifically, we will assume that TC intensity can be described by four dimensions corresponding to four TC scales including V_max, U_max, W_max, and P_min. While it is not known in advance the exact dimension of the PI attractor, Kieu et al. (2022) suggested that a minimum dimension of 4 should be sufficient to capture TC intensity chaos within the deterministic framework. As such, we will treat these four dimensions as input features for several ML models to be presented in Section 3b. With these ML models, we can examine how they forecast intensity at different lead times and estimate the TIP range as expected.

3.1. CM1 model configuration

In this study, the same axisymmetric configuration of the Cloud Model (CM1) (Bryan and Fritsch (2002) was used as in Kieu et al.’s study (2022), which produces a quasi-stationary PI state during a long integration of 100 days. This model configuration has 360 grid points on a stretching grid in the radial direction, with the highest resolution of 2 km in the vortex’s inner core region and stretched to 6 km in the outer core region. Unlike the radial direction, the model was configured with 61 levels in the vertical direction, with a fixed resolution of 0.5 km. This fixed vertical resolution was found to be more numerically stable and also less restrictive when choosing the number of vertical levels. In addition, we applied the open-radiative lateral boundary conditions option to the radial direction, and free (no) slip boundary to the top (bottom) boundary in our simulation. The model was initialized from the tropical Jordan sounding on an f-plane, with fixed sea surface temperature (SST) = 302.15 K.

Similar to the results in Hakim (2011) and Kieu and Moon (2016), a long simulation of a quasi-stationary TC intensity with the CM1 model would require a proper choice of model physics to avoid the gradual change in TC environment inside a box domain, which can cause decaying due to the transport of low angular momentum in the outer core region. A simple treatment for this environment change is to apply a fixed Newtonian cooling relaxation of 2 K day-1, as in Kieu et al. (2022), which can result in a quasi-stationary maximum intensity equilibrium in the CM1 simulation for 100 days. Along with this radiative forcing, a suite of other physical parameterizations were also used, which include the YSU boundary layer scheme, the TKE subgrid turbulence scheme, and the explicit moisture Kessler scheme with no cumulus parameterization.

With the PI equilibrium in our 100-day simulation established, the model was then output at every time step of 36 seconds, producing a dataset of length $\approx \hspace{0.17em}\hspace{0.17em}\mathcal{O}({10}^6)$. All time series of four major TC scales including V_max, U_max, W_max, and P_min were extracted from the model output and further split into three subsets for ML training—the training, validation, and test sets with a ratio of 90%, 5%, and 5%, respectively. To ensure that all data are selected at the PI equilibrium, the first 10 days of simulations were discarded. Other details of this CM1 100-day simulation can be found in Kieu et al. (2022), and so we do not repeat them here.

3.2. Deep-learning models

Given the low dimensionality of feature vectors used for our ML training, we present in this study several deep-learning models for TC intensity prediction. Specifically, three popular ML architectures including a deep neural network (DNN) model, a gated recurrent unit (GRU) model, and a long-short term memory (LSTM) were implemented. The applications of these deep learning models in the weather domain have been rapidly growing due to their capability as well as the availability of computational resources, which help accelerate their execution for practical problems. With four TC scales as input features and one real-value output V_max representing TC intensity, predicting TC intensity with the above ML models thus becomes a familiar supervised regression problem for which these ML models fit very well.

For our TC intensity prediction with DNN, a simple design of 3 hidden layers with layer sizes of 32, 64, and 64 was used, followed by an output layer of size 1 that corresponds to V_max. Each neural layer was applied a standard ReLU activation, which helps ML models capture nonlinear effects and increase the interaction among layers. One could certainly design a deeper neural network for a more complex relationship between input and output layers. However, our experiments with different DNN designs showed very little improvement for more than 3 hidden layers when predicting TC intensity in a low-dimensional input space. As such, a fixed design of 32, 64, and 64 nodes was used.

For LSTM and GRU, these are recurrent neural models that further require a data interval in the past to capture the memory in the training data. Our model architectures for these LSTM and GRU models thus need some additional setup. Specifically for these recurrent network models, we used a range of time slices, i.e., $t_i,\hspace{0.17em}\hspace{0.17em}i\in \hspace{0.17em}\hspace{0.17em}[-M,\hspace{0.17em}\hspace{0.17em}0]$, as input for LSTM/GRU models when predicting TC intensity at any given lead time. Here, M determines the number of time slices in the past that are needed for recurrent networks, which varies from 5–20 in this study. To avoid overfitting during the training process, we also used three layers of size 16, 32, and 64, with a dropout rate of 0.5. Technically, dropout is a type of regularization that can help reduce overfitting in ML models. There is no particular formula to choose the value for this hyperparameter, other than empirical trials. For our intensity chaos problem, this dropout turns out to be important to ensure good model performance.

All of these ML models employed the mean absolute error (MAE) metrics for the accuracy and the root mean squared errors for the loss function, with a fixed number of training epochs set to be 200. The standard optimizer for the gradient search based on the stochastic mini-batch learning method, the so-called Root Mean Squared Propagation (RMSprop), was applied to all training. Because of the different scales of the wind and pressure variables, all input data was normalized by the standard deviation around their mean value, which corresponds to the PI state of the model vortex at the quasi-stationary equilibrium.

4.1. Predictability limit on attractor

Given the current definition of TC intensity in terms of the maximum 10-meter wind, we examine first TC intensity chaos using V_max as a metric for TC intensity in our ML models. Recall here that, if TC dynamics possesses chaos at the PI equilibrium, then TC intensity should have a limited predictability range that dictates the value of TIP. Thus, a finite TIP range can serve as a proxy for the existence of TC intensity chaos that we are searching for.

In this regard, Figure 1 shows the training absolute mean error as a function of epochs (iterations) for three models at three forecast lead times including 3 minutes, 1 hour, and 3 hours. One notices indeed that the training errors rapidly decrease for $\tau =3$ minutes in all three ML models (black curves in Figure 1), reaching a relative minimum error of ≈ 0.05, 0.12, and 0.22 for LSTM, GRU, and DNN models, respectively. Looking at the correlation between the ML forecast and the true TC intensity for the test data in Figure 2 (red points), it confirms that all ML models could predict very well TC intensity variability for the short lead time $\tau =3$ minutes in the 4-dimensional phase space. This result is noteworthy, because these ML models require a minimum number of input features, yet they could produce a good forecast of TC intensity based solely on training data. From this perspective, Figures 1 and 2a help confirm that a low-dimensional phase space suffices to predict TC intensity variability at least for a short lead time, without any physical or governing equations.

Figure 1

Figure 2

At the 1-hour lead time, Figure 1 (red curves) shows, however, that all three ML models start losing their ability to be trained quickly. By 3 hours (blue curves), all ML models can no longer be trained, with their errors roughly the same ≈75–85% relative to the initial error value during the entire training period no matter how many epochs are used. Their predictions for the test set at the 3-hour lead time display almost completely no correlation to the true intensity (Figure 2, blue dots). This result reveals that $\tau =3$ hrs is the longest lead time that these ML models can predict TC intensity at the PI limit. It is of interest is that this estimation is also consistent with the estimation from attractor invariants based on a leading Lyapunov exponent and the Sugihara-May correlation in Kieu et al. (2022), which showed that TC intensity loses predictability in just ≈3–6 hours as soon as TCs reach their PI equilibrium.

The dependence of these ML-based intensity predictions on forecast lead times is best seen when we compare these predictions to a reference (or climatology) forecast, which is taken to be a simple average of V_max at the PI equilibrium. Figure 3a shows the forecast skill of three ML models relative to this reference forecast as a function of lead times. It is apparent from Figure 3a that ML models perform best for $\tau <3$ hours. Beyond this, the ML-based prediction skill is no better than a simple forecast using just the averaged V_max at the PI equilibrium. This short predictability is further supported by the error growth curve (Figure 3b), which displays a typical behavior of chaotic systems with rapid error growth during the initial period and reaching a saturation level after ≈ 3 hours.

Figure 3

We emphasize that the decaying of the ML-based forecast skill with $\tau$ does not hold true for any system. In fact, a simple experiment using purely random data as an input for ML training would result in zero forecast skill at all lead times (not shown). This is because the information from one time step does not have any influence on the next, and so ML models cannot learn anything. On the other hand, for completely periodic systems, the forecast skill will always take a constant value of 1 for all lead time $\tau$ as discussed in, e.g., Sugihara and May (1990). As such, the decaying curve of the forecast skill or the error growth curve shown in Figure 3a–3b is an inherent characteristic of chaotic systems, which is captured by our ML models.

Given the average V_max = 84 ms^–1 with a standard derivation is ≈7.5 ms^–1 at the PI equilibrium, the TIP range obtained from Figure 3a implies further that TC intensity must vary indistinguishably within an interval of 84 ± 7.5 ms^–1 after just 3 hours, even for a perfect TC model. This TIP range may even be shortened if stochastic forcings, asymmetric processes, or model internal errors are taken into account as discussed in Nguyen et al. (2020) or Kieu et al. (2022), which are, however, beyond the scope of our study here. Despite these issues, the results obtained from our ML models herein can at least advocate the existence of intrinsic variability of TC intensity due to chaotic dynamics, which prevents the absolute intensity errors in any TC models from being reduced to zero.

Among the three ML models, it should be mentioned that LSTM and GRU appear to perform best in terms of either the training/validation mean errors due to their use of extra past information (Figure 1–2) or forecast skill (Figure 3a). This past information contains some temporal relationships that help improve future prediction. Thus, LST/GRU presents a very different way of forecasting as compared to the traditional approach based on physical principles. To some extent, recurrent networks improve their prediction in the same way that four-dimensional data assimilation optimizes an initial condition over an interval instead of just a single time slice. Despite this extra information from the past, the predictability of V_max could not be lengthened beyond 3 hours in both LSTM and GRU models as shown in Figures 2, 3.

4.2. Metric dependence

Because predictability is metric-dependent, an apparent question is how the estimation of TIP changes when using P_min for TC intensity instead of V_max. In this regard, Figure 3c–d shows the ML-based forecast skill and error growth curve for P_min as a function of lead time for three ML models, similar to Figure 3a–b. It is of significance to observe that all models could again capture similar decaying of the P_min forecast skill, but with a significantly slower rate and thus a longer TIP range of ≈8–9 hours as compared to 2–3 hours for V_max. Such slower decay of the P_min forecast skill in Figure 3c well accords with slower error growth and a longer time to reach the error saturation as shown in Figure 3d, similar to what was obtained from the phase-space reconstruction method in Kieu et al. (2022).

The fact that these ML models capture such a different predictability range between V_max and P_min is intriguing. Recall that V_max and P_min are highly correlated in terms of temporal variability due to their pressure-wind relationship. However, P_min represents the total mass at the storm center while V_max fluctuates more vigorously due to fine-scale processes at each model grid point. As such, P_min tends to better display a slow component of TC dynamics, which ML models could somehow detect even when training data contains strong fluctuations from the wind field. From this standpoint, using P_min for TC intensity could lengthen the range of intensity predictability for operational forecasts as previously noticed (e.g. Kieu et al., 2022; Klotzbach et al. 2020; Magnusson et al., 2019).

Regardless of intensity metrics V_max or P_min, the above results reiterate the finite range of TIP as obtained from different ML models. Such a finite range is in fact held for all TC scales that we have examined, not just P_min or V_max. Hence, the existence of low-dimensional chaos of TC dynamics as captured in our ML models is well supported, even though we do not currently know any mathematical model that can describe this chaotic dynamics of TC intensity in a low-dimension phase space. In particular, the existence of TC intensity chaos derived from this TIP information confirms that a part of intensity variability in TC models must be inherent to TC dynamics, which cannot be removed from model outputs simply by improving the model physics or initial conditions.

4.3. Sensitivity experiments

To further address the robustness of our results, this section presents some additional experiments in which more input features, different recurrent windows M, and a coarser sampling frequency (30 minutes) for input training data are used. Note that these sensitivity experiments cover just a small part of possible sensitivities that one can examine. For example, one can design deep neural networks with an arbitrary number of layers, nodes, dropout rates, or data augmentation. Within the scope of this study, we will however limit our sensitivity experiments to several key experiments, which suffice to highlight the important points that we want to present herein.

Figure 4a–4b show the forecast skill and the error growth of V_max for the first group of sensitivity experiments for which the recurrent timesteps M = 5, 10, 20 at a sampling frequency of 30 minutes. Recall from Section 3.2 that M determines how many time slices in the past up to a present time t will be used to predict the future state at $t+\tau$. Thus, a larger M (i.e., a longer window of past information) would allow for more input information, and should increase the forecast skill of ML models.

Figure 4

As seen in Figure 4a–4b, forecast skill in these sensitivity experiments all decays quickly during the first 3 hours across the models and recurrent window M, similar to the control settings in Figure 3a–3b. For each model, note however that the longest window M = 20 (dotted line in Figure 4a) appears to be slightly more skillful as compared to the shorter ones, which can be seen in all models. This indicates that longer input windows for recurrent models appear to help improve a forecast skill for chaotic systems, albeit marginally.

The most significant change in these sensitivity experiments is, however, a small but persistent skill of all ML models relative to the reference forecast, even at long forecast lead times $\tau >3$ hours (Figure 4a). Such a marginal skill is due to the coarser sampling frequency of input data, which eliminates fast fluctuation in the training data and results in some skill relative to the constant reference forecast. From the error growth perspective, it is seen however that there is no significant change in the time needed for the model forecast error to reach its saturation (Figure 4b). All ML models display the same error growth saturation after ≈3 hours, which is an intrinsic property of a chaotic system and independent of how reference forecasts are defined.

Our sensitivity experiments with more input features capture a similar behavior of the error growth and forecast skill as in the recurrent time window experiments (Figure 4c–4d). Specifically, adding more features such as RMW, warm core anomaly (T’), or both does not improve the TIP range. That is, both forecast skill and error growth show rapid decay/increase during the first 3 hours and then approach a saturation limit afterward. This result supports the existence of low-dimensional chaos for TC intensity in the sense that more dimensions in the phase space will not help improve any estimation of attractor invariants as discussed in Kantz and Schreiber (2003) or Kieu et al. (2022). The consistent TIP estimation among ML models and sensitivity experiments as shown above, therefore, confirms the robustness of our finding about the existence of low-dimensional chaos for TC intensity.