Combining Monte Carlo and Ensemble Probabilities in Tropical Cyclone Forecasts near Landfall

Accurate forecasts of tropical cyclones (TCs), especially landfalling hurricanes, are crucial for reducing the damage they cause (Smith, 2017). Despite progress in numerical weather prediction (Bauer et al., 2015) and statistical forecast (DeMaria et al., 2009), forecast error mitigation is difficult due to the complexity of TC evolution and the lack of observations (Rappaport et al., 2009; Emanuel and Zhang, 2016). To overcome forecast uncertainty and provide early warning for TCs, a TC probability program was implemented in 1983 (Sheets, 1985) with support from the National Hurricane Center (NHC). With this program, 10-yr samples of official forecast errors are composited relative to the forecast position at the discrete forecast verification times of 12, 24, 48, and 72 h (up to 120 h after 2003). These provide information on the likelihood of a storm of interest moving over a given area. Following this effort, DeMaria et al. (2009) developed a Monte Carlo probability (MCP) model to generate a TC wind probability map at individual locations during a 0–120-h forecast period. In this wind probability model, random samples from the past five years of forecast errors [usually with a sample size of 1000, as suggested by DeMaria et al. (2009)] are required to generate cyclone track, intensity, and corresponding radii realizations [with the simple climatology and persistence model (radii-CLIPER); Knaff et al., 2007] based on an official forecast track from NHC. These generated realizations can then be used to produce a probability map in a single forecast result by counting the number of realizations for which the point is inside the radius of the wind speed threshold of interest (e.g., 34, 50, and 64 kt, where 1 kt = 0.5144 m s⁻¹).

To narrow the forecast track error distribution, after 2010, the NHC introduced the Goerss predicted consensus error (GPCE; Goerss, 2007) into the MCP model to generate a more relevant wind probability map. The different GPCE selections at different forecast times lead to a slight improvement in MCP model forecasts in most cases (DeMaria et al., 2013). Meanwhile, because of well-distributed historical errors with a mean value of about zero (DeMaria et al., 2009), the MCP model generates track realizations that are always distributed along the given forecast track. When there are strong spreads between the forecast and best track [i.e., the poststorm best estimates of the cyclone track and intensity, introduced by Jarvinen et al. (1984) and Knapp et al. (2010)], the 1000 generated track realizations can be identified as inaccurate realizations that lead to poor TC realizations and result in poor wind probability forecasts. In order to reduce the forecast errors, more forecast tracks and errors should be considered and added to the MCP model.

Besides the MCP model, major operational numerical weather prediction centers, such as the United States’ NCEP, the ECMWF, and the United Kingdom Met Office (UKMO), all provide reliable TC ensemble forecasts. While these ensembles can improve probabilistic forecasts for severe weather events (Froude, 2010, 2011; Frame et al., 2011; Swinbank et al., 2016), wind probability maps (Matsueda and Nakazawa, 2015) for hurricane forecasts can also be generated by combining ensemble forecasts from multiple sources (e.g., from ECMWF, UKMO, and NCEP).

In this study, we develop a modified MCP model by combining probability maps from both MCP and multi-model ensemble forecasts to improve TC probability forecasts near landfall. Based on the MCP model using the NCEP Global Forecast System (GFS) forecast as the official forecast track, we first evaluate and adjust the MCP method by limiting the unrealistic TC moving speeds and directions to be more suitable for TCs near landfall. Then, we combine it with ensemble forecast maps from ECMWF, UKMO, and NCEP to further mitigate the unrealistic TC moving directions to improve hurricane 0–120-h cumulative wind speed probability forecasts. It is our purpose to explore a feasible method of combining MCP and ensemble forecasts for accurate TC wind forecasts near landfall. Hurricanes Sandy (2012), Irma (2017), Harvey (2017), and major hurricanes during the 2018 season over the Atlantic Ocean are used as cases for this study. The MCP model, ensemble forecasts, and modifications of the MCP model are described in Section 2. Probability forecast results and statistical evaluation of cases are provided in Section 3. Section 4 includes an evaluation based on major hurricanes during the 2018 season. Concluding remarks are presented in Section 5.

The MCP model estimates the probability of wind speed at individual locations during a TC forecast period (e.g., 0–120 h) based on randomly sampled errors from previous years (historical) of official track errors by adding these errors to the official deterministic forecast positions. The sample error distributions are calculated by comparing the official forecast positions to the NHC “best track” positions, which are the poststorm best estimates of the cyclone track and intensity from observations. Therefore, with the MCP model, random samples of track and intensity forecast errors from operational forecasts, such as those from the NCEP GFS, in the previous five years are added to a single, deterministic forecast (e.g., official GFS forecast) to generate 1000 realizations of cyclone track, intensity, and corresponding radii. The model produces a probability map by counting the number of realizations for which the points are inside the radius of the wind speed threshold of interest (e.g., 34, 50, and 64 kt). Commonly, the NHC generates 0–120-h official wind probability maps based on its official cyclone forecasts (DeMaria et al., 2009, 2013).

TC track errors are defined as the great-circle distance from the forecast track to the best-track positions. According to DeMaria et al. (2009), track errors can be decomposed into along-track (AT; positive when the forecast position is ahead of the best-track position) and cross-track (CT; positive when the forecast position is to the right of the best-track position) errors. Since these AT and CT errors (i.e., ATE and CTE) at a certain time can be affected by forecast errors at earlier times (DeMaria et al., 2009), directly adding these randomly sampled past ATE and CTE to a single forecast can lead to unrealistic tracks. Therefore, a serial correlation is used to eliminate this problem and generate realistic track forecasts. DeMaria et al. (2009) used a linear fit method to determine the influence of an error at a previous time on ATE and CTE:

