The tropical cyclone best track data used in the present research is from the Joint Typhoon Warning Center, http://www.metoc.navy.mil/jtwc, which provides the latitude, longitude, wind speed, and sea level pressure of tropical cyclones for every 6 h (Chu et al., 2002). Our study period is from 1999 to 2010, the common period of the reforecasts of the WMO S2S Prediction Project Database (Table 1). Since tropical cyclones tend to appear in boreal summer in the western North Pacific, for each year we only analyze the period from May 1 to November 30, with the last forecast starting from the end of October. We also only use those tropical cyclones with sustained 10 min wind speeds greater than 17.2 m/s at least once during their life cycle, which is a standard criterion for separating tropical storms from tropical cyclones in the Indian Ocean and western North Pacific Ocean (Holland 1993).

To describe the intraseasonal variability during boreal summer, especially the variability in the south and east Asia monsoon region, Lee et al. (2013) proposed BSISO index using multi-variable Empirical Orthogonal Function analysis of 850 hPa zonal wind and outgoing longwave radiation. The time series of the first four principal components are used to construct two indices, BSISO1 and BSISO2. The BSISO1 and 2 indices used are directly from the website http://www.iprc.soest.hawaii.edu/users/jylee/bsiso/, computed using NCEP/DOE Reanalysis II (Kanamitsu et al., 2002). BSISO1 and BSISO2 describe different aspects of boreal summer intraseasonal variability. These two indices are associated with different spatial patterns of wind and outgoing longwave radiation and temporal scales. The temporal scale of BSISO1 is around 30–60 days, while the temporal scale of BSISO2 is around 10–30 days.

The details of the eleven S2S forecast systems in the S2S Prediction Project Database are summarized in Table 1 which is based on Vitart et al. (2017). The S2S forecast systems have a typical spatial resolution of .25^o to 2^o degrees. Seven of them are coupled systems with a dynamic oceanic component. The other four systems are not coupled. These systems provide S2S reforecast with lead time varying from 33 to 62 days. In addition, each system has different forecast days and frequencies for their reforecasts, with the number of ensemble members varying from as few as 3 to as many as 33. Only two systems (NCEP, CMA) made reforecasts every day. The number of S2S prediction samples will influence the prediction skill evaluation in subtle ways since tropical cyclones are rare events. The operational year of these systems (last column in Table 1) is also different. In comparing their skill for predicting tropical cyclone activity in the western North Pacific, we used a common grid of 1^o × 1^o, a prediction lead time of 30 days, prediction samples of common operational years from 1999 to 2010.

The output from these systems in the S2S Prediction Project Database was archived on a 1.5^o × 1.5^o grid every 24 h. Previous research has shown that tropical cyclones, or close proxies, can be identified and tracked in the output of these systems and used successfully in the evaluation of the genesis (Lee et al., 2018) and weekly tropical cyclone occurrence for 15^olatitude × 20^o longitude boxes (Lee et al., 2020). The tracking algorithm used was developed by Vitart et al. (1997) and further modified by Vitart and Stockdale (2001). The algorithm was adapted for operational use at ECMWF since 2011 (Vitart et al., 2011). The tracking algorithm defines a tropical cyclone center at a local minimum sea level pressure center where 1) a local vorticity maximum ( > 3.5 × 10 − 5   s − 1 ) at 850 hPa is within 2^o latitude, 2) a warm core with 500–250 hPa vertically averaged temperature greater than .5°C is within 2^o latitude, 3) the two locations detected in (1) and (2) are within a distance of 8^o latitude, and 4) a local maximum thickness center between 1,000 and 200 hPa is within 2^o latitude. The above criteria were tested for an atmospheric general circulation model with a resolution of 2.8^o (Vitart et al., 1997) and the ECMWF Seasonal Forecast System with a resolution of 1.875^o (Vitart and Stockdale, 2001). The number of model tropical cyclones tracked using the above criteria in the above ECMWF systems is consistent with the number of tropical cyclones with wind speed greater than 17 m/s (Vitart and Stockdale, 2001).

Different forecast systems usually have their own tracking algorithm to identify tropical cyclones in their output (e.g., Camp et al., 2018), which may render different results from the tracking algorithm used in our research. However, the output in the S2S Prediction Project Database was already archived on a 1.5° × 1.5° grid, and it is not practical to use different tracking algorithm for different forecast systems since the tracking algorithms used in different systems and specific criteria are resolution-dependent. The DTCP we used to quantify the evolution and movement of tropical cyclones, only makes use of the location information of tropical cyclones. The influence of using different tracking algorithms may not be significant. On the other hand, all the S2S forecasting systems are directly evaluated against the observed tropical cyclones. Postprocessing, which is forecast system dependent and can improve the prediction skill of a forecast system (e.g., Camp et al., 2018; Lee et al., 2020), is not applied for any of the systems.

