{'fraction_non_alphanumeric': 0.03566009104704097, 'fraction_numerical': 0.016795420057939027, 'mean_word_length': 4.795922446532081, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tropical Cyclone Intensity Evolution Modeled as a Dependent Hidden Markov Process

A hidden Markov model is developed to simulate tropical cyclone intensity evolution dependent on the surrounding large-scale environment. The model considers three unobserved (hidden) discrete states of intensification and associates each state with a probability distribution of intensity change. The storm's transit from one state to another is described as a Markov chain. Both the intensity change and state transit components of the model are dependent on environmental variables including potential intensity, vertical wind shear, relative humidity, and ocean feedback. This Markov environment-dependent hurricane intensity model (MeHiM) is used to simulate the evolution of storm intensity along the storm track over the ocean, and a simple decay model is added to estimate the intensity change when the storm moves over land. Data for the North Atlantic (NA) basin from 1979-2014 (555 storms) are used for model development and evaluation. Probability distributions of 6-h and 24-h intensity change, lifetime maximum intensity, and landfall intensity based on model simulations and observations compare well. Although the MeHiM is still limited in fully describing rapid intensification, it shows a significant improvement over previous statistical models (e.g., linear, nonlinear, and finite mixture models).


Introduction
Tropical cyclones (TCs) are one of the most destructive natural phenomena on Earth and are the leading cause of great social and economic losses. Accurately assessing TCrelated risk is of great importance. However, while TC intensity plays a large role in the damage potential of a storm, TC intensity modeling is still far from satisfactory (Rappaport et al. 2009). This is mainly because the intensity evolution of TCs involves complicated dynamic and thermodynamic processes that have not been fully resolved Wing 2016, Walsh et al. 2016). Apart from internal variability (e.g., eyewall replacement cycles; Willoughby et al. 1982), the storm environment plays a key role in controlling the storm's intensity evolution. The pre-storm sea surface temperature (SST) has been considered as one of the most important environmental variables (Whitney et al. 1997), since TCs draw their energy from the underlying seawater. Emanuel et al. (2004) used an ocean-coupled model to show that other environmental factors, including maximum potential intensity (MPI) and vertical wind shear, also have great influences on storm intensity. Furthermore, the interaction between TC and the ocean can significantly influence TCs' intensity. The entraining of cold, deep water induced by storm's surface winds may cause strong cooling of SST and thus negatively affect a storm's intensification (Schade and Emanuel 1999;Lloyd and Vecchi 2011).
Under the cooperative work of inner-core, oceanic, and large-scale environmental processes, many TCs undergo rapid intensification (Kaplan and DeMaria 2003) during their lifecycles. Rapid intensification (RI) is the significant strengthening of a TC in a short time (specifically, intensity increases by 30 kts or more over 24 hours), and almost probability distribution of intensity change with each state. Unlike the (statistical) mixture model, which assumes that the state identity for each time step is determined separately, the (stochastic) MeHiM considers the storm's transit from one state to another as a Markov chain.
Both the intensity change and state transition components of the MeHiM are dependent on identified environmental variables. To keep the model relatively simple, we apply only four environmental variables (MPI, vertical wind shear, relative humidity, and an ocean feedback parameter) and two storm variables (current intensity and previous intensity change). Few of the developed ocean feedback parameters (e.g., Schade and Emanuel 1999, Vincent et al. 2012, Balaguru et al. 2015, which indicate ocean-storm interaction, have been evaluated in statistical modeling of TC intensity. Here we adopt the ocean feedback parameter of Schade and Emanuel (1999), which has also been used in a simplified dynamical model (Emanuel 2017), and test its impact on the overall performance of the MeHiM.
The MeHiM is developed to predict the intensity change of storms moving over the ocean.
To simulate the intensity evolution over storm's entire lifecycle, we add a simple decay model to predict a storm's intensity change over land. Here we develop the model for the North Atlantic basin, though the same methodology can be applied to other basins or on the global scale. To evaluate the developed model, we apply it to simulate historical storms using historical tracks and environment taken from reanalysis datasets. The performance of the model is measured by its capability to reproduce the observed TC intensity metrics, including spatial distribution of intensity over the basin, probabilistic distributions of intensification rate, LMI, and landfall intensity, as well as temporal evolution of intensity for sample cases. We focus particularly on the model's capability to capture the extremes related to RI.
The MeHiM improves on the mixture model of Lin et al. (2017), because 1) it accounts for the (possibly large) temporal correlation of the states (e.g., slow, moderate, rapid) of intensity change over storm's lifecycle and 2) it integrates the mixture model of intensity change and the identification model of intensification state within a Markov framework, avoiding statistical inconsistency between the two components. It is interesting to ask whether these improvements in methodology will lead to significant improvements in the model performance, especially in capturing RI related extremes. To answer this question, we build the mixture model (with an added identification model) using the same probability distributions (i.e., Gaussian for intensity change and Logistic for state identification), environmental and storm variables, and observational datasets as in the MeHiM. In addition, we build the (Gaussian) linear regression model using the same variables and datasets to provide a baseline. The linear and mixture models are run in parallel with the MeHiM, and the simulation results are compared among the models and with the observations. This study is organized as follows. Following this introduction, the development of the models (i.e., linear regression model, mixture regression model, MeHiM, and land decay model) is described in Section 2. Model simulation results are shown and compared in Section 3. Section 4 discusses the effects of adding the ocean feedback parameter and correctly identifying the extreme state of intensity change in the MeHiM. Section 5 summarizes the main findings of the study.

