The Impact of Stochastic Physics-Based Hybrid GSI/EnKF Data Assimilation on Hurricane Forecasts Using EMC Operational Hurricane Modeling System

Abstract

The National Oceanic and Atmospheric Administration’s (NOAA) cloud-permitting high-resolution operational Hurricane Weather and Research Forecasting (HWRF) model includes the sophisticated hybrid grid-point statistical interpolation (GSI) and Ensemble Kalman Filter (EnKF) data assimilation (DA) system, which allows assimilating high-resolution aircraft observations in tropical cyclone (TC) inner core regions. In the operational HWRF DA system, the flow-dependent background error covariance matrix is calculated from the HWRF self-cycled 40-member ensemble. This DA system has proved to provide improved initial TC structure and therefore improved TC track and intensity forecasts. However, the uncertainties from the model physics are not taken into account in the FY2017 version of the HWRF DA system. In order to further improve the HWRF DA system, the stochastic physics perturbations are introduced in the HWRF DA, including the cumulus convection scheme, the planetary boundary layer (PBL) scheme, and model surface physics (drag coefficient), for HWRF-based ensembles. This study shows that both TC initial conditions and TC track and intensity forecast skills are improved by adding stochastic model physics in the HWRF self-cycled DA system. It was found that the improvements in the TC initial conditions and forecasts are the results of ensemble spread increases which realistically represent the model background error covariance matrix in HWRF DA. For all 2016 Atlantic storms, the TC track and intensity forecast skills are improved by about ~3% and 6%, respectively, compared to the control experiment. The case study shows that the stochastic physics in HWRF DA is especially helpful for those TCs that have inner-core high-resolution aircraft observations, such as tail Doppler radar (TDR) data.

1. Introduction

Ensemble forecasting has been widely used to take account of the uncertainties in both model initial conditions and model dynamics and physics. Because of the chaotic nature of the atmosphere, a small difference in numerical weather prediction (NWP) model initial states can result in very different model forecasts (Lorenz, 1963, 1965, [1,2]). The source of uncertainties in model initial conditions include the imperfect observation network, observational errors, and model resolutions. On the other hand, highly non-linear model physics, and interactions among the large-scale and sub-grid scale can also cause divergent solutions (Palmer, 2001, [3]). The uncertainties in model dynamics and physics include resolution truncation errors and imperfect model physics parameterization schemes. There are two main applications of ensemble prediction system (EPS) in an NWP model system, predicting the probability density function (PDF) of the model prognostic variables due to model initial and physics uncertainties, and providing the background error covariance matrix for the data assimilation (DA) system. Both applications have been widely utilized in tropical cyclone (TC) NWP modeling systems. While many studies have focused on using uncertainty information from EPS to improve hurricane track and intensity forecasts (Zhang and Krishnamurti 1997, 1999, [4,5]; Krishnamurti 2000 [6]; 2013 Weber [7], Zhang et al., 2014 [8]), there have been fewer studies on the impact of EPS on the TC track and intensity, through taking account of the uncertainties in model physics, which is supposed to provide a more realistic background error covariance matrix for TC inner-core DA.

