For this research, we used the best track dataset for TCs in WNP with wind speeds of 34 kt or higher. This data was provided by the Joint Typhoon Warning Center (JTWC) and spans from 2004 to 2021. Oceanic variables, specifically SST and DAT, were computed using analysis/reanalysis data from the Hybrid Coordinate Ocean Model and the Navy Coupled Ocean Data Assimilation nowcast/forecast system (HYCOM+NCODA), as provided by the Naval Research Laboratory. DAT values were calculated at varying depths ranging from 10 m to 120 m, at 10 m intervals (DAT₁₀ through DAT₁₂₀). These values were used to compute various oceanic components such as MPI (MPI₁₀ to MPI₁₂₀, henceforth referred to as MPIs), POT (POT₁₀ to POT₁₂₀, henceforth referred to as POTs), OC_PI (OC_PI₁₀ to OC_PI₁₂₀, henceforth referred to as OC_PIs), and NGR (NGR₁₀ to NGR₁₂₀, henceforth referred to as NGRs). Atmospheric variables were obtained from the Global Forecasting System analysis, provided by the National Centers for Environmental Prediction (NCEP), with a spatial resolution of 1° x 1° and a temporal resolution of 6 hours. The average radius of gale-force winds in WNP is approximately 200 km (Kim et al., 2022). Wang and Toumi (2021) have identified that the effective radius for TC-induced sea surface cooling is roughly of the same magnitude. To accurately capture the effects caused by a TC, we averaged both oceanic and atmospheric variables within a 200 km radius of the TC center. Furthermore, to isolate and remove the influence of the TC from our data, we analyzed conditions from three days before the storm’s arrival.

NGR is calculated as the difference between the rate of energy generation (G) and the rate of surface frictional dissipation (D), all within the context of Emanuel’s MPI framework. Kim et al. (2022) showed that C_d is the most critical factor in determining the magnitude of NGR. This C_d not only significantly influences D but also plays an important role in vertical mixing within the ocean, which in turn affects the saturation enthalpy determined by SST. Therefore, using a more realistic C_d is crucial for comprehending the mechanisms of TC intensification. In this study, C_d parameterization from Soloviev et al. (2014), based on two-phase parameterization and observations from previous studies, was used.

The NGR is computed using Emanuel’s software package, with a modification: the SST in the original equation is replaced by DAT. Additionally, the model employs a wind-speed-dependent drag coefficient (C_d (V)) rather than using a constant drag coefficient.

The equations are as follows:

where DAT is the depth-averaged temperature (Price, 2009), T _o represents the TC outflow temperature, C_k is the enthalpy coefficient, ρ is the air density, k_o ^* is the sea surface saturation enthalpy, k is the surface enthalpy in the TC environment, C_d is the drag coefficient, and V is the surface wind speed.

Higher NGR values suggest that more energy is available for TC intensification. Given its superior performance in predicting short-term TC intensification, NGR can be used as an ideal predictor for RI events. Its ability to more accurately capture the ocean’s contribution to TC intensity changes within a 24-hour range makes it especially suitable for the RI events prediction.

The TC-induced vertical mixing depth is determined by various parameters such as the size, intensity, and translation speed of TC, the Coriolis effect, and the vertical structure of the upper ocean. The depth of this mixing is crucial because it determines the SST where heat exchange occurs during the intensification of TC. Lin et al. (2013) showed that using an average mixed layer depth of 80 m minimizes the bias in the MPI for TCs that are the Saffir-Simpson scale Category 2 or higher. Price (2009) indicated that 0-100 m DAT can adequately represent the mixing caused by major TCs. Meanwhile, LEE19 conducted a sensitivity analysis using the NGR for various depths of mixing and showed that fixing the depth at 50 m yielded the highest predictive performance for changes in the intensity of the overall TCs. In this study, we took a comprehensive approach to account for the sensitivity of vertical mixing depth and to explore all possible combinations of predictors. We calculated all major predictors, including POT, NGR, and OC_PI, based on DATs. These calculations were done at 10-meter intervals up to a depth of 120 meters and were subsequently included in our predictors pool ( Table 1 ).