Model Development
Three statistical TC intensity models are developed to predict 6-hourly storm intensity change DV: the multiple linear regression model, mixture model (with an added state identification model), and MeHiM. We build these models on the same six variables: last step intensity change (DV p ), the current intensity (V), maximum potential intensity (MPI), vertical wind shear (SHR), relative humidity (RH), and an ocean predictor (OCN). Apart from these intensity models, which are used to simulate a storm's intensification over the ocean, we also add a simple land model to simulate the decay of the storm intensity over time when the storm moves on land.
All models are developed based on historical records from 1979 to 2014 in the North Atlantic basin. The TC dataset is taken from the IBTrACS WMO archive (Knapp et al. 2010). It includes 6-hourly latitude and longitude positions and 10-minute maximum sustained wind speeds at 10 m above the sea surface. The location data are used to calculate the translational velocity of the storm. The storm maximum wind intensity (V) is estimated by subtracting the surface background wind, estimated as a fraction of the storm translation velocity (Lin and Chavas 2012), from the observed maximum wind.
The 6-hour intensity change DV (and DV P ) is then obtained from the time series of the estimated storm maximum wind intensity.  (Dee et al. 2011)]. A storm's maximum potential intensity, MPI, is theoretically derived following Emanuel (1995) and Bister and Emanuel (1998;2002). Similar to that in previous studies (e.g., Kaplan 1994, Lin et al. 2017), the deep layer vertical wind shear, SHR, is defined as the difference between the 850-and 200-hPa level winds, averaged over a 200-800 km annulus around the storm center.
High-level relative humidity, RH, is computed as the averaged relative humidity in the layer between 300-and 500-hPa within a 500-800 km annulus around the storm center.
In addition to the atmospheric predictors, we also incorporate an oceanic variable (OCN), developed by Schade and Emanuel (1999) based on numerical modeling, to represent the ocean's negative impact on storm intensification. OCN is defined as where ! is the storm translation speed, ℎ ! is the ocean mixed layer depth, and Γ is the thermal stratification below the ocean mixed layer. OCN is a number between 0 (strong weakening influence on the storm) and 1 (weak influence); it is closer to 0 when the storm is stronger, moving more slowly, with a shallower mixed layer depth, and/or with stronger thermal stratification. OCN is more sophisticated than other ocean variables such as the ocean temperature averaged over the mixing layer and the ocean heat content, which indicate only the amount of heat stored in the ocean and are incapable of representing the interaction between the storm and ocean. The ocean heat content is used in SHIPS, but the data are available only for a limited time period and over a limited part of the Atlantic basin (Mainelli et al. 2005); the ocean temperature averaged over the mixing layer was found to be an insignificant predictor of TC intensification (Lee et al. 2015 The model coefficients and are estimated by minimizing the squared residual error (i.e., OLS). The estimated model coefficients as well as the standard deviation ( ) are shown in Table 1. The signs of the coefficients in front of all variables are as expected (based on previous studies; e.g., DeMaria and Kaplan 1999, Schade andEmanuel 1999, Lin et al. 2017), i.e., DV is positively (negatively) correlated to DV p (V); it is expected to increase when MPI and RH are larger, SHR is smaller, and/or ocean cooling is weaker (OCN larger). The storm's last-step intensification DV p is shown to be the strongest predictor, followed by the current intensity (V). Among the environmental variables, MPI and SHR are more important than RH and OCN. The intercept ( ) is nearly zero, which means that the expected value of 6-h intensity change DV is nearly 0 when all five predictors are at their means (i.e., standardized values being zeros).

