A Prototype Regional GSI-based EnKF-Variational Hybrid Data Assimilation System for the Rapid Refresh Forecasting System: Dual-Resolution Implementation and Testing Results

Studies have shown that background error covariance (BEC) plays an important role in atmospheric data assimilation (DA), and especially for the mesoscale and convective scale where weather systems are more transient and intermittent (Zhang et al., 2006; Meng and Zhang, 2007; Liu and Xue, 2008; Ancell et al., 2014). In typical 3D variational (3DVar) systems, static BEC is usually used in the cost function, which does not contain flow-dependent spatial covariance (Lorenc, 1986); any cross-covariances usually only reflect static, large-scale balances (Barker et al., 2004; Descombes et al., 2015). For smaller scale flows, flow-dependent spatial covariance and cross-covariances tend to have increasing importance.

The ensemble Kalman filter (EnKF), originally proposed by (Evensen, 1994), uses flow-dependent BEC derived from an ensemble of forecasts. EnKF has become increasingly popular in recent years and has been applied to operational (e.g., Buehner et al., 2010a, b) or prototype numerical weather prediction (NWP) systems (e.g., Hamill et al., 2011; Wang et al., 2013). However, due to the relatively small ensemble size dictated by available computing resources, sampling error in the ensemble-derived BEC is usually significant (Hamill et al., 2001; Miyoshi et al., 2014; Anderson, 2016). Covariance localization is typically used to help alleviate the problem, which can affect flow balances (Hamill et al., 2001; Greybush et al., 2010). As an alternative way of alleviating the problem, (Hamill and Snyder, 2000) proposed a hybrid algorithm, in which a weighted average of the static and flow-dependent covariances is used within a 3DVar framework; they found that the hybrid system outperforms EnKF when the ensemble size is small or when the model error is large. An alternative computationally more efficient implementation of the hybrid idea, called the extended control variable (ECV) method, was later proposed by (Lorenc, 2003). The ECV method extends the original control vector in the 3DVar cost function by adding an additional term——the extended control variable term preconditioned upon the square root of the ensemble covariance. (Wang et al., 2007) showed that the method of (Hamill and Snyder, 2000) using a simple weighted average and the ECV method are mathematically equivalent. With the ECV approach, a hybrid 3D ensemble variational (3DEnVar) algorithm is relatively easy to implement within an existing 3DVar framework. Recent studies have demonstrated that forecasts initialized from analyses of a hybrid algorithm are usually better than those from traditional 3DVar, and are comparable to or better than those from pure EnKF (Wang et al., 2009; Li et al., 2012; Zhang and Zhang, 2012; Wang et al., 2013; Zhang et al., 2013; Pan et al., 2014; Schwartz and Liu, 2014). This approach has been used in several global NWP (Kuhl et al., 2013; Wang et al., 2013; Kleist and Ide, 2015) and regional mesoscale systems (e.g., Wang et al., 2008a, b; Schwartz and Liu, 2014; Schwartz et al., 2015; Wu et al., 2017).

Most of the hybrid DA studies mentioned above employed single-resolution (SR) configurations, in which the deterministic EnVar analysis and the ensemble analyses and forecasts are performed at the same resolution. For a high-resolution (HR) DA system, SR configuration can be computationally expensive, especially for operational purposes. For these reasons, dual-resolution (DR) hybrid schemes have been proposed, which combine HR background forecasts and hybrid analysis with BECs derived from low-resolution (LR) ensemble forecasts that are typically provided by a cycled EnKF system (Buehner et al., 2010b; Hamill et al., 2011; Clayton et al., 2013; Kleist and Ide, 2015).

For regional models, one DR hybrid example was recently reported in (Schwartz et al., 2015), based on Weather Research and Forecasting (WRF) DA systems. (Schwartz et al., 2015) compared the WRF SR and DR hybrid analyses and forecasts over a ∼3.5-week period. They found that the 45/15-km DR analyses were completed around three times faster and required about one quarter of the disk space of the 15-km SR analyses. Moreover, forecasts launched from the 15-km SR hybrid analyses had no significant differences from those launched from the 45/15-km DR hybrid analyses, although forecasts from the SR system captured more fine-scale features. The DR hybrid system consumes much less computational resource than the 15-km SR hybrid system.