Besides NGR, this study incorporates other well-established predictors commonly employed in statistical-dynamical models for TC intensity forecasting by various organizations (DeMaria et al., 2005; Knaff et al., 2005; Kim et al., 2018). In this study, we utilized a total of 65 potential predictors, encompassing a diverse range of factors. These include 5 static predictors, 7 atmospheric synoptic predictors, SST, MPI and the POT derived from it, OHC, and 49 predictors based on the DAT. All these predictors are summarized in Table 1 . To evaluate the impact of NGR on the ML-based prediction of RI events, our study uses two distinct sets of these predictors. The first set consists of commonly used predictors related to TC intensity change, as identified in numerous statistical models (in Table 1 , excluding NGRs). The second set incorporates NGR into the first set in Table 1 (as illustrated in Figure 1 ).

In this study, we employ a diverse ensemble of well-known classifiers to predict RI events in the WNP. The ensemble includes a Decision Tree (DT), Logistic Regression (LR), SVM, k-nearest Neighbors (KNN), and Feed Forward Neural Network (FNN). DT serves as a comprehensive data-mining tool, which is adept at generating decision-making rules, identifying patterns, and uncovering knowledge embedded in archived databases (Quinlan, 1987). More specifically, the DT algorithm evaluates the conduciveness of environmental conditions for RI by systematically checking whether specific environmental predictors satisfy the thresholds set by the trained tree model. LR is used to predict a categorical variable such as the class label (Walker and Duncan, 1967). It is an extension of linear regression, where the classification problem is converted into a regression problem by estimating the log (odds) of each class in place of probability itself. The model uses the logistic function to squash the output of a linear equation between 0 and 1, making it interpretable as a probability. This method is prized for its simplicity, interpretability, and effectiveness in various domains. The SVM is designed to discover a hyperplane that best separates the data classes (Cortes and Vapnik, 1995). It achieves this by maximizing the margin between different classes in the feature space. The KNN algorithm makes predictions by storing all training data and identifying the classes of the k closest neighbors to each test sample (Keller et al., 1985). It aims to classify an unknown sample based on the known classifications of its nearest neighbors. Finally, FNN is a straightforward artificial neural network composed of an input layer and an output layer. The flow of crucial input information in FNN moves strictly from its input layer to its output layer, making the model especially well-suited for parameter identification tasks (Zhang et al., 2022). Each of these classifiers brings its own set of strengths to the ensemble, contributing to a more robust and reliable RI prediction for the WNP.

To enhance predictability, we employ a hard-voting ensemble method that combines the predictions of the individual classifiers. In this approach, each classifier ‘votes’ for a class when presented with a test instance. The ensemble then selects the class that receives the majority of votes as its final prediction. By employing this hard-voting scheme, we aim to benefit from the complementary strengths of each classifier, thereby achieving a more robust and accurate model for predicting RI events in the WNP. Given that RI events are not commonly observed, as shown in Table 2 , there is a clear class imbalance in our dataset. To effectively address and rectify this imbalance, we made use of the Synthetic Minority Over-sampling Technique, commonly known as SMOTE (Chawla et al., 2011). This method effectively tackles the issue of class imbalance in datasets, which is critical in many ML applications. Generating synthetic data for the minority class and creating new data points between existing ones helps balance the dataset. This balance is crucial for training unbiased models and ensuring they effectively learn the characteristics of all classes. This approach ultimately leads to an enhancement in the overall accuracy and performance of the model, making it more reliable for real-world applications. In this study, Principal Component Analysis (PCA) was applied to the pool of predictors to combat multicollinearity within the model. The integration of PCA into our ML-based classification model brought several advantages. It effectively streamlined the dataset by reducing dimensionality, which helped mitigate issues related to the curse of dimensionality and overfitting. This feature reduction also led to improved computational efficiency. Additionally, by focusing on the primary directions of data variance, PCA successfully filtered out noise and irrelevant information, resulting in a more refined dataset (Tefas and Pitas, 2016). Notably, we only used those principal components that represented at least 99% of the cumulative explained variance as predictors in our model. Given the limited size of our dataset, we employed a 10-fold cross-validation approach, ensuring the selection of the most effective models and preventing overfitting during training. The dataset from 2004 to 2018 was designated for training, while data from 2019 to 2021 was reserved for testing. Within the training data, there were 627 RI cases and 3,388 non-RI cases. Meanwhile, the testing dataset comprised 103 RI cases and 581 non-RI cases, as detailed in Table 2 .