where a, b, c, and d are fitting coefficients determined by the input historical ATE and CTE, and the subscript t is the forecast time. Then, based on Eq. (1), the historical ATE (CTE) can be divided into two parts: the influence of previous errors (i.e., calculated errors), which is calculated from Eq. (1) with coefficients a, b, c, and d; and the residuals, which are the historical ATE (CTE) minus the calculated values from Eq. (1). For t = 12, the coefficients a and c are equal to zero, and the ATE (CTE) calculated from Eq. (1) is equal to coefficient b (d), which is the mean of the historical ATE (CTE) at this time. The residuals are generated from the differences between the historical ATE (CTE) and the coefficient b (d). From t = 24 to t = 120, the residuals are generated from the differences between the historical ATE (CTE) and calculated errors from Eq. (1). These residuals at each time are the samples of the following MCP realizations.

After the historical ATE (CTE) is transformed into residuals as discussed above, the residuals are randomly sampled at each forecast time to generate the final ATE (CTE) of MCP and make the MCP forecast. At t = 12, The residuals are randomly sampled and added to Eq. (1) with coefficients b₁₂ and d₁₂ as the sample mean AT and CT track error biases at 12 h (as the error is equal to zero at t = 0; a₁₂ = c₁₂ = 0). Then these values of ATE (CTE) are input into Eq. (1) to obtain the calculated errors and added to the randomly sampled residuals at 24 h to generate the final ATE (CTE) at 24 h. Following a similar step, the final ATE (CTE) at 24 h is again input into Eq. (1) and combined with the residuals to generate the final ATE (CTE) at the subsequent forecast time. Finally, with the forecast errors from the past five years, a total of 1000 × 10 final ATE (CTE) are generated. These final errors are added to a single forecast track, thus producing 1000 hurricane tracks for the 120-h forecast.

Similar to track realizations, TC intensity errors in the previous five years are required as random samples for intensity realizations. In the MCP model, intensity is represented by the maximum wind speed at 10-m height. Therefore, the intensity error is the maximum wind speed error (SE). Meanwhile, the SE at any given time can be influenced not only by earlier SE but also by forecast cyclone maximum wind speed (SM) and the distance from the cyclone center to land (D, which is positive when the cyclone is over the ocean) at a given time. Considering the small influence of land when a cyclone is sufficiently far away from land, DeMaria et al. (2009) set D = 500 km when D > 500 km. Therefore, a linear relationship is assumed between the SE 12 h previous and the forecasts of SM, D, and SE at time t (DeMaria et al., 2009), thus:

where e, f, g, and h are least squares fitting coefficients, determined by the historical SE, SM, and D, and the subscript t is the forecast time. The historical SE at 12–120 h is transformed into a residual SE by subtracting the calculated SE (which represents the influence from the error in the previous 12 h) using Eq. (2). Similar to the track realizations, at t = 12, the final SE is equal to the sum of the residual and the coefficient h. The final SE at t = 12 is then input into Eq. (2) to calculate the SE at t = 24. Then, this calculated SE is combined with randomly sampled residual SE to generate the final SE at t = 24. Similar steps continue and generate the final SE at subsequent forecast times, and these final SEs are then added to the official (deterministic) forecast intensity to produce 1000 hurricane intensity realizations for the 120-h forecast at 12-h intervals.

Furthermore, when the underlying surface, either ocean or land, is different in the realization than in the forecast, the difference can lead to an incorrect SE, and consequently to poor intensity realizations. Therefore, two methods are used to modify these poor intensity realizations produced by the mismatched underlying surfaces between realizations and forecasts. When the deterministic forecast position is inland and the MCP realizations are over the ocean, SM_t in Eq. (2) will be replaced by the value over the ocean at the nearest forecast time; when the forecast position is over the ocean and the MCP realizations are inland, SM_t will be replaced by the value from an empirical inland wind decay model, as suggested by Kaplan and DeMaria (1995):

where SM₀ is the MCP forecast intensity at the first time when the TC moves inland.

The radii-CLIPER model suggested by Knaff et al. (2007) is used to generate radius realizations in the MCP model:

where V is the wind speed, r is the radial distance from the cyclone center, θ is the azimuth measured counterclockwise relative to 90° to the right of the cyclone’s motion direction, V_m is the maximum wind speed, r_m is the radius of the maximum wind, A is an asymmetry factor, and θ₀ is a constant that allows the maximum wind speed to be located at an azimuth other than 90° to the right of the motion direction. The values of A (= 1.06) and θ₀ (= 17) are both suggested by Knaff et al. (2007).

The cyclone size parameter x and r_m can be estimated by a simple climatic model (Knaff et al., 2007):

where φ is the latitude and the subscript c denotes a climatic estimate. With historical V, r, V_m, r_m, and θ, the best fit x for each historical hurricane forecast record can be produced with Eq. (4). After that, the best fit x is transformed to residual x by subtracting the climatic-estimated x_c from Eq. (5). Then, similar to the track and intensity realizations, the final x for each forecast time can be produced by adding the randomly sampled residual x to the climatic-estimated x_c in Eq. (5). Meanwhile, since the climatic-estimated x_c is often unequal to the real x at t = 0, when there is no residual x for realizations, this mismatch between climatic-estimated x_c and real x may transmit to subsequent realizations. Therefore, as suggested by DeMaria et al. (2009), exponential decay with an e-folding time of 32 h is used in radius realizations for each forecast time to cover the mismatch of climatic-estimated x_c and simulated x at the forecast start time. Finally, after these transformations and modifications, the final x for MCP is produced with climatic-estimated x, randomly sampled residual x, and the start-time x error. The azimuth-dependent wind profile for each realization is produced through Eq. (4) with this final x, hurricane maximum wind speed V_m, maximum wind radius r_m, and wind speed thresholds of interest (e.g., 34, 50, and 64 kt).