Finite Mixture Regression (FMR)
The Finite Mixture Regression (FMR) extends the OLS by incorporating a heterogeneity structure to the model. In this model, each observation is considered to be a random sample generated by one of its k components (groups/clusters), while the identity of the generating component is not observed. In a Gaussian setting, each of the k components is a linear regression, and the probability distribution of the predictand Y is mixed Gaussian, where ! is the component mixing proportion/weight, ! and ! are regression coefficients, and ! is the standard deviation, for each component . These model parameters can be estimated from data using the maximum likelihood estimation method with the generalized expectation-maximization (EM) algorithm (McLachlan and Krishnan 2007). In addition to model parameters, the EM algorithm also estimates the probability that each observation was generated from each of the k groups, providing an optimal group assignment for each observation.
Based on the same dataset as for the OLS, the estimated model parameters for FMR with three groups (k=3) are shown in Table 2(a). The obtained group weights and features are similar to those from the FMR analysis of Lin et al. (2017). Specifically, over 60% of the observations are optimized by the EM algorithm into Group 2, nearly 30% into Group 1, and only less than 10% into Group 3. Groups 1, 2, and 3 again feature "static," "normal," and "extreme" intensity change (DV), respectively, in the sense that, first, the mean of the intensity change as a response to the environment is modeled to be weak for Group 1, relatively strong for Group 2, and the strongest for Group 3 (as indicated by the magnitude of the model coefficients in front of the environmental variables); and, second, the model standard deviation (representing the width of the subgroup Gaussian distribution) increases significantly from Group 1, to Group 2, and to Group 3. As a result, the sample mean of DV increases moderately in absolute value from Group 1 (-0.07), to Group 2 (0.82), and to Group 3 (-3.86) and the sample standard deviation of DV increases substantially from Group 1 (0.47), to Group 2 (5.89), and to Group 3 (16.73), indicating that (regardless of the environmental conditions) small intensity change observations are assigned into Group 1, extreme intensification (and deintensification) into Group 3, and the rest into Group 2 (according to the optimal assignment based on probability). Such grouping features will be further discussed with the Markov model analysis. It is also noted that some unexpected signs of the model coefficients are obtained in the FMR modeling, e.g., negative response to MPI as a statistical compensation for the large negative response to SHR in Group 3, which will be shown to be improved in the Markov model.
Group 1 vs. belonging to Group 3, and, similarly, the linear combination of covariates for Group 2 (second row) describes the log odds of belonging to Group 2 vs. belonging to Group 3. For example, a one-unit, or one standard deviation, increase in the variable DV p is associated with a decrease in the log odds in the amount of 0.359, which is equivalent to a decrease in the odds by approximately 30%. At the same time, an increase in the variable DV p is also associated with a decrease in the odds of belonging to Group 2 vs.  Table 1) for larger DV, the probabilities that it belongs to Groups 1, 2, and 3 are 0.39, 0.60, and 0.01, respectively. When all variables of the new observation are at the 95 th percentile unfavorable for larger DV, the probabilities that it belongs to Groups 1, 2, and 3 are 0.23, 0.74, and 0.03, respectively. The probabilities obtained in the three scenarios are relatively close, indicating that the model may be limited in capturing the change in the membership probability as a response to the change of the predictor variables. Moreover, the probability of belonging to the extreme group (Group 3) is always significantly smaller than the group weight from the optimization (~ 0.07; Table 2a). This limitation in capturing the extreme state by the FMR classification model was also noted by Lin et al. (2017). This limitation is reduced in the Markov model.