In the realm of binary classification models, the evaluation of predictor significance is pivotal for model accuracy and interpretability. Among various statistical measures, Cohen’s d is an effective tool for quantifying the discriminative power of predictors. Originally designed to measure the standardized difference between two means in psychological research, Cohen’s d can be adapted to assess how individual predictors differentiate between the two classes of the model, typically labeled positive and negative. By calculating the difference in means of a predictor for each class and dividing it by the pooled standard deviation, Cohen’s d provides a standardized effect size, facilitating direct and quantifiable comparison of the predictor’s impact across different models and datasets. Cohen’s d calculated as:

where M₁ and M₂ are the means of the predictor values for each of the two classes. SD_pooled is the pooled standard deviation of the predictor values across both classed. It is computed as:

where SD₁ and SD₂ are the standard deviations for each class, and n₁ and n₂ are the sample sizes. A higher Cohen’s d value indicates greater separation between the classes based on the predictor, signifying its importance in the classification task. Typically, Cohen’s d values around 0.2 are considered small, around 0.5 medium, and 0.8 or higher, large. This gradation helps in pinpointing the predictors with the most significant roles in distinguishing between classes.

In binary forecasts where models predict an RI event or non-RI event for a training and test set, evaluation metrics comprise elements from a confusion matrix that compare observations to model forecasts ( Table 3 ). True Positive (TP) is the number of correct forecasts of RI events, whereas False Positive (FP) is the number of incorrect forecasts. False Negative (FN) is the number where the model did not forecast RI but, RI was observed. True Negative (TN) is the number where the model did not forecast RI and RI was not observed.

Accuracy (ACC) is used to measure the overall performance of a binary classifier and is measured as

FAR is the number of incorrect forecasts of RI divided by the total number of RI forecasts. FAR is calculated as

POD) is the ratio of the correct forecasts of RI occurrences to the actual number of RI occurrences and is calculated as

Precision measures the accuracy of positive predictions in classification problems. It’s the ratio of the correct forecasts of RI occurrences to the total number of positive predictions (which includes both TP and FP). Precision is calculated as

Peirce skill score (PSS), also known as the Hanssen-Kuipers skill score measures skill relative to an unbiased random reference forecast and is calculated as

The F-1 score is a way of combining the precision and POD of the model, and it is defined as the harmonic mean of the model’s precision and POD. The F-1 score is calculated as

A perfect forecast model would achieve an ACC, POD, and PSS score of 1 and a FAR score of 0. In general, higher values of ACC, POD, Precision, PSS and F-1 score, coupled with a lower FAR, indicate superior model performance.

In this section, we examine the classification performance of potential predictors for RI events before developing a classification model. The mean distribution for RI and non-RI classes for each predictor, the effect size of the mean differences between these classes, and the correlation coefficients with the 24-hour intensity change were analyzed ( Table 4 ; Figures 2 , 3 ). Excluding DAT-based predictors, the ocean temperature and MPI theory-based predictors (SST, MPI, POT, OHC) exhibited the highest Cohen’s d and correlation coefficients. Following these, static predictors such as DVMX and LAT displayed the next highest values of Cohen’s d. Synoptic predictors, apart from wind shear-related predictors (SH200, SH500, U200), generally demonstrated lower predictive performance.

DAT-based predictors demonstrated higher Cohen’s d values compared to those derived from traditional SST ( Figure 2 ). For DAT-based predictors, excluding NGR_d , Cohen’s d between the two classes increased progressively with greater mixing depths, peaking at depths of 100-110 meters. NGR_d , in contrast, displayed a steadily increasing Cohen’s d value with depth, reaching a peak at 60 meters and demonstrating a higher Cohen’s d value that overshadowed the other potential predictors. Figure 3 illustrates distinct patterns of correlation for each predictor as a function of mixing depth, indicating that the relationship between predictors and TC intensity change is sensitive to the mixing depth. Notably, NGR_d emerges as a superior predictor, with its maximum correlation coefficient occurring at a mixing depth of 60 meters ( Figure 3 , red line). This is not only higher than those of other predictors but also aligns with the depth where Cohen’s d—a statistical measure of effect size—reaches its peak ( Figure 2D , blue line). The consistency of the NGR_d peak with the maximum of Cohen’s d at 60 meters suggests a strong and possibly causal relationship between DATs of this depth and TC intensification rates, as well as RI events. This underscores the value of NGR_d , based on 60-meter DAT, as a potentially powerful single predictor for anticipating changes in TC intensity, which is crucial for early warning systems and preparedness measures in vulnerable coastal regions.