Inner-core data assimilation (DA) plays an important role in cloud-resolving, high resolution hurricane model forecasts (Pu et al., 2009 [9]; Zhang et al., 2011 [10], Weng and Zhang, 2012 [11], Li et al., 2012 [12], Zhu et al., 2016 [13], Christophersen et al., 2017 [14]). In the National Centers for Environmental Prediction (NCEP), a self-cycled hybrid Ensemble Kalman Filter (EnKF) variational DA system, which uses the Hurricane Weather and Research Forecasting (HWRF)-based ensemble rather than the Global Data Assimilation System (GDAS)-based ensemble to provide the forecast error background covariance matrix, was developed and implemented based on the operational grid-point statistical interpolation (GSI) system to assimilate high-resolution data in the TC inner core regions in 2017, for the operational HWRF system (Tallapragada et al., 2012 [15], Tallapragada et al., 2014 [16]). In the operational HWRF, inner core DA is run after vortex initialization (VI), which includes vortex relocation and initial vortex size and intensity adjustment (Liu et al., 2012 [17]). Therefore, in order for an ensemble-based DA system to provide both dynamically and observationally consistent initial model conditions in the TC inner core region, two components need to be considered: the accurate TC vortex initial position and the structures that are provided by assimilating inner-core observations, and a set of ensemble forecasts that provide a self-consistent and flow-dependence background error covariance matrix. Tong et al. (2018) [18] evaluated the impact of assimilating high-resolution inner-core reconnaissance observations on tropical cyclone initialization and prediction. They demonstrated the benefits of inner-core DA in the HWRF system. They also noted that the mismatch between the simulated cyclone structure by VI and observations could cause initial spin-down, especially for strong storms, even though the inner-core data assimilation results in analysis that better fit the observations. The results suggested the importance of both the model VI process and model physics. On the other hand, the inner core DA system relies on model ensembles to determine the weights given to model fields by the use of an ensemble background covariance matrix. Pu et al. (2016) [19] proved that the use of HWRF-generated ensembles can further improve the initial TC structure analysis, and hence improve the model track and intensity forecasts. In the 2017 version of the operational HWRF system, the self-cycled DA was first introduced for the priority Atlantic storms identified by the National Hurricane Center (NHC). One of the important upgrades in the 2018 HWRF system was to further add stochastic perturbations to model physics on top of initial field perturbations to account for the stochastic characteristics of model physics. The paper is organized as follows. Section 2 briefly describes the 2017 version of the operational HWRF system followed by the description of the 2017 version of the HWRF DA system. The experimental design, dataset, the DA upgrade and model physics perturbations are described in Section 3. The impact of the stochastic physics-based hybrid GSI/EnKF data assimilation on the TC track and intensity forecasts, using the operational HWRF system, are analyzed and evaluated in Section 4. Concluding remarks are provided in Section 5.

2. The HWRF Model and DA System

2.1. The HWRF Model

The HWRF is a high resolution, cloud-resolving operational TC forecast system at the NCEP, which was developed based on the Non-Hydrostatic Meso-Scale Model (NMM) dynamical core of the WRF model, with the model physics specifically tuned and the components specifically designed for TC prediction, including the ocean model, the wave model, the NCEP coupler, which couples the atmospheric model with ocean/wave models, the vortex initialization, data assimilation, and the vortex tracker. Currently, the Princeton Ocean Model (POM) is used for North Atlantic (NATL), Eastern Pacific (EPAC) and Central Pacific (CPAC) oceans. HWRF has undergone continuous improvements and upgrades every year since it became the operational TC forecast system at the NCEP in 2007 (Tallapragada et al., 2012 [15], Tallapragada et al., 2014 [16], Mehra et al., 2018 [20]). HWRF provides the TC track and intensity as well as the TC structure forecast guidance, which are widely used by the operational TC forecast centers in the world over all global oceanic basins. In this study, we used the 2017 version of the HWRF (referred to H217 hereafter) with two major upgrades. One is the horizontal resolution increase from 18/6/2 km to 13.5/4.5/1.5 km. The other is the introduction of stochastic model physics into the HWRF ensemble forecast to provide a background error covariance matrix for the inner-core DA. The vertical resolution is 75 hybrid pressure-sigma levels. A detailed HWRF system description can be found in Bernardet et al., 2011 [21]; Tallapragada et al., 2014 [16], Tallapragada et al., 2012 [15]; Mehra et al., 2018 [20].

2.2. The DA in H217

The HWRF DA system is one of the most important components in the operational HWRF. A special vortex initialization procedure was utilized to provide a first guess for data assimilation. The vortex initialization extracts the TC vortex from the HWRF 6 h forecast, and makes an adjustment to its size and intensity according to the storm message files or TC vitals, which is used by the NCEP and many other modeling groups to begin the process of running a hurricane model. The vortex in the 6 h forecast from the Global Data Assimilation System (GDAS) is then removed and replaced the adjusted vortex in the observed location. The technical details of the vortex initialization can be found in Liu et al. (2012) [17] and Tong et al. (2018) [18]. The TC inner-core DA assimilates the observations on top of the results of vortex initialization. Figure 1 describes the dual-resolution self-cycled DA system used in H217. The DA system used in the operational HWRF is a GSI-based hybrid EnKF Variation system, which uses either the 6 h forecast 80-member ensemble from the GDAS, or the 6 h forecast 40-member ensemble generated by the HWRF (self-cycled DA) to calculate the model background error covariance matrix. Due to the computer resource limit, the HWRF 40-member self-cycled DA system is only run for high priority TCs targeted at the National Oceanic and Atmospheric Administration (NOAA) P3 Tail Doppler Radar (TDR) missions. Details of the HWRF DA system can be found in Lu et al., 2017 [22] and Tong et al., 2018 [18]. An ideal TC inner-core DA system should include: 1. flow-dependent background error covariance matrix that is generated using the same dynamic model system; 2. the mass-wind balanced initial conditions to avoid model initial shock; 3. An ensemble that represents the realistic model uncertainties and forecast errors. The first two issues listed above have been addressed in H217 by implementing the dual-resolution self-cycled DA system and blending DA increments with vortex initialization in the inner-core area for storms with a maximum 10 m wind speed greater than 65 kt (Biswas et al., 2017 [23]). A realistic ensemble forecast error sampling will lead to a better estimate of model forecast background error covariances, hence a more realistic TC initial analysis. In H217, only uncertainties from model initial conditions, such as the TC initial intensity, initial position and initial large-scale flows, are accounted for through EnKF, while the uncertainties due to model physics are not considered.