The Rapid Refresh (RAP) forecasting system is an hourly-updated operational DA/prediction system of the U.S. National Weather Service using a 13-km horizontal grid spacing, and it replaced the operational Rapid Update Cycle (Benjamin et al., 2004) system as a regional operational analysis and forecast system in 2012. RAP uses the WRF-ARW model (Skamarock and Klemp, 2008) for forecasting and the Gridpoint Statistical Interpolation (GSI) analysis system (Wu et al., 2002; Kleist et al., 2009) for DA. GSI 3DVar was used until February 2014, when it was replaced by a hybrid 3DEnVar using flow-dependent covariances derived from the Global Forecast System (GFS) EnKF system. An updated version was implemented in August 2016 using 75% flow-dependent covariances derived from a GFS 80-member ensemble run at an ∼30-km grid spacing (Benjamin et al., 2016). In a sense, the operational RAP hybrid DA system is using a DR algorithm, although the ensemble forecasts are from a completely different model, and the GFS EnKF ensemble forecasts are only available four times a day; therefore, many hourly RAP hybrid analyses share forecasts from the same cycle. Such a configuration is clearly not optimal. Thus, a self-consistent EnKF DA system for RAP itself, running at LR, is desirable for maximum consistency and computational cost saving.

In fact, a regional GSI-based EnKF system using the EnSRF (ensemble square-root filter) (Whitaker and Hamill, 2002) had already been established for RAP in (Zhu et al., 2013), and tested at a grid spacing (40 km) three times that of the operational RAP. The same RAP operational data stream, with 3-h intervals over a 9-day period, was assimilated, and 18-h deterministic forecasts were evaluated. Results showed that the EnKF system was consistently better than the parallel GSI 3DVar system. Building on the EnKF system tested in (Zhu et al., 2013), (Pan et al., 2014) implemented a coupled EnKF-3DEnVar hybrid system and evaluated it at the same 40-km resolution as (Zhu et al., 2013) for the same test period. With equal weights given to the static and flow-dependent covariances, the hybrid system outperformed GSI 3DVar, and was either comparable to or slightly better than the EnKF system (Pan et al., 2014).

The EnKF and 3DEnVar hybrid systems implemented for RAP in (Zhu et al., 2013) and (Pan et al., 2014) were tested at a reduced resolution with a 40-km grid spacing because running the ensemble analyses and forecasts at the full ∼13-km resolution is expensive. With a DR hybrid system, a single hybrid analysis is performed at the native ∼13-km grid spacing to provide initial conditions for the RAP deterministic forecast (the operational RAP is deterministic at this time), using ensemble forecasts produced at the 40-km grid spacing at a much lower cost. Testing and evaluating an ∼13-km hybrid system coupled with EnKF cycles at the 40-km grid spacing is the main goal of this paper.

The rest of the paper is organized as follows: section 2 describes the implementation of the coupled DR EnKF-3DEnVar hybrid system for RAP; the design of the testing experiments is given in section 3; the results and comparisons of the experiments are discussed in section 4; section 5 provides conclusions and additional discussion.

The GSI-based DR EnKF-3DEnVar hybrid system for RAP is based on the EnKF system reported in (Zhu et al., 2013) and the 3DEnVar hybrid system reported in (Pan et al., 2014); both were tested at a reduced SR with a 40-km grid space. Two-way and one-way coupling can be applied in a DR hybrid system (Clayton et al., 2013; Schwartz et al., 2015). For one-way coupling, EnKF only provides the LR ensemble covariance for HR hybrid 3DEnVar analyses. In the two-way DR hybrid, an additional step is performed that re-centers the LR ensemble about the HR analysis. However, (Clayton et al., 2013), (Wang et al., 2013) and (Schwartz et al., 2015) found relatively small impacts of performing such a re-centering. Furthermore, keeping the DR system one-way coupled with the LR EnKF system makes the intercomparisons cleaner. Thus, one-way coupling is employed in this study.