The Markov environment-dependent Hurricane Intensity Model (MeHiM)
The MeHiM is based on the idea of the hidden Markov model (HMM; Rabiner 1989, Cappé et al. 2006, an effective tool for modeling stochastic processes. In an HMM, the visible observations are distributed according to the underlying hidden states, where the states are discrete and are assumed to form a first-order Markov process. With a selected number (k) of states, an HMM can be fully determined by the following three sets of parameters: 1. Initial state probability ! : the probability of starting from state .
2. Transition probability !" : the probability of moving from state to .
3. Response distribution ! ! : the distribution of the response variable in state .
The HMM can be extended to the dependent HMM with the above model parameters dependent on covariates (Zucchini et al. 2016). Then, the likelihood of the response observation sequence !:! = ! ! … ! , given the covariates sequence !: where ! denotes the state at time , and all three model components (i.e., initial state probability, transition probability, and response probability) depend on the covariates.
For the entire dataset, the total likelihood is the product of the individual likelihood functions for all time series.
The basic idea of a dependent HMM is applied in the MeHiM to model the sequence of TC intensity change (DV) over the storm's lifetime dependent on the selected storm and environmental covariates (i.e., DV p , V, MPI, SHR, RH, and OCN). Motivated by the FMR analysis, we assume the storm in one of three hidden states ( = 3) at each time step. Given the state, the response distribution is assumed to be Gaussian, i.e., dependent on the 6 predictor variables as in the OLS and FMR models.
Similar to the FMR analysis, we apply the multinomial logistic regression to model the transition probability that the storm moves from the present state to the next state, i.e., The transition model (3×3 probability matrix) contains three multinomial logistic regressions for each of the three states, and the probabilities are dependent on all six variables as in the FMR classification analysis. The initial state probabilities are also similarly determined by a multinomial logistic regression.
The parameter fitting process involves the backward-forward algorithm for likelihood calculation and the Baum-Welch method, an iterative procedure for optimization similar to the EM method. In addition to model coefficients, the optimal state sequence can also be estimated through the Viterbi algorithm (Zucchini et al. 2016). A more detailed description of the algorithms can be found in the statistical R-package depmixS4 (2010: depmixS4: an R-package for hidden Markov models. 36.). To prepare the temporal sequences for model fitting, it is intuitive to assume that each storm contributes to a single time sequence. However, since a TC behaves differently when over land, we filtered out these observations by dividing the storm's intensity evolution into multiple temporal sequences. That is to say, the time series may begin with either the genesis of the storm or the first observation after the storm makes landfall and moves back to the ocean again. The time series ends when the storm dies or hits land. However, the temporal sequences in between can be very short (sometimes less than 3 observations) when a storm hits land twice within a short period of time. These short sequences are considered as the "noise" of the training data and thus are not used. After removing sequences that are shorter than 12 observations (i.e., 3 days), 10072 observations, which consist of 357 time series, remain for model training.  Table 3a, the model standard deviations of the three states show similar patterns with the three subgroups in FMR (Table 2a). Similar explanations can be applied here so that the three underlying states represent storm's "static" (State 1), "moderate" (State 2), and "extreme" Other estimated parameters of the response component of the MeHiM are also shown in Table 3(a). Improved over the FMR results (Table 2a)

Land Decay Model
TCs weaken quickly over land mainly due to the lack of fuel sources and the entrainment of dry air. Estimating the decay of storms after making landfall or when moving across small islands is also a key component of TC intensity models. Georgiou (1985)  shown that the inclusion of such decay models can improve intensity forecasts in multiple basins (DeMaria et al. 2006).
Here, we develop a land model based on the same dataset applied for the over-ocean intensity models. The model structure is similar to that in Demaria and Kaplan (1995), except that we set the reduction factor to be 1 for simplicity, ignoring the assumed sudden intensity decrease when storms approach land. In this model, the storm's intensity exponentially decays as a function of time , as shown in Eq. (6).
where ! is the storm's original intensity prior to landfall, which decays with a constant temporal rate , to a constant background wind, ! , that a storm may maintain over land.
To estimate and ! in this model, we collect all TC intensity sequences when the storm's center is over land. While most of these segments are quite short, with an average length of about 5 observations, we first exclude segments that are shorter than 2 observations (12 hours). Also, as the model does not apply for storms weaker than ! , we assume the storms will not decay if the landfall intensity is smaller than ! . Accordingly, we exclude segments with ! smaller than 20 kt with try and error. The model parameters and ! are then fitted based on 229 segments, and the optimal decay rate and background wind are estimated to be 0.049 kt/6h and 18.82 kt, respectively, by nonlinear least squares regression. The fitted model predicts the intensity decay over land that compares well with the observation for storms with various initial intensities, as shown in Fig. 3.

Model Evaluation
To evaluate the developed models, we perform MC simulations of the intensity evolution over storms' lifetime for historical tracks and compare model simulations with observations. We examine the models' capability to reproduce the observed TC intensity metrics, including spatial distribution of intensity over the basin, probabilistic distributions of intensification rate, LMI, and landfall intensity, as well as temporal evolution of intensity for sample cases.

Monte Carlo (MC) simulation
For all simulations, the storm is initialized by the first two steps of observed intensity to obtain initial and DV p . Then, for the OLS, the intensity change for each step is obtained by sampling from a Gaussian distribution (Eq. 2 and Table 1). For the FMR, the classification model is first applied to compute the mixture weights (Eq. 4 and Table 2b), based on which group membership is determined with the discrete acceptance-rejection method. The intensity change is then obtained by sampling the Gaussian distribution from the selected group (Eq. 3 and Table 2a).
Simulation with the MeHiM has one more step in the initialization. Given the environmental parameters (i.e., MPI, SHR, and RH) at the storm's occurrence, we use the initial state model with the acceptance-rejection method to determine the state for the first two steps. Then, we continue the simulation by repeatedly choosing the next state to visit according to the transition probabilities (Eq. 7 and Table 3b) with the discrete acceptance-rejection selection and updating storm intensity with the intensity change sampled from the state's corresponding response distribution (Eq. 6 and Table 2a).
During all simulations, the land model is applied when the storm's center hits land, and the over-the-ocean intensity models are resumed if the storm moves back to the ocean.
Storms with ! smaller than ! are simulated to keep their original intensity over land, while storms with ! greater than ! are simulated to decay according to the filling model (Eq. 8). For the MeHiM, the storm's states over land are assumed to be the same as the last state when it is over ocean, before hitting the land. Simulations stop at the end of the historical record of the track or when the storm' intensity becomes lower than 10 kt.
Each historical storm is simulated 100 times with each of the three intensity models, and the comparisons with observations are based on all realizations. We also compare the histograms of 24-h intensity change (24-h DV) between observations and simulations (Fig. 6) Particularly, although the MeHiM does not capture the outlier extreme 6-h intensification (Fig. 5b), it produces larger, more extreme 24-h intensifications and higher probability densities for the bins of 24h DV greater than 30 kts (i.e., RI) than IBTrACS (Fig. 6b). On the other hand, the MeHiM underestimates both the rates of most extreme 6-h (Fig. 5b) and 24-h (Fig. 6b) deintensifications. The large 24-h intensifications are likely generated by the "lock-in" mechanism in the MeHiM, which supports continuous large intensification once the storm enters the extreme state and the environment is favorable.

Spatial distributions
The MeHiM lacks such a mechanism to support continuous large deintensification, as discussed earlier.
Next, we compare the probability distribution of lifetime maximum intensity (LMI) in the simulations and observations (Fig. 7). Among the three models, the MeHiM compares with IBTrACS the best. It captures the shape and peak of the LMI distribution quite well, although it is smoother and misses the "fluctuations and bumps" in the tail, and as a result it underestimates the extreme LMIs to some extent.  (Fig. 6b), but it still underestimates the tail of the LMI distribution ( Figs. 7 and 8). This underestimation of extreme LMIs also explains the underestimation of spatial extremes in Fig. 4.
Then, we compare the simulated and observed probability distribution of the landfall intensity for the North Atlantic coastline (Fig. 9). Both the OLS and MeHiM capture the landfall probability distribution well, with the MeHiM slightly better in capturing the tail.
FMR underestimates the tail more than the MeHiM and OLS and significantly overestimates the peak. Also, while the 15 th -85 th percentile bond of the OLS and MeHiM estimates can generally cover the distribution from IBTrACS, the percentile bond of FMR estimates is quite narrow, and it often misses IBTrACS. It is noted that for the Northeastern US coastline, the distribution from IBTrACS is quite rough and the percentile bond of the model estimates is relatively large, due to data limitations for both the observations and simulations for this region.

Sample cases
Apart from capturing various statistical distributions, we compare the performances of the OLS, FMR, and MeHiM in simulating the intensity evolution of individual storms. However, even for such an extreme case, we find that the ensemble of the MeHiM is able to cover the observed LMI -the maximum (over the 100 realizations) of the simulated LMI is 175 kt from the MeHiM, while it is 135 kt from OLS and 147 kt from FMR. In addition, about 40% of the realizations from the MeHiM can reach a LMI over 100 kt, which is much more than that from OLS and FMR (~15%). Indeed, an upward shift over the early development stage in the MeHiM realizations (indicated by the shadings) indicates a larger simulated intensification rate in the MeHiM compared to the OLS and FMR. Hurricane Wilma weakened greatly in the later stage of its lifecycle, when its intensity decreased by 35 kts within 24 hours. The extreme cases of the MeHiM realizations (shading) are able to represent this rapid deintensification and large decay.
In Fig. 9(g-i), an unnamed Tropical Depression in 2003 is shown as an example of weak storms with LMI less than 45 kt. In this case, both FMR and the MeHiM are capable of simulating the intensity evolution, with the ensemble mean matching the observation.
The OLS, on the other hand, significantly overestimates storm's intensity, with mean as well as the 20 th -80 th percentile envelope shift upward away from the observation. This overestimation by the OLS is typical for weak storms, as discussed before (see Figs. 7 and 8).
From the above results, we find that both FMR and the MeHiM are able to simulate relatively weak storms, while the OLS tends to overestimate weak storms. The MeHiM has significant advantage over both FMR and the OLS when simulating RI storms or storms with large intensification or deintensification. It is also noted that the models tend to overestimate the storm's intensity after making landfall (e.g., Figs. 10a-c), which indicates that the simplified land model may not be sufficient to simulate a storm's rapid decay over land.

Effect of the ocean feedback variable
The three (over-ocean) intensity models are developed with four environmental covariates including three atmospheric and one ocean feedback variable. As the ocean feedback variable (OCN) is not used in previous statistical models, here we test its significance and examine how OCN represents TC-ocean interaction. We use two storms as examples: Tropical Storm Larry (2003) and Hurricane Felix (1995). The upper two panels in Fig. 11 show simulated intensity evolution for these two storms from the original MeHiM (as discussed above), and the lower two panels show simulations from  Fig. 13 (a-c). Fig. 13 ( Future studies may also include a comparison between the MeHiM and dynamic models, such as the Coupled Hurricane Intensity Prediction System (CHIPS; Emanuel et al. 2004) and its simplified algorithm (Emanuel 2017). The ultimate goal of our study is to integrate the MeHiM with a genesis model and a track model into a complete hurricane climatology modeling system, similar to Emanuel et al. (2008) and Lee et al. (2018).   To