The DTCP in the present research is defined as the probability of tropical cyclone occurrence within 500 km of a grid point during a particular day. The temporal and spatial criteria of 1 day and 500 km from a specified location are selected for their simplicity in terms of interpretation and application. Different criteria may also be used, like 1 week and 15^olatitude × 20^o longitude (Lee et al., 2020). It will render different (and higher) prediction skills. However, with a narrow temporal window and specific spatial information that is straightforward to interpret, the metric DTCP can pinpoint the genesis and movement of tropical cyclones.

The DTCP is computed on a common 1° × 1° grid over the western North Pacific (0°-30°N, 100°-180°E) for both observation and model outputs. The computation of DTCP is slightly different for ensemble model output and observation. For a particular ensemble forecast system, DTCP is computed for each day as the number of tropical cyclones present in all the ensemble members within 500 km of a location divided by the total number of ensemble members. If this value is greater than 1, it is set as 1. As such, the computed DTCP on each grid is the probability that a tropical cyclone would occur within a radius of 500 km in the ensemble forecast in 1 day. Strictly speaking, the above definition is the average of the number of tropical cyclones capped by one as pointed out by one reviewer. The DTCP from ensemble model output should be the number of members that have tropical cyclones within 500 km of a location divided by the total number of members. Since it is very rare to have two tropical cyclones within 500 km in a model with course resolution of ∼1^o degree, we ignore the difference of the two computations.

For observations, the DTCP is computed as the number of tropical cyclones within 500 km of a location at 0GMT for each day. When the number of tropical cyclones exceeds 1, it is set as one. Different from the definition of tropical cyclone density used in Vitart and Robertson (2018), which is the number of tropical cyclones within 200 km of a location, DTCP defined in the present research is always less than or equal to one, and can be interpreted as the probability a tropical cyclone (or tropical cyclones) would occur within 500 km of a particular location. Following the above definition, the occurrence of a tropical cyclone can be treated as a dichotomous event. The DTCP from an ensemble forecast and observation can be directly compared and the forecast skill of tropical cyclones can be measured by debiased Brier Skill Score (Weigel et al., 2007) which can account for the difference in the number of ensemble members compared to the Brier Skill Score proposed by Brier (1950).

As an example, Figure 1 presents DTCPs from eleven forecast systems (Figures 1A–K) and from the observation (Figure 1L) on 20 May 2004. Forecast lead times for this case are selected from 11 to 20 days to focus on the S2S timescale and simultaneously have a minimum difference in lead time among models. On this day, two tropical cyclones are observed in the western North Pacific (Figure 1L), located in the east of the Philippines, one around 12^oN, the other around 21^oN. In observation (Figure 1L), the DTCP is 1 around the neighborhood of these two tropical cyclones. The eleven forecast models exhibit varying skills in predicting these two tropical cyclones (Figures 1A–K), partly related to different lead times and partly related to their skills in general. For this case, several models, except JMA and HMCR, predict the existence of at least one tropical cyclone. However, these models have various errors in the location of these two tropical cyclones. It should be noted that some of the forecast models tend to overpredict the number of tropical cyclones (ECCC, METFR, NCEP, UKMO). Overall, the forecast results are encouraging and it is plausible to directly compare the DTCP from observation and forecast models to assess the prediction skill of tropical cyclones in these models on the S2S timescale. If the metric to evaluate the skill of tropical cyclone prediction is the number of tropical cyclones within a region (say, 110E-150E, 0-30N), eight of these 11 systems predict one or more tropical cyclones within the region and would have some skills based on this metric. However, the location information of the tropical cyclone might be too broad to be useful for societal applications (Camargo et al., 2019).

With the above definition of DTCP, the daily occurrence of tropical cyclones can be treated as dichotomous events with 1 if the event (tropical cyclone) occurs and 0 if it does not. The Brier Score can then be used as a measure of the accuracy of the tropical cyclone prediction (Brier 1950), B S = ∑ k = 1 2 ∑ j = 1 N   Y k j − O k j 2 / N (1) where · is the ensemble mean over prediction samples (j = 1,…, N) and k = 1, 2 is the number of categories, Y k j is the prediction for category k and sample j. O k j is the observed DTCP for sample j. The value of B S changes from 0 to 2, with 0 as the perfect prediction and 2 as the situation when the prediction is wrong every time.