3. Experiment Design

Two experiments were designed to study the impact of self-cycled inner-core DA, with model physics uncertainties added to the HWRF ensemble system, on the initial HWRF TC and large-scale environment flow analysis, and hence on the track and intensity forecasts. The first experiment served as a control experiment (referred to H40C hereafter), and used the same configuration used in H217 except that the experiment always uses the self-cycled DA system. Another experiment (referred to H40P hereafter) is the same as H40C, but adding stochastic physics model perturbations in the self-cycled DA system. The stochastic perturbation approach provides a methodology to represent the uncertainties due to model physics in ensemble forecasts, and also has the capability of reducing systematic error through the concept of nonlinear noise-induced rectification. In this study, the three important perturbations and uncertainties are added in the HWRF model physics, including stochastic convective trigger perturbations, stochastic perturbations on planetary boundary layer (PBL), and stochastic surface drag coefficient perturbations. The simplified Arakawa–Schubert (SAS) (Arakawa et al., 1974 [24]) convection scheme was used, in which the convective trigger function is perturbed by a uniformly distributed random number generator in the magnitude of extrema of ±50 hPa, while the PBL heights and surface drag coefficients are stochastically perturbed by a uniformly distributed random number with an extrema of ±20% of their values. All the stochastic perturbations are added at every model time step. It should be noted that these perturbations by no means include all the internal uncertainties related to the model physics. Our choice of generating the ensemble with the convective trigger function was based on our experiments and experience with different model physics changes, which showed that TC intensity forecasts are sensitive to the convective trigger function, PBL heights, and surface drag coefficient, as compared to other physical changes (Zhang et al., 2014 [8]).

The detailed descriptions of these perturbations can be found in Zhang et al., 2014 [8]. Both experiments are conducted for all TCs in the 2016 Atlantic hurricane season.

4. Results and Discussions

4.1. Ensemble Spread Comparison with and without Stochastic Model Physics

One of the challenging issues for hurricane EPS is that the ensemble spread of many predicted variables, especially wind fields (Zhang et al., 2014 [8]), is generally under-dispersed. In other words, the hurricane intensity forecast errors are under-represented. There will be notable consequences when the under-dispersed TC intensity ensemble is used in the Hybrid GSI/EnKF DA system, which assigns weights to both first guess fields and observations based on the background error covariances provided by the ensemble forecast and the observational error variances, respectively. The narrower ensemble spread causes the DA system to be overconfident of the first guess fields and implies less impact by the observations, and assigns less weights to the observations, known as filter divergence. Therefore, in order to evaluate if the ensemble system benefits from adding stochastic physics to the model, the ensemble spread of various model variables were compared between the control experiment (H40C) and the stochastically perturbed physics experiment (H40P) for all the hurricanes in 2016. In this study, the standard deviation was used to represent the ensemble spread. For brevity, we used Hurricane Hermine 2016 as an example, which was named a hurricane on August 31, 18Z and dissipated on September 06, 12Z. Figure 2 clearly demonstrates that adding stochastic physics to the HWRF ensemble system increases the ensemble spread of domain-averaged 10 m wind speed, 2 m temperature, and the mean sea level pressure fields. The ensemble spread increases can be as large as 0.5 m/s for the 10 m wind speed, 0.15 °C for T2 m, and 0.5 hPa for the mean sea level pressure. The vertical distribution of the spatially averaged ensemble wind speed field, and the spatially averaged ensemble temperature field spread are also examined in Figure 3, which further confirms that ensemble spread produced by H40P is indeed larger than that of H40C. The maximum spread difference between H40C and H40P occurs the around model sigma level 20, or about 850 hPa. Figure 4 shows an example of the comparisons between the ensemble spread of 10 m wind speed and 2 m temperature for Hurricane Matthew on 29 September 2016, 00:00 Universal Time Coordinated (UTC). Both 10 m wind speed and 2 m temperature fields show larger ensemble spread magnitudes and broader areas of ensemble spread, which indicates that the hurricane position spread was also increased by adding stochastic physics in H40P.