In this study, we use the NCEP GFS deterministic forecasts as official forecasts in the MCP model to generate the 0–120-h probability forecasts near and during landfalls of three major hurricanes: Hurricane Sandy, as an extratropical cyclone during its landfall from 0000 UTC 27 to 0000 UTC 29 October 2012; Hurricane Harvey from 0000 UTC 25 to 0000 UTC 29 August 2017; and Hurricane Irma from 0000 UTC 8 to 1200 UTC 10 September 2017. The GFS forecasts from the previous five years for each case, provided by the Automated Tropical Cyclone Forecast (ATCF) database, are used to provide the historical forecast errors for the MCP model.

For track realizations, Table 1 shows the coefficients a, b, c, and d in Eq. (1) obtained from the previous five years, namely, 2012–2016 TC data for Hurricanes Harvey and Irma (2017). A determination coefficient (R²) is also defined to measure the correlations between ATE_t (CTE_t) and ATE_t+12 (CTE_t+12). We found that there are considerable small determination coefficients (R²) at some lead times, for example, R² < 0.5 at 24, 84, and 108 h in our case. This indicates that there might be some track data that fail to satisfy Eq. (1). Namely, in some cases, early ATE (CTE) may not affect the following ATE (CTE). For instance, when a TC undergoes a large change in moving direction during 24 h, the rapid change in its track means that there is no significant relationship between ATE (CTE) in the first 12 h and those in the subsequent 12 h. Considering this problem, we divide the TC tracks into two groups based on direction changes (dγ) between a certain time t and t + 12, namely, the cases of dγ ≤ 45° or dγ > 45°.

Table 2 shows the parameters obtained from the TC data during 2012–2016 for cases of dγ ≤ 45° and dγ > 45°. The determination coefficient for dγ ≤ 45° is larger than that in Table 1, while it is close to zero for dγ > 45°. The results here confirm that Eq. (1) cannot be well applied in the track realizations when the subsequent forecast track has large directional changes. For instance, there are 21 times with dγ > 45° from 5 × 120 h GFS track forecasts, and the track realizations using Eq. (1) based on these GFS tracks will lead to poor wind speed probability forecasts. Therefore, we made an adjustment to the MCP model to ensure that the final ATE and CTE are randomly sampled directly from the past five years of GFS ATE and CTE without applying Eq. (1) when dγ > 45°.

Additionally, the moving speed of a TC at a certain time is usually related to its speed at an earlier time, and the forward speed will not increase quickly during a short period. Thus, some randomly sampled ATE and CTE could lead to track realizations with an unrealistic forward speed. To mitigate these unrealistic track realizations, we introduce a coefficient r_dis, defined as the motion distance ratio, namely:

where dis_t and dis_t+12 are the distance traveled between forecast intervals t − 12 and t, and between t and t + 12, respectively; r_dis denotes the ratio of TC motion distance dis_t+12 and dis_t. Based on data from the previous five years, we found that r_dis has a maximum and changes with dis_t. The original MCP model generates some unrealistic track realizations with r_dis over the maximum in historical data. Therefore, in the adjusted MCP model, a maximum r_dis limit is introduced into track realizations to eliminate these unrealistic tracks. For each forecast time, a series of maximum r_dis in the past five years is determined at different dis_t. The archived best track can be accurate only at 0.1° for latitude and longitude, namely, 0.1° of the great circle for the cyclone around the equator. Therefore, the dis_t interval for which the maximum r_dis is determined is set as 0.1° of the great circle on the earth (11.13 km). Then for each MCP forecast, we calculate the forecast r_dis and compare it with the corresponding threshold of the maximum r_dis to eliminate excessive points. For example, for dis₁₂ of 1.5°, the realizations with dis₂₄ of 6.0° are eliminated by the maximum r_dis of 2.4.

With the mentioned adjustment, the final ATE (CTE) is determined through the MCP model and added to the GFS deterministic track forecasts to generate the 0–120-h track realizations for each hurricane with time intervals of 12 h. For the intensity realizations, there is no adjustment to the MCP model. The final SE is determined by the MCP model with GFS forecast data from the previous five years. It is then combined with the GFS deterministic track, intensity forecasts, and track realizations to produce the intensity realizations.

For radius realizations, there are negative x residuals in the GFS data from the past five years that can lead to a small final x and large unrealistic radii. To eliminate these excessive radii, in the modified model, radius realizations are limited by maximum radii of 34, 50, and 64 kt at different forecast times based on cyclone samples from the past five years. Furthermore, due to the underlying surface friction of land, the wind speeds of major hurricanes such as Sandy, Irma, and Harvey decrease very quickly during landfall. Data for standard radii of 34, 50, and 64 kt are not available, resulting in a lack of data for evaluation. To better assess the MCP forecast, we also forecast the impact radii for lower wind speeds, such as 20 and 27 kt. Similar to 34–64 kt, the radii for 20 and 27 kt are limited by the maximum radii that are determined by the regression with radii data for 34, 50, and 64 kt. With this correction, the final x is applied to the GFS track and intensity realizations to generate radius realizations for each hurricane. For convenience, the adjusted MCP model is referred to hereafter as the AMCP model.