In order to assess the skill of a probability forecast, the Brier Skill Score compares the forecast with a reference forecast, such as a climatological forecast, and is defined as 1 − B S B S C L (2) in which B S C L is the Brier Score when a reference climatological forecast is used to forecast DTCP. If the climatological DTCP is p for one location, then the B S C L is 2p(1-p) in theory. In our evaluation of the eleven forecast models, the theoretical B S C L computed from climatological p as shown in Figure 2L is used instead of its estimate from individual model prediction samples.

When ensemble forecasts are involved, small number of ensembles tends to cause a negative bias in the value of the Brier Skill Score as defined by Eq. 2. To remove the bias, a debiased Brier Skill Score (   B S S D ) was proposed to evaluate the categorical event forecast (Weigel et al., 2007). In our study,   B S S D used to evaluate the tropical cyclone probability forecast is B S S D = 1 − < B S > < B S C L > + D (3)

In Eq. 3, D is D = 1 M p 1 − p (4) in which M is the number of ensembles and p is the observed climatological DTCP defined above. The term D is used to correct the negative bias when the number of ensemble members is small. As evident from Eq. 4, the correction term D is reciprocal to the number of ensemble members, which becomes small for a large ensemble size. The debiased Brier Skill Score B S S D is used to compare the forecast skill of DTCP of the eleven forecast models. The B S S D is computed for each grid point, which is different from the practice such as in Yamaguchi et al. (2015), in which the N in Eq. 1 is for grid points in a region not for prediction samples. The area average B S S D is used as a measure of tropical cyclone forecast for a specific region. Though the B S S D is computed for the northwestern Pacific (0°-30°N,100°-180°E), the area averaging is conducted only for the oceanic points of the region in northern South China Sea and east of the Philippines (10°-30°N,105°-150°E), where tropical cyclones tend to appear more frequently.

To quantify the skill of the S2S models for reproducing the climatological DTCP and its modulation by BSISO, the Taylor Score (TS) index proposed by Taylor (2001) was used, T S = 4 1 + R 4 σ m σ o + σ o σ m 2 1 + R 0 4   , (5) where R is the pattern correlation coefficient between the model results and observations, and R ₀ is the maximum attainable pattern correlation coefficient between the model results and observations, which was taken as 1 in the present research. The pattern correlation coefficient of DTCP was computed using DTCP from each oceanic grid point. The variables σ m and σ o are the standard deviations of the modeled and observed spatial fields, respectively. The values of the TS index were between 0 and 1. A larger value for this index indicates a more similar model spatial pattern to the observations. The TS indices for models in reproducing the climatological DTCP and its modulation by BSISO1 and BSISO2 are shown in Table 2 and used in Section 3.5.

In order to assess the tropical cyclone prediction in these eleven models, we first evaluate their skill in reproducing the climatological DTCP from 1999 to 2010. To have a robust comparison among eleven models, the days from May 1 to October 30 of each year from 1999 to 2010 are all filled with the closest reforecasts. Thus, all the eleven forecast models have the same number of reforecast of DTCP to compare with observations. Figure 2 compares the climatological DTCP from these forecast models (Figures 2A–K) and the observation (Figure 2L). As expected, a large DTCP appears in the northern part of the South China Sea and east of the Philippines. These forecast models can more or less capture the spatial distribution of the observed DTCP, although fidelity varies. Most of these models (Figures 2D, F, G, I–K) can reproduce the DTCP east of the Philippines, but underestimate DTCP in the northern part of the South China Sea. The maximum observed DTCP is around .1 (Figure 2L), which means one would expect one tropical cyclone every 10 days.

To qualitatively compare the skill of these models in reproducing the observed climatological DTCP, the Taylor diagram is used following Taylor (2001). The ECMWF model can represent the spatial distribution of climatological DTCP the best, followed by NCEP, KMA. and JMA, as shown in the Taylor diagram (Figure 3) and Table 2. The capability of these models in reproducing climatological DTCP will influence their forecast skills.

Tropical cyclone variability in the western North Pacific is modulated by intraseasonal oscillation (e.g., Vitart 2009; Zhao et al., 2015; Nakano et al., 2021). To evaluate the performance of the eleven forecast models in reproducing the modulation of tropical cyclone variability by intraseasonal oscillation, the DTCP composite for eight phases of BSISO1 and BSISO2 are compared with that from the observation. The same method used in Section 3.1 is also used here to have a 2184-day output from 1999 to 2010 between May 1 and October 31. It should be noted that the composite for model output is based on the observed BSISO index computed by Lee et al. (2013) using NCEP/Department of Energy Reanalysis Ⅱ (Kanamitsu et al., 2002). When constructing composite for BSISO phases, only the days when the magnitude of the BSISO index is greater than 1.5 are used as in Lee et al. (2013).