4.2. Impact of Stochastic Model Physics on Hurricane Analyses

Improved HWRF ensemble system in the hybrid GSI/EnKF data assimilation system will make the better use of the available observations and result in the model analysis that are close to the observations through the improved model background error covariance matrix. The values of the observation-minus-first-guess (OmF) and observation-minus-analysis (OmA) are generally used to assess the overall impact of observations on model initial conditions and analysis errors. Two 2016 hurricanes, Hermine and Matthew, were used to evaluate the impact of the proposed new data assimilation, i.e., the introduction of stochastic physics into the HWRF-based ensemble system. The reasons to choose these two storms are two-fold: both were long lasting storms that posed threats to the US coasts (Figure 5), and the storm inner-core TDR radial velocity observations were available when storms were approaching US coastal areas. Hurricane Hermine developed in the Florida Straits on 28 August 2016, from a long-tracked wave. It intensified into a 130 km/h category hurricane just before making landfall in the Florida Panhandle on September 2. Hurricane Matthew formed on 28 September 2016, and was the first Category 5 hurricane since 2007, which caused catastrophic damage. The TDR data from the NOAA P-3 aircraft has been assimilated in the operational HWRF system since the 2013 hurricane season, and has demonstrated a positive impact on the hurricane track and intensity forecasts (Tong et al., 2018 [18]). The detailed quality control processes of the TDR data can be found in Gamache 2005 [25].

Figure 6 shows comparisons of OmF and OmA between H40C and H40P for the cycles of hurricane Hermine in which the TDR radial velocity data is available and assimilated. The advantages of H40P over H40C are clearly manifested in that for the values of OmA from H40P are smaller or improved than that of H40C throughout most of the TDR cycles, except for the first cycle, while for the values of OmF are comparable between H40C and H40P. The cycle- and domain-averaged vertical OmA profile further confirmed that the improvement of the initial conditions was shown in vertical levels between 800 hPa and 300 hPa (Figure 6c). Similar results are shown in Figure 7 for Hurricane Matthew with larger improvements on the OmA vertical profile and the improvements now extending to the surface. It is worth noting that the analysis fields of other variables such as wind components u, v, and temperature are also improved in H40P compared to that in H40C (not shown).

4.3. Impact of Stochastic Model Physics on Model Forecasts

The impact of the proposed self-cycled HWRF stochastic physics DA system on the track and intensity forecast skills are examined for all 2016 hurricanes with a total verifiable sample of 291 at the initial time and 91 at day 5. The forecast track and intensity skills of H40P relative to H40C is defined as follows:

Skill Score = (E_H40C−E_H40P)/E_H40C

where E_H40C and E_H40P denote the forecast errors from the H40C and H40P experiments, respectively. A positive (negative) value represents a smaller (larger) forecast error and a better (worse) performance compared to HWRF. Figure 8 indicates that the addition of stochastic physics in the HWRF DA had a neutral to slightly positive impact on the track forecast skills, with one exception of degradation at 18 h. There was about an 8% improvement at all lead times on the intensity forecast skills. The forecast skill improvements in H40P in general merited the proposed scheme to be included in the 2018 version of the operational HWRF upgrades. In order to ensure the robustness of the track and intensity improvements in H40P, Figure 9 compares the frequency of the superior performance (FSP) between the H40C and H40P, which further confirms that H40P generally outperforms H40C in terms of track and intensity forecasts, especially in the later forecast hours. The improvements of the track and intensity forecasts of H40P for the individual storms/cycles are examined. Figure 10 is an example of the 5-day track and intensity forecast comparisons between H40C and H40P for Hurricane Matthew 14 L, 00:12 UTC 28 September 2016. It shows that the track forecasts of H40P are comparable to that of H40C, but the H40P track is smoother and realistic. The intensity forecasts of H40P shows a clear advantage over H40C, especially during the Rapid Intensification period forecast hours of 48–72 h. The stratified verification was also performed for the initially weak (<50 knots) and strong (≥50 knots) storms. The results show that compared to the control experiment (H40C), the stochastic physics DA experiment (H40P) had a more positive impact on the track forecast skills (~3%) and the intensity forecasts skills (average ~8%) for the weak storms (Figure 11c,d) than that for the strong storms (Figure 11a,b), whose impacts were basically from a neutral to slightly negative impact.