Considering that hurricanes can move several hundred kilometers within 12 h, directly assessing wind speed probability with 10 AMCP forecasts (12–120 h with 12-h intervals) could miss the influence of cyclone motion during the 12-h time period. For instance, the radii for 64-kt wind usually measure less than 40 km from the cyclone center. If we directly use the 12-h interval realizations, such as 12–24 h, to calculate the probabilities of 64-kt wind for any point outside the radius of 40 km, the probability will be zero and these points will be omitted in the calculation. Therefore, to reduce this omission, DeMaria et al. (2009) linearly interpolated the realizations into 2-h intervals, during which the cyclones, especially after landfall, usually do not move far. Following DeMaria et al. (2009), the realizations here including track, intensity, and radii are linearly interpolated at 2-h intervals within the realizations at 12-h intervals. For instance, the track, intensity, and radii at 12 and 24 h are used to linearly interpolate the track, intensity, and radii at 14, 16, 18, 20, and 22 h. Then, AMCP at any given location and for any forecast interval to 120 h can be estimated by counting the number of realizations for which the point is inside the radius of the wind speed threshold of interest, relative to the total number of realizations. For example, the wind probability for 34 kt at a given location from 0 to 120 h is the proportion of the number of realizations in the area covered by the radius of 34 kt during 0–120 h, to the total number of realizations.

Because MCP and AMCP are based on official forecasts (e.g., GFS) with random samples, a poor track forecast may introduce large errors to wind speed realizations. Therefore, additional wind speed probabilistic forecasts are needed to supply MCP or AMCP forecasts and generate a combined wind speed probability forecast. Given the availability of ensemble forecasting at large operational centers and its usefulness in improving operational TC forecasts, it seems natural to combine AMCP with ensemble members for improved forecasts. In this study, we attempt to evaluate whether including ensemble probability will improve MCP forecasts.