Modulation of the DTCP by the BSISO is evident in observation (bottom row of Figure 4). It is obvious that the DTCP is much more prominent when BSISO1 is in phases 5, 6, 7, 8, and 1 than that in phases 2, 3, and 4. Eleven models can reproduce this modulation feature with various fidelity. Based on the Outgoing Longwave Radiation (OLR) and 850 hPa wind anomaly patterns of different phases of BSISO1 (Lee et al., 2013), it is evident that the enhanced DTCP is associated with active convection in the northern part of the South China Sea and western North Pacific and the low-level cyclonic circulation that gradually moves northeastward when BSISO1 is in phases 5, 6, 7, 8, and 1 (Figure 9 in Lee et al., 2013). From the observed DTCP, the northeastward movement of large DTCP fields is quite apparent for BSISO1 phases 5, 6, 7, 8, and 1. The observed average DTCP over the northern part of the South China Sea and east of the Philippines (10°-30°N, 105°-150°E) is ∼3.5 times larger during the active phases of the BSISO1 (phases 5, 6, 7, 8, and 1) than that during the non-active phases (phases 2, 3, and 4).

Similarly, BSISO2 can also modulate the activity of tropical cyclones. Figure 5 illustrates the composite DTCP during the eight phases of BSISO2. When BSISO2 is in phases of 8, 1, 2, 3, and 4, the DTCP is enhanced in the northern part of the South China Sea and north of Taiwan (the bottom row of Figure 5). The average DTCP of the northern part of the South China Sea and east of the Philippines (10°-30°N, 105°-150°E) is ∼2.5 times larger during the phases 8, 1, 2, 3, and 4 than that during the phases 5, 6, and 7. The modulation of DTCP by BSISO, which would have been reproduced in forecasting models, would certainly provide a source of predictability of tropical cyclone activity at the S2S timescale.

To qualitatively compare the skill of these models in reproducing the modulation of tropical cyclone activity of BSISO, the DTCP fields for eight phases of BSISO over the northern part of the South China Sea and east of the Philippines (10°-30°N,105°-150°E) are combined to compare with the observation in Taylor diagrams (Figure 6A for BSISO1; Figure 6B for BSISO2). From Figures 6A, B, it is evident that the ECMWF model stands out as the best one to reproduce the modulation of tropical cyclone activity by BSISO, followed by NCEP, KMA and other systems. ECMWF has the highest Taylor Score in reproducing the climatological DTCP and its modulation by BSISO1 and BSISO2 (Table 2). The relationship between the model performance in representing the climatological DTCP and the modulation of tropical cyclone activity by BSISO and the skill of S2S tropical cyclone forecast is discussed in Section 3.5 using these Taylor Score indices.

From Figure 2, the climatological DTCP has large values around the northern part of the South China Sea and east of the Philippines (10°-30°N,105°-150°E). To evaluate the performance of these models in terms of their skills in forecasting tropical cyclones, the debiased Brier Skill Score B S S D averaged for the region is shown in Figure 7. With the increase of lead time from 1 to 30 days, the B S S D progressively drops to zero with a time scale of 3–12 days with the longest one from the ECMWF model. The averaged B S S D from 11 to 30 days is below zero, as shown in Table 2 for these eleven models since our evaluation criteria is strict compared with the one used in Lee et al. (2020). Among the eleven models, ECMWF has the highest averaged B S S D , −.002, from 11 to 30 days. A negative B S S D indicates that the skill of the ECMWF model is slightly worse than the reference climatological forecast. The skill score of ECMWF is followed by HMCR, BoM, METFR, and JMA (Table 2). The HMCR is a special case since it has the lowest Taylor Score for reproducing climatological DTCP and its modulation by BSISO1 and BSISO2. From Figure 2, the HMCR model severely underestimates climatological DTCP and the model makes predictions occasionally (Table 1). The case of HMCR shows that if a model does not predict tropical cyclones at all and only makes prediction occasionally, it might achieve a fake high skill. By investigating the climatological DTCP of the model, this type of forecast models should be easily picked up and not be confused with a model that actually has skill. This practice might be applicable when evaluating the forecast skill of other rare extreme events at S2S timescale.

Regarding the performance of ECMWF, the results above can be compared with that of Lee et al. (2018), Lee et al. (2020). Using the same S2S Prediction Project Database but only for reforecasts from six models, Lee et al. (2018) shows that the ECMWF model has the best performance in predicting tropical cyclone genesis at weekly time windows for all the tropical cyclone basins. Using weekly tropical cyclone occurrence for 15^o latitude × 20^o longitude boxes, Lee et al. (2020) shows that the ECMWF model has the highest prediction skill in all the tropical cyclone basins. The temporal window and spatial distance selection can impact significantly prediction evaluation results.