In the previous section, we have shown that when the stochastic physics are applied to the HWRF self-cycled DA system, the storm initial wind field errors represented by OmA were reduced for Hurricanes Hermine and Matthew, for which the TDR radial wind observations are available. Improved model initial conditions generally will lead to improved model forecasts. Verification for the two TDR storms, Hurricane Hermine and Matthew, is shown in Figure 12, which manifests that compared to H40C, the H40P track forecast skills are indeed improved at nearly all forecast lead times. The H40P intensity forecast skills show variable to slightly negative for earlier forecast hours, and are improved at the later forecast hours, ~12% at day 4. The negative intensity impacts at earlier forecast hours are mainly caused by model spin down issue, due to the inconsistence between the HWRF vortex initialization and the in-core DA. The issue is discussed in Tong et al., 2018 [18]. It is worth noting that compared to H40C, the H40P track and intensity forecast skills are more improved for Hurricane Hermine than that for Hurricane Matthew. For example, the averaged H40P forecast error for Hurricane Hermine is reduced ~60 km for track and ~3.1 m/s for intensity at day 4 compared to H40C, while for Hurricane Matthew, the error reduction at day 4 is ~2 km for track and 0.5 m/s for intensity (Figure is not shown here). This could be related to the TDR flight observation sampling pattern. In general, the stochastic physics DA has a positive impact on both model analyses and forecasts in terms of track and intensity forecast skills.

5. Concluding Remarks

In this study, a set of stochastic physics perturbations, including perturbations in the cumulus convection scheme, in the PBL scheme, and in surface layer scheme, was introduced into the self-cycled hybrid GSI/EnKF DA system in NOAA’s operational HWRF. The impacts of stochastic physics in DA on the hurricane track and intensity forecast skills were investigated in several aspects by conducting two experiments for all 2016 Atlantic hurricanes. The control experiment (H40C) was the H217 configuration with the self-cycled hybrid GSI/EnKF DA system that utilizes a 40-member HWRF-based ensemble with no stochastic model physics. The other experiment (H40P) used the same configuration as H40C but with stochastic physics perturbations added in the HWRF-based ensemble system. The ensemble spreads from the two experiments were first compared to ensure that H40P configuration resulted in the desired ensemble spreads and statistical characteristics. The results showed that H40P produced a larger ensemble spread of domain-averaged wind speed, temperature, and mean sea level pressure fields than those of H40C. The results were further confirmed by the domain- and cycle-averaged vertical profiles of wind speed. The OmF and OmA from the two experiments were then examined. The two hurricanes, Hermine and Matthew of 2016, which posed a threat to US coastal areas, were selected for case studies involving high-resolution, inner-core reconnaissance TDR velocity observations taken by NOAA’s P3 aircraft. The approach used in H40P most often produced values of both OmF and OmA that were smaller than H40C, indicating improved initial conditions for tropical storm forecasts. The OmF and OmA reduction of using stochastic physics in HWRF DA is clearly demonstrated for the two TDR storms.

Three sets of verification were performed to study the impact of stochastic physics in HWRF DA on the track and intensity forecasts. The forecast skills were first compared between H40C and H40P for all 2016 Atlantic hurricanes. The verification was then conducted for initially strong (≥50 kt) and weak (<50 kt) storms. Finally, we also compared the forecast skills between H40C and H40P for hurricanes for which TDR observations were assimilated. It was evident that the H40P configuration improved the track and intensity forecast skills which averaged ~3% and ~6% over H40C, respectively, for all 2016 hurricanes. It was also found that the additional stochastic physics perturbations introduced to the DA system have more positive impacts on the track and intensity forecasts for weaker storms than for stronger storms. This could be attributed to the initial unbalance between the model first guess field and the observation data, or the spin-down issues, especially for stronger storms, explained by Tong et al., 2018 [18], while H40P provides more realistic background error covariance matrix. The verification also showed that the track and intensity forecasts for the two hurricanes with TDR data assimilated and also benefited from introducing stochastic physics in the generating ensemble for the DA system.

Based on the analysis and evaluation shown in this study, the stochastic DA system was then included in the FY2018 HWRF upgrades.