Figure 1.  Flowchart of a full GSI-based EnKF-3DEnVar DR hybrid DA cycle with one-way coupling between the EnKF (upper portion at 40-km horizontal grid spacing) and 3DEnVar DR hybrid analysis (lower portion at 13-km horizontal grid spacing). The arrows pointing from EnKF downward to 3DEnVar indicates the prevision of 40-km ensemble perturbations from EnKF to hybrid analyses. GSI is used as the DA system for 3DEnVar.

Figure 1 shows the flowchart for the one-way-coupled DR EnKF-3DEnVar hybrid cycles used in this paper. As with the SR hybrid system documented in (Pan et al., 2014), the GSI system is used for observation processing, including data quality control, thinning and calculation of the innovations. The EnKF system directly ingests observation innovations processed by GSI and produces an ensemble of analyses. The EnKF analyses and forecasts are run on the 40-km grid. In the DR system, a single deterministic forecast and 3DEnVar hybrid analysis are produced on the HR grid of ∼13-km grid spacing in each cycle, using the 40-km ensemble forecasts from the EnKF cycles for flow-dependent BECs. In the rest of this paper, LR refers to the 40-km grid spacing, and HR refers to the ∼13-km grid spacing.

Specifically, the DR hybrid algorithm is described below, with the presentation of the general algorithm mostly following (Pan et al., 2014) and (Schwartz et al., 2015). The analysis increment vector, δx, of the DR hybrid can be expressed as \begin{equation} \label{eq1} {\delta}{x}={\delta}{x}_{1}+{SD}{\alpha} , \ \ (1)\end{equation} where δx₁ is the HR analysis increment vector associated with the static background covariance with a length of n _h, and SDα is the increment associated with flow-dependent covariances from the K-member ensemble. Equation (2) includes an interpolation operator S, that interpolates the Dα from LR to HR space, while it is an identity matrix for an SR 3DEnVar algorithm. D is a n_l×(Kn_l) matrix defined as \(D=[\rm diag(x'_l1) \rm diag(x'_l2) \ldots \rm diag(x'_lK)]\), where n_l is the dimension of the state vector on the LR grid, and l indicates the lower resolution. diag() is an operator that converts vector x'_li into a diagonal matrix with size n_l× n_l (Wang, 2010; Schwartz et al., 2015). The subscript i indicate the i th member of the ensemble. α is an extended control variable with column vector length of Kn_l formed by concatenating extended control variables for each ensemble member α_K as \begin{equation} \label{eq2} {\alpha}=\left[ \begin{array}{l} {\alpha}_{1}\\ {\alpha}_{2}\\ \vdots\\ {\alpha}_{K} \end{array} \right] . \ \ (2)\end{equation} The corresponding cost function to obtain the analysis increment δx₁ and α is expressed as \begin{eqnarray} \label{eq3} J({\delta}{x}_{1},{\alpha})&=&\beta_1J_{\rm b}+\beta_2J_{\rm e}+J_{\rm o}\nonumber\\ &=&\dfrac{1}{2}\beta_1{\delta}{x}_{1}^{\rm T}{B}^{-1}{\delta}{x}_{1}+\dfrac{1}{2}\beta_2{\alpha}^{\rm T}{A}^{-1}{\alpha}\nonumber\\ &&+\dfrac{1}{2}[{H}{\delta}{x}-{y}']^{\rm T}{R}^{-1}[{H}{\delta}{x}-{y}'] , \ \ (3)\end{eqnarray} where J _o is observation term, H is the tangent-linear observation operator, y' is the observation innovation vector, B is the static BEC, and R is the observation error covariance. A is a Kn_l× Kn_l block diagonal matrix that prescribes the covariance localization scale for the flow-dependent covariance derived from the ensemble forecasts. β₁ and β₂ in front of J _b and J _e are the weights given to static and ensemble BECs, respectively, and they are constrained by \begin{equation} \label{eq4} \frac{1}{\beta_1}+\frac{1}{\beta_2}=1 . \ \ (4)\end{equation} The gradients of the cost function with respect to δx₁ and α have the form \begin{eqnarray} \label{eq5} \nabla_{{\delta x}_{1}}J&=&\beta_1{B}^{-1}{\delta x}_{1}+{H}^{\rm T}{R}^{-1}({H\delta x}-{y}') ,\ \ (5)\\ \label{eq6} \nabla_{{\alpha}}J&=&\beta_2{A}^{-1}{\alpha}+{D}^{\rm T}{S}^{\rm T}{H}^{\rm T}{R}^{-1}({H\delta x}-{y}') , \ \ (6)\end{eqnarray} where D ^T and S ^T are the adjoints of D and S in Eq. (2), respectively. S ^T is applied to H ^TR^-1(Hδ x-y'), while S is applied to Dα from the LR space to HR space.

In the RAP hybrid DA system, WRF is used as the forecast model. The Model Evaluation Tools (Brown et al., 2009) developed by the Developmental Testbed Center is used for forecast verification.

Figure 2.  Example of the horizontal distributions of (a) sounding, profile and VAD (Velocity Azimuth Display), (b) surface stations over land and for ships, (c) GPS-PW (Precipitable Water) and GPS-RO (Radio Occultation), and (d) aircraft observations at 0000 UTC 8 May.

As stated earlier, the DR hybrid system uses a 40/∼13-km horizontal grid spacing combination. The LR domain at 40-km grid spacing for the ensemble covers North America with 207× 207 horizontal grid points (bold box in Fig. 2a), while the HR domain at ∼13-km grid spacing has 616× 616 horizontal grid points covering roughly the same domain (the HR domain is not plotted in Fig. 2). Both domains have 50 vertical levels extending up to 10 hPa at the model top using terrain-following hydrostatic-pressure-based vertical coordinates that stretch with height (Skamarock et al., 2008). The static background error statistics calculated based on the NCEP North American Model forecasts using the National Meteorological Center method, as provided in the GSI system (Hu et al., 2016), are used in this study. The error statistics are latitude and sigma-level dependent only; they are interpolated to the analysis grid within GSI. The flow-dependent BECs are derived from the ensemble forecasts provided by the EnKF system at the 40-km grid spacing.

As in (Pan et al., 2014), the operational data stream of RAP excluding satellite radiance data is assimilated in the DR hybrid system. The distributions of the data at 0000 UTC May 8 are shown in Fig. 2. Eighteen-hour forecasts from the analyses are verified against surface and sounding data in the HR domain. The surface data verified include surface pressure, 2-m relatively humidity (RH), 2-m temperature (T), 10-m zonal and meridional wind components (U and V, respectively); the sounding observations include RH, T, U and V.

The root-mean square error (RMSE) is used to evaluate the forecasts, and the bootstrap resampling method (Candille et al., 2007; Buehner and Mahidjiba, 2010; Schwartz and Liu, 2014) following (Pan et al., 2014) is used to determine the statistical significance of error differences. The RMSEs are calculated against the observations at certain levels and forecast hours first, and then aggregated over all cycles for specific forecast hours for skill evaluation.

To assess the statistical significance, bootstrap resampling is performed. New samples are created by randomly drawing from the dataset 3000 times, allowing the same data to be drawn more than once. With the resample, we calculate the aggregated RMSEs along with a two-tailed confidence interval from 5% to 95%. As in (Pan et al., 2014), RMSE differences are calculated between a specific experiment and its benchmark first. The bootstrap method is then applied to the RMSE differences with confidence intervals from 5% to 95% to determine the significance of improvement. When all confidence intervals of the RMSE differences are below/above zero, the experiment is significantly better/worse than the benchmark experiment at a 90% confidence level. Additional discussion on the use of the bootstrap method for calculating the statistical significance of forecast differences can be found in (Pan et al., 2014), (Schwartz and Liu, 2014), and (Xue et al., 2013).

The Gilbert skill score (GSS) (Gandin and Murphy, 1992) is used to evaluate the precipitation forecast skills of 12-h deterministic forecasts against the 4-km NCEP Stage IV precipitation data in the CONUS (Conterminous United States) domain in which the data are available (Lin and Mitchell, 2005).

As in (Zhu et al., 2013) and (Pan et al., 2014), the same test period from 8-16 May 2010 is used, which contained active episodes of convection. All DA experiments start at 0000 UTC 8 May and end at 2100 UTC 16 May with continuous three-hourly cycles. The initial fields and boundary conditions are interpolated from operational GFS analyses and forecasts. Random perturbations are created by the random CV3 option in the WRF DA system (Barker, 2005; Barker et al., 2012) and added to the GFS analysis initial condition at 0000 UTC 8 May 2010 to start the ensemble forecasts for the EnKF and the GFS forecasts to create perturbed ensemble boundary conditions (Torn et al., 2006).

(Pan et al., 2014) suggested that the BEC weighting factor, as one of the important tuning parameters in a hybrid algorithm, has a great impact on the performance of the hybrid DA system. To examine the performance of the DR hybrid system and its sensitivity to the weighting factor, three DR hybrid experiments——namely, HyDR05, HyDR09 and HyDR10——are run using 50%, 90% and 100% weights given to the ensemble BEC, respectively. Because HyDR09 produces the best analyses and forecasts among all DR hybrid experiments, it is also called the DR hybrid control experiment, named HyDR_Ctl (Table 1).

The well-tuned 3DEnVar experiment (Hybrid1W_Ctl) at the 40-km grid from (Pan et al., 2014) is adopted as a benchmark, and the same well-tuned EnKF experiment with 40 members from (Pan et al., 2014) is also used in this study to provide the LR ensemble perturbations for the DR hybrid experiments. Hybrid1W_Ctl from (Pan et al., 2014) employed half and half static/flow-dependent covariance, a horizontal covariance localization scale of 300 km, and a vertical covariance localization scale of 0.3 in terms of natural logarithm of pressure. Experiment HyLR_HRF (Table 1) involves forecasts initialized from interpolated fields from the analyses of Hybrid1W_Ctl every cycle. A HR GSI 3DVar DA experiment, named VarHR_HRF, is also run at the ∼13-km resolution (Table 1). In other words, HR forecasts from the LR hybrid control experiments and the HR 3DVar analyses are used as references for evaluating forecasts from DR DA experiments.

The analyses of DR DA may benefit from the HR background forecasts, which contain more detailed flow structures. The configurations of HyDR05 are the same as HyLR_HRF, except that the deterministic background forecasts of HyDR05 are performed on the ∼13-km grid instead of the 40-km grid. In HyLR_HRF, forecasts are run at ∼13-km grid spacing from interpolated 40-km analyses of Hybrid1W_Ctl. The comparison between HyDR05 and HyLR_HRF isolates the impact of the increased background forecast resolution.

For variational minimization in either 3DVar or 3DEnVar, two outer-loop iterations and 100 inner-loop iterations are used. Evaluations are mainly based on forecasts on the ∼13-km grid, launched from either the HR analyses or fields interpolated from LR 40-km analyses. All experiments are listed in Table 1.

In (Pan et al., 2014), the skills of EnKF, 3DEnVar hybrid and 3DVar at a single reduced 40-km resolution were compared. The results showed 3DEnVar hybrid using 50% ensemble covariance significantly outperformed GSI 3DVar for all variables through the entire 18-h forecast period, while 3DEnVar hybrid had a comparable performance to EnKF overall. In this section, the performance of the DR hybrid system and its sensitivity to the weighting factor of BECs and the resolution of background forecasts are examined. In section 4.4, precipitation forecast skills are evaluated.

In (Pan et al., 2014), the lowest RMSEs were obtained when using 50% ensemble BEC in their 40-km control hybrid experiment, Hybrid1W_Ctl. At an ∼13-km grid spacing, smaller-scale features can be captured, which tend to be more transient and hence more flow-dependent. The analysis may benefit from a higher weight for the ensemble covariance. Experiments HyDR05, HyDR09 (also named HyDR_Ctl) and HyDR10 are compared to examine the impact of flow-dependent covariance in the DR hybrid system.

Figure 3.  Aggregated 3-h forecast RMSEs along with confidence error bars at different height levels verified against sounding data for (a) RH, (b) T, (c) U, and (d) V for experiments HyDR05, HyDR09, and HyDR10. The error bars represent the two-tailed 90% confidence interval (5% on the left and 95% on the right) using the bootstrap distribution method.

The aggregated 3-h forecast RMSEs verified against sounding data are shown in Fig. 3. The RMSEs at each pressure level were obtained by averaging values within a layer 50 hPa above and below that pressure from all cycles, except for the topmost and lowest levels. The 3-h forecasts are also used as the background in each DA cycle, and their errors can be used as a proxy for measuring the DA quality. The results show that HyDR09 has the smallest RMSEs for RH, U and V at almost all levels. For T, the RMSEs from HyDR09 are higher than those from HyDR05 above 800 hPa. The performance of HyDR10 is comparable to or worse than HyDR05 for RH below 600 hPa, and for T, U and V at all levels. These results indicate that, with the DR hybrid 3DEnVar system, when the grid spacing of the hybrid analysis as well as the background forecast is decreased from the 40-km used in (Pan et al., 2014) to ∼13-km, optimum results are obtained when the weight for the ensemble BECs is 90% (among the weights examined), instead of the 50% for the SR LR case. This may be because of the increased level of flow dependency of the background errors at HR, as suggested earlier. Raising the weighting factor for the flow-dependent covariances means that more mesoscale information can be involved in the DA. However, the forecasting skill of T to the weighting factor is opposite to RH, U and V at the middle to upper levels.

In this section, we compare the performance of experiments HyDR_Ctl (i.e., HyDR09) using hybrid 3DEnVar with VarHR_HRF, which uses the pure 3DVar DA method run at HR (see Table 1).

The 9-day aggregated RMSEs of the 3-h forecasts verified against sounding data at all levels are shown in Fig. 4. As shown in Fig. 4, VarHR_HRF underperforms HyDR_Ctl, with its errors being significantly larger for RH, U and V at most levels, while the errors for T are comparable.

The overall domain- and level-aggregated RMSEs verified against sounding and surface data are shown in Fig. 5 and Fig. 6, respectively, for analyses (hour 0) and forecasts at 3-h intervals up to 18 hours. HyDR_Ctl significantly outperforms VarHR_HRF at the analysis and forecast for all variables throughout the entire forecast period. The RMSEs of all variables are noticeably lower in the analyses than in the forecasts, and forecast errors increase quickly in the first three hours before becoming more stable thereafter; such rapid error growth is likely associated with fast small-scale error growth.

Overall, the DR coupled EnKF-3DEnVar hybrid scheme significantly outperforms the 3DVar scheme for all variables at all forecast hours when verified against soundings and surface observations. The results suggest the efficacy of using a DR configuration for a hybrid DA system.

Figure 4.  As in Fig. 3 but for experiments HyLR_HRF, HyDR_Ctl, VarHR_HRF and HyDR05. The error bars indicate the two-tailed 90% confidence interval using the bootstrap method with 5% on the left and 95% on the right.

The impacts of the HR background forecast are investigated by comparing HyLR_HRF with HyDR05, in which the only differences lie with the resolution of the background forecasts. HyDR_Ctl is also included in this section to assess the impacts of HR background forecasts and flow-dependent covariance.

The 3-h forecast RMSEs verified against sounding data (Fig. 4) show that HyDR05 underperforms HyLR_HRF for RH and performs comparably for T, U and V at most levels. When using a higher weighting factor of 90% for flow-dependent covariances in HyDR_Ctl, the RMSEs are smaller than those from HyDR05 for RH at all levels, except for T at 1000-800 hPa, and U and V at 500-300 hPa. These results suggest that HyDR_Ctl benefits from the HR with 90% flow-dependent covariances.

The comparisons among HyLR_HRF, HyDR_Ctl and HyDR05 of the domain-aggregated RMSEs for the analyses and forecasts up to 18 hours against sounding data are shown in Fig. 5. The RMSEs of HyDR05 are comparable or slightly worse than those from HyLR_HRF, while the RMSEs of HyDR_Ctl are significantly smaller than those of HyLR_HRF for all variables except T. At the lower resolution of 40 km, the best analyses and forecasts were obtained in (Pan et al., 2014) when equal weights were given to the static and flow-dependent covariances in the hybrid DA. As the grid resolution increases, the smooth static covariance becomes less appropriate, which explains why a higher ensemble covariance weight of 90% used in HyDR_Ctl is beneficial.

Figure 5.  The bar chart in each frame shows the RMSEs of forecasts verified against sounding data, aggregated over the entire domain and over the nine-day period. The lower panel shows the 90% confidence interval of the RMSE differences between HyDR_Ctl and VarHR_HRF or HyDR05 and HyLR_HRF for (a) RH, (b) T, (c) U, and (d) V, for different forecast hours. If the interval does not include zero, the difference is statistically significant at the 90% confidence level. The error bars in the histograms represent the two-tailed 90% confidence interval with 5% at the bottom and 95% on the top using the bootstrap distribution method.

Figure 6.  The bar chart in the upper panel of each frame shows the RMSEs of forecasts verified against surface station observations, aggregated over the entire domain and over the nine-day period, for (a) surface pressure, (b) 2-m RH, (c) 2-m T, (d) 10-m U, and (e) 10-m V for different forecast hours. Confidence error bars represent the two-tailed 90% confidence interval (5% at the bottom and 95% on the top) using the bootstrap distribution method. The lower panel of each frame shows the 90% confidence interval of the RMSE differences between HyDR_Ctl, VarHR_HRF or HyDR05 and HyLR_HRF.

The impacts of the HR background forecasts are further examined by verifying analyses and forecasts up to 18 hours against surface data (Fig. 6). HyDR_Ctl has significantly smaller RMSEs than HyLR_HRF for 2-m RH and 10-m U and V from the analysis time and throughout the entire forecast period, and for surface pressure except at a few forecast hours. Large differences are found in the RH errors between the HR 3DVar/DR hybrid and the LR 3DEnVar hybrid (Fig. 6b), suggesting that for the surface moisture field DA can benefit significantly from the increased background resolution, given the better resolution of terrain and mesoscale boundary layer structures. For 2-m T, smaller errors at the analysis time in HyDR_Ctl and VarHR_HRF than those in HyLR_HRF indicate a better fit of the analyses to surface T observations. However, the forecast errors of T in HyDR_Ctl and VarHR_HRF become larger after three hours of forecasting than those in HyLR_HRF.

The results seem to suggest that the humidity and wind fields benefit more from the higher background resolution with increasing flow-dependent covariance, while this is not necessarily the case for the temperature forecasts, at least when verified against conventional data in terms of the RMSEs. Experiments with various combinations of resolutions used in the analysis and forecasting steps shed some light on such complex behaviors (not shown), but are not enough to fully answer the questions. Results also imply a need for multi-scale DA algorithms that explicitly treat observations and background errors representing different scales (Li et al., 2015) and use scale-dependent (Buehner and Shlyaeva, 2015) and/or multi-scale covariance localization (Miyoshi and Kondo, 2013).

In this section, the precipitation forecasts from HyLR_HRF, HyDR_Ctl, VarHR_HRF and HyDR05 are verified against the 4-km NCEP Stage IV precipitation data. The GSS (Gandin and Murphy, 1992), also known as the equitable threat score, is calculated, as in (Pan et al., 2014), for the 0.1, 1.25 and 2.5 mm h^-1 thresholds.

The GSSs are shown in Fig. 7. That HyDR_Ctl outperforms VarHR_HRF for all thresholds and all forecast hours suggests that the analysis method is important for precipitation forecasting skill. The results are consistent with those of (Schwartz, 2016), who examined DR hybrid DA with a 20/4-km grid combination.

With the HR forecasts, HyDR05 has better skill than HyLR_HRF after five hours at the threshold of 0.1 mm h^-1, and at one to eight hours at the threshold of 2.5 mm h^-1. With more flow-dependent covariance being used, HyDR_Ctl shows the best skill among all experiments after 10 hours at the 0.1 mm h^-1 threshold, and generally all hours at the 1.25 and 2.5 mm h^-1 thresholds.

The results indicate that, for precipitation, especially heavier precipitation, there is a clear benefit to running the hybrid DA at HR (relative to the LR hybrid), and to using ensemble-derived flow-dependent covariance (relative to 3DVar). The improved precipitation forecasts are consistent with reduced errors in the analyses and forecasts of humidity.

Figure 7.  Aggregated precipitation GSSs of 13-km forecasts as a function of forecast length for thresholds of (a) 0.1 mm h^-1, (b) 1.25 mm h^-1 and (c) 2.5 mm h^-1.

Based on the NCEP operational GSI variational framework, a coupled EnKF-3DEnVar DR hybrid DA system, using a 40/∼13-km combination, is tested for the RAP model configuration. The LR ensemble is provided by a single LR EnKF system developed and tested in (Zhu et al., 2013), and the single LR EnKF-3DEnVar hybrid system tested in (Pan et al., 2014) is used as a benchmark to evaluate the performance of the coupled DR hybrid system. As in (Zhu et al., 2013) and (Pan et al., 2014), a nine-day period from 8-17 May 2010, containing active convective systems, is used to examine the performance of the DR hybrid DA system. The conventional data stream used by the operational RAP system is assimilated. The analyses and forecasts are verified against surface and sounding observations, while HR precipitation forecasts are verified against Stage IV precipitation data.

A 90% weight for the ensemble covariance produces the best forecasts with the DR hybrid system (HyDR_Ctl), while a 50% percent weight given to the flow-dependent ensemble covariance is found to have the best performance with the 40-km LR hybrid system in (Pan et al., 2014). The comparisons among HyLR_HRF, HyDR_Ctl and HyDR05 suggest that the hybrid analyses and forecasts can benefit from the HR background (3-h deterministic HR forecast) when the weight given to the ensemble covariance is larger. In comparison, the impacts of the HR background forecasts are limited when 50% static covariance is used. By increasing the weight for the ensemble covariance to 90% within the DR hybrid system, the humidity and wind fields are improved throughout the 18 hours of the forecast, and the improvements to the humidity fields are the largest. However, the response of temperature forecasting skill to the weighting factor is opposite to other variables at the middle to upper levels. The exact reasons require further investigation.

The overall benefits of performing the hybrid analyses at HR while still keeping the EnKF cycles at LR to reduce the computational cost are clear, especially for the precipitation forecasting skill. Such benefits can be larger when the grid resolution becomes convection-allowing or convection-resolving. Corresponding to the largest improvement of RH, the precipitation forecast skill from the DR hybrid system is higher than that from the HR 3DVar method, suggesting that the analysis method is as important as the analysis resolution for convection-allowing predictions. These results are consistent with (Schwartz, 2016), in which 4-km convection-allowing forecasts using 20/4-km DR hybrid DA based on the GSI system were examined. Their precipitation forecasts over the first 12 hours from 4-km 3DVar and hybrid 3DEnVar were better than the forecasts from corresponding downscaled 20-km analyses. All precipitation forecasts from their 4-km hybrid analyses were more skillful than those from their 4-km 3DVar analyses.