In terms of the spatial distribution of B S S D , the eleven forecasting models exhibit different regions where their forecast skill is better than that of the reference climatological forecast. Figure 8 presents the   B S S D distribution at a lead time of 17 days for all the eleven forecasting models, which shows the typical spatial feature of   B S S D distribution and presumably excludes the influence of pre-existing tropical cyclones because it is very rare for a typical tropical cyclone to last 17 days in the open ocean after its formation. The regions where negative   B S S D exist, indicating no forecasting skill relative to the reference climatological forecast, are different for each of these models. If we compare Figure 8 with Figure 2, the regions where a forecast system has negative   B S S D are also the regions where the climatological DTCP of the system shows large discrepancies compared with the observation. The CMA system (Figure 2B) shows a systematic bias in reproducing the observed DTCP (Figure 2L). The climatological DTCP of the CMA system is located south of 18°N, which is systematically shifted southward compared with that of the observation. The negative   B S S D region of the CMA system is also located south of 20°N (Figure 9B). Compared with the observation, the NCEP forecast system tends to underestimate the DTCP east of 130°E (Figure 2J vs. Figure 2L), the negative   B S S D of the NCEP system also appears in this region. The ISAC and JMA systems tend to underestimate the DTCP in the northern part of the South China Sea (Figures 2F, G). The negative   B S S D regions of these two systems are also present in the northern part of the South China Sea. Interestingly, the positive   B S S D regions of the HMCR system are places where tropical cyclones do not occur with a considerable probability, indicating that a system which does not forecast extreme events at all may achieve a higher score than a reference climatological forecast if it makes no-event forecast each time. The   B S S D of ECMWF system generally has positive values in the northern part of the South China Sea, around Taiwan, and east of the Philippines, where tropical cyclones often appear. However, some of its skill is from the region east of 140°E, where tropical cyclones do not occur very often.

The modulation of DTCP by BSISO can also influence its forecast skill. Here we consider the DTCP forecast with lead time 11–30 days and compare the skill between the active phases of tropical cyclones (phases 8, 1, 2, 3, and 4) and non-active phases (phases 5, 6, and 7) for BSISO1. We specifically examine the influence in the ECMWF model given that it demonstrates better forecast skills than the reference climatological forecast for lead time longer than 10 days. Figure 9 compares the B S S D for active and non-active phases of BSISO1 when the lead time is 11–30 days to exclude the influence of pre-existing tropical cyclones. Significant differences in the spatial distribution of B S S D are readily discerned. The ECMWF model actually has better skill for the non-active phases of tropical cyclones on the S2S timescale. In the active phases, the region that exhibits better skill is located in central Philippines. The averaged   B S S D for the active phases is negative, while the averaged   B S S D for the non-active phases is positive. The averaged   B S S D for the active and non-active tropical cyclone phases of the other models also exhibit higher skill for non-active phases (Table 2). The averaged   B S S D is higher during the non-active tropical cyclone phases of BSISO1. Lee et al. (2020) have shown that the prediction skill is higher when the predictions start from the active phases, who were concerned with the existence and movement of tropical cyclones inside a basin rather than its presence and movement on a more finer spatial and temporal scale. The above analysis indicates that the current S2S models have better skills for situation of low tropical cyclone probability.

To understand which factors can better explain the S2S prediction skill difference among models, we compute the correlation coefficient between the averaged B S S D of the northern part of the South China Sea and east of the Philippines (10°-30°N,105°-150°E) and the Taylor Score indices of models for reproducing the climatological DTCP and its modulation by BSISO1 and BSISO2 (Figure 10). Since the HMCR model is a particular case, this model is not included in the correlation coefficient computation. The JMA model has B S S D starting from lead time day 2. Thus, the correlation coefficient calculation starts from lead time day 3. The significant correlation coefficient at 90% level is also labeled with a red circle. In general, the correlation coefficients between the model’s capability in reproducing the climatological DTCP and its modulation by BSISO1 and BSISO2 are significant at lead time day 3 and decrease as the lead time increases. Before lead time day 10, the correlation coefficients between Taylor Score indices of BSISO1 and B S S D and Taylor Score indices of BSISO2 and B S S D are below zero. However, the correlation coefficient between Taylor Score indices for reproducing climatological DTCP and B S S D is always positive until lead time day 30. This might indicate that the model’s capability in reproducing climatological DTCP can better explain the B S S D difference among models.