As outlined in Figure 4 and Table 5 , our study includes a comprehensive summary of the performance metrics — POD, PSS, FAR, Precision, ACC, and F-1 score — for individual ML models. These were evaluated during the training period running from 2014 to 2018. A modest change emerged when we incorporated NGRs into the predictor pools: the metrics of POD, Precision, PSS, F-1 score and ACC generally increased across the models, while the FAR metric correspondingly decreased. The only exception to this was observed in the DT model. This underscores the relevance and value of incorporating NGRs into the feature set, as models with NGRs consistently outperformed those without. Individually, the NGR-based FNN exhibited the highest predictive performance overall, closely followed by the SVM model.

During the test period of 2019 to 2021, summarized in Figure 5 and Table 6 , the favorable impact of NGRs was further corroborated. NGRs-based models once again outperformed their counterparts that lacked this feature. This improvement is related to the increase in the number of samples in TN. NGR-based models detected non-RI cases relatively better ( Figure 5 ). Among the individual models, the NGR-based SVM emerged as the best performer, particularly in terms of PSS. The consistency of this impact across both training and test periods reaffirms the generalizability and reliability of our methodological approach. An interesting point of divergence between the training and test periods was in the performance indicators. In addressing the class imbalance, oversampling was applied during the training phase. However, this method artificially inflates the TP count. The notable increase in POD and a corresponding decrease in FAR during the training period is attributable to the oversampling technique employed. Because oversampling is not performed in the testing phase, the ratio of RI cases decreases significantly compared to the training phase. This results in a relatively large decrease in TP, which in turn inflates FAR. This highlights the distortion in model performance metrics due to the uneven application of oversampling across training and testing datasets.

To generate an ensemble forecast, we employ a hard-voting method based on the collective performance metrics of five distinct classifiers. Our ensemble performance metrics for the training and test are shown in Table 5 and Table 6 , respectively. What becomes evident is that integrating NGRs into our ensemble model substantially augments its predictive capabilities. This improvement is noticeable during the test period. Notably, the PSS and F-1 score saw a 10% increase when NGRs were included in the ensemble model, demonstrating a more skillful forecast ( Table 6 ).

To contextualize our ensemble’s performance, it is useful to compare it with other contemporary ML-based RI forecasting models. Wei et al. (2023) presented a deep learning network model called TCNET, which they compared against two Statistical Hurricane Intensity Prediction Schemes (SHIPS)-based models (COR-SHIPS and LLE-SHIPS), along with other models from Yang (2016); henceforth referred to as Y16) and Kaplan et al. (2015); henceforth referred to as KRD15). Ko et al. (2023) explored the application of a consensus machine learning (CML) model in TC intensity change forecasting and indicated the CML exhibits better performance on RI predictions compared to the operational models such as SHIPS, GFS. Narayanan et al. (2023) proposed a simple deterministic binary classification model based on the co-occurrence of environmental parameters (MCE) to predict an RI event. Their results indicated that MCE shows improved skill over the decision tree and logistic regression models, with more accurate RI predictions in the overall testing dataset. The PSS values for these models, displayed in Figure 6 , show that our ensemble model (NGR-ENS), with a PSS of 0.56 and a POD of 0.79, surpasses all these competing models including TCNET (0.48), MCE (0.40) and CML (0.50). TCNET has the lowest FAR (0.43) followed by CML (0.50), LLE-SHIPS (0.56) and NGR-ENS (0.62). This holds even when considering different target periods or datasets. In essence, our ensemble approach fortified by the inclusion of NGRs offers superior predictive accuracy for RI events with an advantage of the noticeably high POD rate and the relatively low FAR rate.

Lee et al. (2016) suggested that the bimodal distribution of lifetime maximum intensity in TCs can be attributed to two distinct types of TC: those that experience RI (RI storms) and those that do not (non-RI storms). They showed that a significant majority—79%—of major TCs, those classified as category 3 or above, belong to the RI storm. Conversely, only a small fraction—6%—of non-RI storms ever escalate to become major TCs. Therefore, RI prediction performance in major TCs can represent the overall prediction performance. During the test period (2019 –2021), our ensemble model showed noticeable performance improvements when NGR was included as a variable ( Table 7 ). A recent C_d parameterization study showed that C_d decreases after saturating at 33 m s^-1, which leads to an increase in NGR, which can induce RI (Kim et al., 2022). These findings suggest that accurately simulating flux exchanges, especially in storms ranging from categories 1–3, can substantially enhance the model’s ability to predict RI accurately.