Major operational centers make probability forecasts through their ensemble forecast systems with multiple forecasts from numerical weather prediction models with either different initial conditions or varying numerical representations of the model physics (Toth and Kalnay, 1997; Palmer et al., 2009; Barkmeijer et al., 2013; Buizza, 2014). In this study, we use wind speed probability maps from ensemble forecasts, obtained from The International Grand Global Ensemble (TIGGE; https://confluence.ecmwf.int/display/TIGGE), for three major operational centers: ECMWF (with 50 ensemble members), UKMO (with 23 members for Hurricane Sandy and 17 members for Hurricanes Irma and Harvey), and NCEP (with 20 ensemble members). While these ensemble forecasts of TCs provide the uncertainties in the forecasts (Swinbank et al., 2016), we can also make probability maps by counting the occurrence percentage of a certain forecast variable among all ensemble forecasts. The ensemble wind speed probabilities are calculated similar to the method used in the AMCP model, namely, by counting the total number of ensemble members for each of the forecasts that exceed the prescribed threshold (20–64 kt) and dividing by the total number of ensemble forecasts.

By combining the probability maps from both AMCP and the ensemble forecasts as in the following, a modified AMCP (AMCP_mf) is produced:

where P_ECMWF, P_UKMO, and P_NCEP are wind speed probabilities of the ECMWF, UKMO, and NCEP ensemble forecasts, respectively. With TIGGE, the ECMWF, UKMO, and NCEP ensemble forecast data are all available at a grid size of 0.5° × 0.5°. Therefore, the AMCP_mf probability result is also calculated at a grid size of 0.5° × 0.5° for easy processing in Eq. (7).

As shown in Section 3, AMCP overestimates wind speed probability over land (Splitt et al., 2010). Based on tests of various weights for the ensemble and AMCP combinations shown in Table 3, we use a larger weight coefficient (0.9, which leads to a multiplicative bias that is closer to 1) for the ensemble forecast probability results in Eq. (7). Over the ocean, we use weight coefficients of 0.5 for both AMCP and ensemble forecasts to compensate the confidence for each. For consistence, the NCEP GFS final analysis (FNL) data, which are publicly available at a grid size of 0.5° × 0.5°, are used to obtain the analysis (or “truth”) of wind speed frequency by counting the frequency of wind speeds greater than a certain threshold of interest at a given location for the forecast time period. For instance, the frequency of wind speed of 34 kt for a given location is the number of impacted times with wind speed over 34 kt divided by the total number of analysis times of FNL. Namely, for a certain time that has wind over 34 kt over a given area, the frequency should be 1/60 (120-h analysis interpolated into 2-h intervals). The cyclone center is determined by the minimum surface pressure in both ensemble products and NCEP FNL. The ensemble mean track is generated by averaging the cyclone centers from all members. The AMCP_mf track is the average of the AMCP track and the mean ensemble track. Similar to the AMCP_mf model, the track, intensity, and radii from ensemble forecasts and FNL analysis are also linearly interpolated at 2-h intervals from their 6-h data.

Because the change in forecast wind speed probability from MCP to AMCP is small (commonly less than 2%), we emphasize the comparison between AMCP and AMCP_mf in the evaluation. Statistical methods are applied. First, a multiplicative bias is assessed:

where F_i is the forecasted probability, O_i is the frequency from the FNL analysis, and N is the total number of forecasted/observed points over the entire domain.

To assess the accuracy of the above binary deterministic forecast, namely, whether the specific wind speed (e.g., 20 or 34 kt) of a hurricane does or does not affect a certain area based on the AMCP/AMCP_mf probability forecasts, the threat score (TS) is calculated:

where H is the number of hits (forecasted and observed), M is the number of misses (not forecasted but observed), and FA is the number of false alarms (forecasted but not observed). To determine the TS, a threat probability (TP) is defined based on the probability map of whether an event is forecasted with accuracy varying from 1% to 100%. The forecasted event is defined as a location with the wind speed probability over the given TP. Otherwise, it is not a forecasted event. With this definition of TP and observations, the hits, misses, and false alarms are generated. From the perspective of the TS, a larger maximum implies a better deterministic forecast. Thus, TS and TP are used to assess the improvement in AMCP_mf. To assess the coverage of the wind speed warning map, the missing score (MS) and false score (FS) are calculated:

The increased MS (FS) indicates an increasing number of missed alarms (false alarms). Since missed alarms and false alarms are also critical for public disaster warnings, the MS and FS are used to assess the improvement in AMCP_mf from ensemble forecasts, which likely lead to large MS because of the limited ensemble members. Furthermore, a bootstrap method (Efron and Tibshirani, 1994) is used to evaluate the significance of the modification from AMCP to AMCP_mf.

In this section, probability products for Hurricanes Irma, Harvey, and Sandy are discussed to characterize the performance of AMCP_mf compared with AMCP. Figures 1a–d illustrate the 34-kt cumulative wind speed probability forecast generated by AMCP, mean ensemble forecasts (global average forecasts from ECMWF, UKMO, and NCEP probability forecasts), AMCP_mf, and FNL analysis, respectively, for Hurricane Irma for 120-h forecasts at 0000 UTC 9 September 2017. Figures 1e–g show the corresponding 120-h and 34-kt probability differences between the forecasts and FNL analysis cumulative wind speed probability map for Hurricane Irma at 0000 UTC 9 September 2017.

Figure 1.  The 34-kt cumulative wind speed probability for Hurricane Irma from (a) AMCP, (b) the ensemble mean (ENS), (c) AMCP_mf, and (d) FNL analysis for 120 h at 0000 UTC 9 September 2017. (e–g) The 34-kt cumulative wind speed probability difference for Hurricane Irma between AMCP, the mean ensemble forecast, and AMCP_mf, respectively, and FNL analysis for the 120-h forecast at 0000 UTC 9 September 2017. The red, green, magenta, cyan, and black lines (points) represent the hurricane track (positions every 12 h; there are only 10 position points for the best track because of the dissipation of Irma) from GFS, AMCP, the mean ensemble model, AMCP_mf, and NHC best track, respectively.

Figure 2.  As in Fig. 1, but for 27-kt cumulative wind speed probability and probability difference for Hurricane Harvey during 120 h at 0000 UTC 27 August 2017.

From the probability map, there is a significant decrease in wind speed probability in AMCP_mf inland that is closer to FNL, compared with AMCP. The results imply that the combination with the ensemble forecasts will improve wind speed probability forecasts, especially for TCs inland. To qualify this improvement, Figs. 1e–g also give the probability differences between the forecasts and FNL. Since there are areas with zero probabilities for the forecast maps from each probability forecast model, directly subtracting the FNL analysis from these probability maps will generate a lot of invalid points for the assessment. Therefore, these points (or locations) with zero probabilities in all maps, including AMCP, the ensemble, and AMCP_mf, are ignored during the comparison in this section. After eliminating these zero points, the mean (standard deviation) probability difference is 1.34% (3.32%) for AMCP, 0.72% (2.26%) for the ensemble, and 0.66% (1.94%) for AMCP_mf. From the perspective of the mean and standard deviation probability difference, a mean and standard deviation probability difference close to zero implies a small forecast bias. It is obvious that after modification, AMCP_mf provides a better wind probability forecast with minimum bias. For forecasts over the ocean, the mean (standard deviation) probability difference is 0.37% (1.96%) for AMCP, 0.66% (2.07%) for the ensemble, and 0.51% (1.72%) for AMCP_mf. For inland forecasts, the mean (standard deviation) probability difference is 0.98% (2.82%) for AMCP, 0.06% (0.95%) for the ensemble, and 0.15% (0.98%) for AMCP_mf. With a lower mean probability difference, we found that AMCP provides a better probability forecast over the ocean compared to the ensemble. For the probabilities inland, AMCP overestimates the probabilities despite the wind speed decay model contained in the AMCP model (Kaplan and DeMaria, 1995), and the ensemble provides a better forecast inland. During this forecast period, a single model, either AMCP or the ensemble, is not good enough. In contrast, in combination they provide better probability forecasts compared to a single model.

Figures 2a–c illustrate the 27-kt cumulative wind speed probability forecast generated by AMCP, mean ensemble forecasts (global average forecasts from ECMWF, UKMO, and NCEP probability forecasts), and AMCP_mf, respectively, for Hurricane Harvey for 120-h forecasts at 0000 UTC 27 August 2017. Figure 2d shows the corresponding 120-h and 27-kt FNL analysis cumulative wind speed probability map for Hurricane Harvey at 0000 UTC 27 August 2017. Similarly, the false forecast of AMCP inland is significantly reduced after combined with the ensemble forecast. To qualify this improvement, Figs. 2e–g show the probability difference between the forecast and FNL. Similar to the analysis with Hurricane Irma, after ignoring the zero probability points as discussed above, the mean (standard deviation) probability difference is 5.73% (11.67%) for AMCP, 0.17% (1.27%) for the ensemble, and 1.18% (3.24%) for AMCP_mf. While the GFS track moves northwest inland and the best track moves northeast, AMCP produces a poor track realization and a poor wind speed probability forecast inland with a large bias. After modification, AMCP_mf provides a better wind probability forecast with a lower bias compared to AMCP. This indicates that combining AMCP with the ensemble forecasts improves the forecast ability for Hurricane Harvey. For forecasts over the ocean, the mean (standard deviation) probability difference is 1.30% (5.26%) for AMCP, 0.15% (1.24%) for the ensemble, and 0.72% (3.12%) for AMCP_mf. For inland forecasts, the mean (standard deviation) probability difference is 4.43% (10.96%) for AMCP, 0.01% (0.27%) for the ensemble, and 0.46% (1.18%) for AMCP_mf. We found that the ensemble can improve not only inland forecasts but also forecasts over the ocean. Considering that there is a larger error for the GFS forecast track, the ensemble model can modify the AMCP forecast based on that poor track from the GFS.

Figures 3a–c show the 34-kt cumulative wind speed probability maps for Hurricane Sandy from AMCP, mean ensemble forecasts, and AMCP_mf, respectively, for the 108-h forecast at 0000 UTC 28 October 2012. Figure 3d shows the 108-h and 34-kt FNL analysis cumulative wind speed probability map for Hurricane Sandy at 0000 UTC 28 October 2012. We found that the combination not only eliminates the false inland forecast of AMCP but also improves the forecast over the ocean, as AMCP_mf provides better location forecasts of maximum probability. Figures 3e–g show the differences between the forecasts and FNL. After ignoring the zero points, the mean (standard deviation) probability difference is −3.29% (8.54%) for AMCP, 1.75% (4.36%) for the ensemble, and −0.97% (4.77%) for AMCP_mf. For forecasts over the ocean, the mean (standard deviation) probability difference is −3.96% (7.90%) for AMCP, 1.59% (4.36%) for the ensemble, and −1.19% (4.65%) for AMCP_mf. For inland forecasts, the mean (standard deviation) probability difference is 0.67% (2.32%) for AMCP, 0.17% (0.78%) for the ensemble, and 0.22% (0.80%) for AMCP_mf. Compared to AMCP, similar to Hurricanes Irma and Harvey, the AMCP_mf probability forecast is also improved by the ensemble forecast with lower bias. Meanwhile, similar to Fig. 1, AMCP_mf provides the best wind probability forecast with minimum bias compared to AMCP and the ensemble. This again confirms that this combination (AMCP_mf) performs better than either single model, namely, AMCP or the ensemble, in some cases.

Figure 3.  As in Fig. 1, but for 34-kt cumulative wind speed probability and probability difference for Hurricane Sandy during 108 h at 0000 UTC 28 October 2012.

Based on the probability maps above, we can conclude that the AMCP model may produce poor forecasts with a poor essential cyclone track. The ensemble forecast results are useful in modifying the AMCP results, not only for track direction but also for inland probability distribution. To further explore the characteristics of the AMCP, AMCP_mf, and ensemble forecast wind speed probability, Fig. 4 shows three scatter diagrams that indicate the relationship of AMCP, ensemble forecasts, and AMCP_mf to the FNL analysis for the 34-kt and 120-h (108 h for Sandy) cumulative probability forecasts for Hurricanes Irma, Harvey, and Sandy. Considering the small intercept (< 1%), we focus only on the slope of linear fitting in this section. From the perspective of the slope and the determination coefficient (R²) in fitting, a slope close to 1 implies a better forecast and a larger R² implies a small dispersion for the wind speed probability forecast. For Hurricane Irma, the slope (R²) for AMCP, the ensemble, and AMCP_mf is 0.83 (0.75), 1.66 (0.89), and 1.29 (0.89), respectively. Compared to the FNL analysis, AMCP overestimates while the ensemble underestimates the wind speed probability. After modification, AMCP_mf produces a more concentrated probability distribution with increased R². For Hurricane Harvey, the slope (R²) for AMCP, the ensemble, and AMCP_mf is 0.14 (0.22), 0.74 (0.63), and 0.5 (0.53), respectively. Compared to the FNL analysis, both AMCP and the ensemble underestimate the wind speed probability. Similar to Hurricane Irma, after modification, AMCP_mf produces a more concentrated probability distribution with increased R² and a lower bias with slope closer to 1. For Hurricane Sandy, the slope (R²) for AMCP, the ensemble, and AMCP_mf is 1.47 (0.70), 0.96 (0.94), and 1.27 (0.93), respectively. Again, combined with the ensemble, AMCP_mf provides a more concentrated probability forecast and is more likely to distribute along the “x = y” line. Through the above analysis, it can be concluded that combining AMCP with the ensemble can decrease the dispersion of the AMCP probability forecast. In some cases, AMCP_mf is better than the ensemble forecast, proving that the combination of AMCP and ensemble forecasts can provide a better forecast than any of these single models. To further clarify the reliability of AMCP_mf, statistical evaluation is provided in the next section.

Figure 4.  Scatter diagram of AMCP (red points and lines), ensemble (black points and lines), and AMCP_mf (blue points and lines) forecasts of 34-kt cumulative wind speed probability contrasted with the frequency from the FNL analysis for Hurricanes (a) Irma, (b) Harvey, and (c) Sandy. The fitting functions are also added.

In this section, AMCP and AMCP_mf are quantitatively compared to assess the improvement by introducing the ensemble model. First, the multiplicative bias based on cumulative wind speed probability is examined as a function of lead time, as shown in Fig. 5. According to Eq. (8), multiplicative bias > 1 (< 1) implies an overestimation (underestimation) of the averaged forecast probability. Figure 5a shows the averaged multiplicative biases of AMCP and AMCP_mf for Hurricane Irma at different forecast cycles (0000 UTC 8 September 2017 for the 0–120-h forecast; 0000 UTC 9 September 2017 for the 0–108-h forecast; and 1200 UTC 10 September 2017 for the 0–72-h forecast; the forecast times are limited by hurricane dissipation time from the GFS forecast). According to Fig. 5a, the differences between AMCP and the analysis are significantly reduced by adding the ensemble forecasts. After modification, for the 36–72-h forecast, the multiplicative bias for each wind speed threshold is closer to 1, meaning that AMCP_mf produces a wind speed probability closer to the analysis. During this forecast period, the hurricane is in its landfall, and the wind speed in the inland part of the hurricane does not decrease. This causes AMCP to overestimate probabilities for inland areas. At the same time, the ensemble forecast gives a more accurate wind speed and probability forecast. Therefore, the combination with the ensemble can provide better probability forecasts with reduced multiplicative bias. The significant decrease in multiplicative bias for AMCP at the 84–120-h forecast might be caused by wind speed decay through Eq. (3) after hurricane landfall (based on the GFS track, Hurricane Irma’s landfall occurred at 0000 UTC 11 September 2017). The decreased wind speed leads to smaller probability and reduces the difference between AMCP and the analysis, especially for a wind speed of 64 kt, in which the multiplicative bias decreases from 2.57 (for the 72-h forecast) to 1.36 (for the 84-h forecast).

Figure 5.  Averaged multiplicative biases of AMCP (bars) and AMCP_mf (black points) at different wind speeds for Hurricanes (a) Irma, (b) Harvey, and (c) Sandy based on cumulative probability. The error bars for AMCP and AMCP_mf are shown as gray and black whiskers, respectively.

For Hurricane Harvey, the averaged multiplicative biases of AMCP and AMCP_mf for different forecasts (0000 UTC 25–28 August 2017 for the 0–120-h forecast, and 29 August 2017 for the 0–108-h forecast; the forecast times are limited by the hurricane dissipation time from the GFS forecast) are shown in Fig. 5b. Because of the lower wind speeds of Hurricane Harvey at landfall (see Fig. 2d), only the multiplicative biases for AMCP and AMCP_mf at 20–34 kt are displayed here. Errors in both the GFS forecast track and inland AMCP intensity realizations produce higher probabilities inland compared with the FNL analysis. These higher probabilities from the AMCP forecast lead to a very large multiplicative bias, with a maximum of about 9.66 for 20-kt wind, 26.21 for 27-kt wind, and 103.89 for 34-kt wind. After modification, all multiplicative bias decreases, and the largest multiplicative bias at 20, 27, and 34 kt is reduced to about 2.76, 6.25, and 23.00, respectively.

For Hurricane Sandy, the averaged multiplicative biases of AMCP and AMCP_mf for different forecast cycles (0000 UTC 27 October for the 0–108-h forecast; 28 October for the 0–84-h forecast; and 29 October 2012 for the 0–60-h forecast; the forecast times are limited by the hurricane dissipation time from the GFS forecast) are shown in Fig. 5c. Considering that there are many lower probability forecasts in AMCP (see Fig. 3a), the multiplicative bias is less than 1 for wind speeds of 20, 27, and 34 kt. The probability at 50 and 64 kt has a multiplicative bias over 1 because of the lower frequency of the FNL high wind speed. That is, AMCP overestimates high wind speeds and leads to large multiplicative bias. After modification including the ensemble forecasts, the multiplicative bias increases at 20, 27, and 34 kt, while it decreases at 50 and 64 kt, leading to a better probability forecast for multiplicative bias close to 1.

Based on the cumulative forecast data depicted in Fig. 5, Fig. 6 shows the maximum TS and corresponding TP variation before and after the modification for the combined datasets of Hurricanes Irma, Harvey, and Sandy. Because the maximum TSs for Hurricanes Irma, Harvey, and Sandy are all higher in the AMCP_mf than in the AMCP forecast, the maximum TS shows an increase of about 0.18–0.33 at 20–34 kt, as shown in Fig. 6a. This indicates an improvement in binary deterministic forecasts for the modified method. Meanwhile, the corresponding TP is lower except at 27 kt for the 12-h forecast. The larger TS and smaller TP indicate a better binary (TC affects or does not affect) forecast map for disaster warnings.

Figure 6.  (a) Maximum threat score (TS) and (b) corresponding threat probability (TP) of AMCP (bars) and AMCP_mf (gray points) at different wind speeds for the combined cumulative forecast sets of Hurricanes Irma, Harvey, and Sandy. The black points represent the significance level from the bootstrap t-test for the TS and TP change between AMCP and AMCP_mf.

To determine the significance of these changes in TS, a bootstrap t-test is used here. The bootstrap t-test is a useful method to assess the difference between two sets of small-sample data, as it randomly samples from the small sample to verify whether or not a change in the small-sample data is significant (Efron and Tibshirani, 1994). At a specific lead time from different forecast cycles, the threat score is first transformed according to Efron and Tibshirani (1994) and then randomly sampled to generate the samples for the bootstrap t-test. The probability of T (test statistic from new samples) larger than t (test statistic from original data) is compared to the significance level to determine the improvement. Through the bootstrap t-test, we found that the maximum TS improvement for forecast wind speed of 20–27 kt is significant at the 0.05 significance level at each forecast time. Some forecasts (at 24 and 48–120 h) at 34 kt produce significant improvements at the 0.15 significance level, and at the 0.05 significance level for some forecast times. The improvement is not significant for wind speeds of 50–64 kt at a significance level greater than 0.3. The maximum TS indicates the skill of the deterministic forecasts, that is, the binary forecast of whether or not the hurricane can affect a specific location. Therefore, the increased maximum TS at 27–34 kt implies the skill of AMCP_mf in deterministic forecasts, namely, the coverage area for the AMCP_mf probability map is closer to that from the FNL analysis. As indicated by the discussion in Section 3.1, the AMCP coverage area is controlled mainly by the basic forecast track. Therefore, the introduction of the ensemble model results in the AMCP model modifies the poor track forecast and improves the wind speed probability deterministic forecast at low wind speeds. For TP, the bootstrap t-test implies that almost all of the change in TP is significant at significance levels less than 0.1, and less than 0.05 except for the 64-kt forecast at 72, 84, and 108 h. Usually, the lower probabilities always distribute along the edge of the probability map. Therefore, the decreased TP implies that a more accurate probability map can be used to produce binary deterministic (TC affects or does not affect) forecasts. Overall, the AMCP_mf model can provide a better deterministic forecast wind speed probability map at 20–34 kt and is more accurate at all wind thresholds.

In this section, 120-h wind speed probability forecasts are conducted for two major landfalling hurricanes over the Atlantic Ocean in the 2018 season. Forecasts at 120 h for Hurricane Florence (first and last realization case initialized at 0000 UTC 10 and 1200 UTC 16 September 2018) and Hurricane Michael (first and last realization case initialized at 0000 UTC 8 and 0000 UTC 11 October 2018) are conducted at 0000 and 1200 UTC every day. The improvements of AMCP_mf over the ensemble forecasts and AMCP model are quantified.

Similarly, the multiplicative bias for major hurricane cases in 2018 based on the cumulative wind speed probability is examined as a function of lead times (Fig. 7). Similar to the case study for Hurricanes Irma, Harvey, and Sandy, the multiplicative bias is significantly reduced by adding the ensemble forecasts. After modification, especially for the 36–120-h forecast, the multiplicative bias for each wind speed threshold is closer to 1, meaning that AMCP_mf produces a wind speed probability closer to the analysis, consisting of the above case study.

Figure 7.  As in Fig. 5, but for major hurricanes in the 2018 season. The dotted line indicates that the multiplicative bias is equal to 1.

Figures 8a and 8b show the TS and TP for the major hurricane cases in the 2018 season based on the cumulative wind speed probability at each forecast lead time. As in the case studies in the previous section, for the major hurricanes in the 2018 season, AMCP_mf provides larger TS (except that of 64-kt wind speed at 120 h) and smaller TP at each lead time. The bootstrap t-test implies that the increased TS is significant at the 0.05 significance level (0.3 significance level) at the 36–120-h forecast for the wind speed threshold of 20–34 kt (50 kt). Simultaneously, the decreased TP at 36–120 h is significant at the 0.05 significance level except for 64-kt wind at 36 and 120 h. Significantly increased TS and/or notably decreased TP indicate a strong enhanced binary deterministic forecast for 20–64-kt wind. Overall, adding the ensemble forecast can provide a more accurate probability map that can be used to produce binary deterministic (TC affects or does not affect) forecasts.

Figure 8.  As in Fig. 6, but for the combined cumulative forecast sets of the major hurricanes in the 2018 season.

According to the case studies, to some extent, the ensemble forecast is possibly more accurate than AMCP and AMCP_mf. However, for a public disaster warning system such as MCP/AMCP/AMCP_mf, missing alarms are more disastrous than false alarms. Considering that the MS and FS of high wind are more critical for public disaster warnings, Fig. 9 shows the false and missing scores (at the TP that generates the largest TS) of the ensemble and AMCP_mf based on the cumulative wind speed probability at high wind (50–64 kt) thresholds. Figure 9 indicates that MS is always larger, while FS is often smaller in the ensemble forecasts. Compared to the strong reduced MS of AMCP_mf, especially for the wind speed of 64 kt, FS often increases slightly. The combination of AMCP and the ensemble forecasts could extend the wind speed probability map (see Figs. 1–3) and reduce MS from the ensemble forecasts. Figure 9 also indicates that the reduced missed alarms will not all become false alarms after the combination, and some will become hits. Therefore, with decreased MS and slightly increased TS, the extended wind speed probability map from AMCP_mf can provide useful hit warning information that the ensemble forecasts miss. Due to the sample size issue, namely, the limited size of ensemble members, the probability map from the ensemble forecasts is often smaller than that from AMCP_mf. The smaller disaster warning map could lead to many missed alarms and thus result in a severe disaster. Therefore, the combination of AMCP and ensemble forecasts would be better than either one alone for providing a public disaster warning (see Fig. 2).

Figure 9.  (a) False score (FS) and (b) missing score (MS) of AMCP_mf (bars) and ensemble (black points) at wind speeds of 50–64 kt for the combined cumulative forecast sets of the 2018 major hurricane season. MS and FS are collected at the threat probability corresponding to the maximum threat score in each model.

In this study, we first made a slight adjustment to the Monte Carlo hurricane wind speed probability model by adding limits to the direction angle and motion distance. Then, we introduced ensemble probability results derived from ECMWF, UKMO, and NCEP ensemble forecasts to this model to improve the predictions based on GFS model forecasts. Hurricanes Irma (2017), Harvey (2017), and Sandy (2012) and the major hurricanes in the 2018 season are used to evaluate the effectiveness of the modification.

The probability map indicates that the ensemble forecasts can improve not only inland forecasts but also forecasts over the ocean, and this combination (AMCP_mf) is better than either single model result. With 1000 random samples, the realizations generate a large number of possible tracks and cover larger areas, and this poor forecast, especially inland, can be corrected by the ensemble forecasts. Meanwhile, the AMCP model can be a good complement to the ensemble forecasts due to the limited sample size of the ensemble members.

The significant improvement in the maximum threat score, as well as the significantly decreased corresponding threat probability, implies that the AMCP_mf model can provide a better deterministic forecast wind speed probability map at 20–64 kt. The scatter diagram shows that AMCP_mf for all cases is close to the “x = y” line, which indicates a better wind speed probability forecast during the 120-h forecasts.

Overall results from this study prove the benefits of combining the MCP model with ensemble forecasting in potential applications for improved TC forecasts. For disaster warnings to the public, the increased accuracy of AMCP_mf could provide accurate wind speed probability forecasts. The increased deterministic forecast ability means more information from AMCP_mf can be used to provide a disaster warning. The larger threat score and smaller threat probability imply a larger map to determine which places TCs can or cannot affect. Although the ensemble is often more accurate than AMCP_mf, missed alarms (not forecasted but observed) are more disastrous than false alarms (forecasted but not observed) for public warnings. Direct application of the ensemble forecasts would miss some early warnings for high wind, which would go against the purpose of the AMCP. Therefore, the combination of MCP and the ensemble could produce better wind speed probability forecasts.

Acknowledgments. Archived hurricane data and ensemble forecast data are obtained from the public website hosted by the University Corporation for Atmospheric Research (UCAR; https://rda.ucar.edu/) and the National Hurricane Center (https://ftp.nhc.noaa.gov/atcf/archive/).