Abstract

At the National Centers for Environmental Prediction (NCEP), a new three-dimensional variational data assimilation (3DVAR) analysis system was implemented into the operational Global Data Assimilation System (GDAS) on 1 May 2007. The new analysis system, the Gridpoint Statistical Interpolation (GSI), replaced the Spectral Statistical Interpolation (SSI) 3DVAR system, which had been operational since 1991. The GSI was developed at the Environmental Modeling Center at NCEP as part of an effort to create a more unified, robust, and efficient analysis scheme. The key aspect of the GSI is that it formulates the analysis in model grid space, which allows for more flexibility in the application of the background error covariances and makes it straightforward for a single analysis system to be used across a broad range of applications, including both global and regional modeling systems and domains.

Due to the constraints of working with an operational system, the final GDAS package included many changes other than just a simple replacing of the SSI with the new GSI. The new GDAS package contained an upgrade to the Global Forecast System model, including a new vertical coordinate, as well as new features in the GSI that were never developed for the SSI. Some of these new features included changes to the observation selection, quality control, minimization algorithm, dynamic balance constraint, and assimilation of new observation types. The evaluation of the new system relative to the SSI-based system was performed for nearly an entire year of analyses and forecasts. The objective and subjective evaluations showed that the new package exhibited superior forecast performance relative to the old SSI-based system. The new system has been shown to improve forecast skill in the tropics and substantially reduce the short-term forecast error in the extratropics. This implementation has laid the groundwork for future scientific advancements in data assimilation at NCEP.

1. Introduction

Data assimilation is a key component in producing quality forecasts in an operational numerical weather prediction (NWP) system. Many operational centers currently utilize variational schemes (e.g., Derber et al. 1991; Lorenc et al. 2000; Rabier et al. 2000; Rawlins et al. 2007). Variational methods are used to create a best estimate of the state of the atmosphere at a given time using information from observations as well as a background, typically a previous model forecast. This is achieved through the minimization of an objective function that measures the weighted distance of the analysis from observations and the background. The respective weights given to each term are functions of what is known about the error characteristics of the observations and background (forecast). For example, it would not make sense to generate an analysis that draws very closely to an observation that typically has a large error associated with it and should not be trusted. Likewise, it makes sense to draw closely to observations in regions where the forecast model generates large forecast error.

The background term in the objective function allows for the formulation of the analysis problem in different ways, such as in gridpoint space (Wu et al. 2002, hereafter referred to as WPP02), observation space (Daley and Barker 2001), and spectral space (Parrish and Derber 1992; Rabier et al. 2000). Formulating the analysis problem in spectral space has the advantage that reasonable homogeneous, isotropic background error covariances are easily defined and applied (Fisher and Courtier 1995; Derber and Bouttier 1999). Additionally, many global forecast models are formulated in spectral space allowing for the analysis variables to be closely related to the actual model variables. However, in constructing the analysis problem in such a way, it becomes difficult (though not impossible, perhaps through the application of wavelets or other methods) to construct inhomogeneous and anisotropic error covariances (WPP02). This is not necessarily the case if the analysis problem is performed on a grid, where anisotropic, inhomogenous locally defined background error application can be achieved through the use of recursive filters (Purser et al. 2003a,b). A thorough review of various methods of modeling and applying background error covariances for variational data assimilation, including a summary of the advantages and disadvantages in the choice of formulation, can be found in Bannister (2008b).

The use of more robust background error covariances allows the analysis to spread information along such things as terrain- or synoptic-scale weather features like cold fronts. It is also easier to develop a comprehensive system that can be used for both regional and global applications if it is formulated in gridpoint space. This provided the motivation at the National Centers for Environmental Prediction (NCEP) to combine the advantages of its regional analysis scheme [Eta Model three-dimensional variational data assimilation (3DVAR)] with the global Spectral Statistical Interpolation (SSI) system, and develop a 3DVAR analysis system formulated in gridpoint space (WPP02). Although the current operational configuration and application of the Gridpoint Statistical Interpolation (GSI) algorithm do not apply fully anisotropic, inhomogeneous background error statistics, the implementation of a system that utilizes recursive filters will make the testing and transition of such an application much easier in the future.

The GSI analysis system is the realization of a streamlined effort among many collaborators to create a unified, efficient, modular 3DVAR scheme formulated in grid space. It is designed to run on various computing platforms, create analyses for different numerical forecast models, and remain flexible enough to be able to handle future scientific developments such as new observation types, improved data selection, new analysis variables, anisotropic background error covariances, and expansion to four-dimensional variational data assimilation. The GSI has already demonstrated its flexibility through implementations into the NCEP Global Data Assimilation System (GDAS), the regional assimilation system (Wu 2005), and real-time mesoscale analysis (De Pondeca et al. 2007), with a future implementation planned as part of the Rapid Refresh system (Devenyi et al. 2005, 2007).

Although the GSI has already been implemented into many different systems, this paper will focus on details related specifically to the implementation of the GSI into the GDAS that occurred on 1 May 2007. The next section will give an overview of the GSI system and its formulation, highlight differences from the SSI, and provide details specific to the GDAS. Section 3 will show a comparison between the GSI and SSI for single-observation and individual analysis tests. This section is followed by a summary of the impacts of including the GSI as part of an upgrade to the GDAS on the forecast skill of the NCEP Global Forecast System (GFS), derived from more than a year of retrospective testing. A discussion of recent developments and future plans will then be presented.

2. Global GSI configuration

As was the case for the SSI (Derber et al. 1991; Parrish and Derber 1992), the GSI is a 3D variational analysis system. In the GDAS, a new estimate of the atmospheric state (analysis) is required every 6 h to initialize a new 9-h global model forecast. Although the background used for each analysis is the previous 6-h forecast, a 9-h forecast is necessary to allow for time interpolation of asynoptic observations that fall within the 6-h analysis time window (i.e., time interpolation of the background is done between the 3-, 6-, and 9-h forecasts that covers the 6-h data window centered on the analysis time). The analyses are then used as the initial conditions for subsequent forecasts and the cycle continues. There is also an analysis created with an earlier data cutoff time (jobs are submitted 2 h and 45 min after each synoptic time), which is referred to as the GFS analysis, which is then used to run out the global model forecast to 15 days. Each GFS analysis uses the 9-h forecast from the GDAS for its background, as the GDAS cycle has the advantage of the later data cutoff (GDAS jobs are run a full 6 h after synoptic time). This allows for the assimilation of more observations (such as those with delayed arrival times) to create a better analysis and, presumably, a better 9-h forecast.

It is important to have initial conditions with an appropriate mass-wind balance relationship, so as to minimize the amount of spinup or spindown time in the model forecast. To maximize efficiency, streamfunction was chosen as the analysis control variable (WPP02). Although the definition of balanced temperature, divergent wind, and surface pressure increments from the streamfunction increments does an adequate job of defining a quasi-balanced increment, testing showed that this in itself was not enough to get the GSI to perform as well as the SSI analyses. Alternatives exist regarding the choice of analysis variable for the wind but were not pursued because of the difficulty in applying multivariate covariances (wind components), efficiency (potential vorticity), or error characteristics that differed from the other analysis variables (e.g., vorticity, which has much shorter length scales then do temperature and pressure).

A balance equation can be applied easily in spectral space, such as with the SSI and the linear balance equation (Parrish and Derber 1992), as the inverse Laplacian can be inexpensively solved for. Therefore, SSI analysis increments, through the use of a linear balance constraint, are better balanced when compared to increments resulting from the statistically derived relationship used in the GSI. It is difficult to enforce fully nonlinear dynamic relationships because of computational aspects as well as the linearity that is part of the minimization algorithm. To attempt to compensate for the deficiency in the balance and to reduce the noise in the analysis increment, various attempts were made to add a dynamical constraint penalty term to the GSI. This eventually led to the development of a tangent linear normal mode constraint on the analysis increment (Kleist et al. 2009). The addition of this balance operator to the GSI led to a significant improvement in the quality of the analyses and subsequent forecasts, and was a necessary addition prior to implementation and is utilized in all of the GSI experiments described in this manuscript.

With the exception of a reformulation of the moisture variable, all of the remaining analysis variables (streamfunction, unbalanced velocity potential, temperature, and surface pressure, as well as ozone and cloud water mixing ratio) and the corresponding multivariate regression matrices are the same as those outlined in WPP02. There are two choices for the analyzed moisture variable in the GSI, including a pseudo–relative humidity (Dee and da Silva 2003) and normalized relative humidity (Holm 2003). For the purposes of the global application of GSI, the normalized relative humidity option has been chosen, which subsequently allows for a multivariate coupling of the moisture, temperature, and pressure increments as well as flow dependence. For example, the new analysis variable allows for a moisture increment from a temperature observation, in addition to flow dependence, as the background error for the humidity variable becomes a function of the humidity field itself (Fig. 1).

Fig. 1.

Background RH [shaded, contour interval (CI) is 10%] and specific humidity analysis increment at 850 hPa (white contours, CI is 0.001 g kg−1) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for a GSI analysis valid at 1200 UTC 20 Apr 2007.

Fig. 1.

Background RH [shaded, contour interval (CI) is 10%] and specific humidity analysis increment at 850 hPa (white contours, CI is 0.001 g kg−1) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for a GSI analysis valid at 1200 UTC 20 Apr 2007.

The latitude-dependent background error statistics for the GSI implementation were estimated using the National Meteorological Center (NMC, now known as NCEP) method (Parrish and Derber 1992; Rabier et al. 1998). Though other methods exist for deriving estimates for the background error covariances (Bannister 2008a), the NMC method was chosen for its computationally efficiency as well as familiarity with the method. Due to the fact that this implementation involved model changes including a new vertical coordinate, the derivation of the background error statistics was an iterative procedure, starting first with forecasts using the new model that were initialized from operational analyses. Eventually, the error statistics were estimated from a sample of over 200 cases, spanning all seasons and therefore smoothing out the seasonal dependence, from fully cycled forecasts using the new system. Utilizing the difference between the 24- and 48-h forecasts valid at the same time from this sample, the estimates for the balance projection matrices, standard deviations, and correlation length scales were computed. Although the new statistics are qualitatively similar to those shown in WPP02, some differences were found due to the much higher horizontal resolution and larger sample size. Likewise, the normalization necessary for the new moisture variable resulted in a much different relative humidity background error standard deviation structure.

The background error statistics computed by the NMC method are an estimate and only provide a starting point for a realistic background error applicable to a 6-h forecast. As such, it is necessary to make the proper adjustments to the statistically derived background error standard deviations and length scales. The multivariate correlation matrices were the only components of the background error statistics that were not tuned. The tuning of the background error was accomplished through running single-observation experiments, single-analysis cases, and several fully cycled runs. The goal was to obtain observational fits of the GSI analyses similar to those fits to SSI-generated analyses, while keeping in mind that the structure of the increment would be represented more locally for the GSI analyses. This was generally achieved despite the fact that the analysis variable definitions and structure functions of the background error terms for the GSI and SSI are very different. In general, this tuning resulted in approximately a 30%–40% reduction in the standard deviations and length scales estimated by the NMC method, very similar to the reduction applied to the error standard deviations for the SSI.

The use of recursive filters for the application of the background error in GSI is the same as is outlined in WPP02 and Purser et al. (2003a). Unlike in WPP02, however, where a single horizontal length scale was chosen, three sets of horizontal scales are applied in order to achieve the “fat tailed” feature in the spectrum of the error covariance. All three horizontal scales were derived from statistically estimated horizontal length scales from the NMC method, and applied through a linear combination of three recursive filters, one for each length scale. As in WPP02, the horizontal length scales vary with model level, latitude, and variable. Although it is possible to introduce nonseparability by estimating a different vertical scale for each horizontal scale, testing showed this did not improve the results. As such, a single vertical recursive filter is applied for the vertical correlations (like other aspects of the background error, the vertical length scales are a function of the latitude and height), applicable to all three horizontal scales. The original tests of the global GSI at low resolution in WPP02 required the horizontal background error for the streamfunction and velocity potential to be applied in spectral space at the upper levels of the model, due to the insufficient resolution of the recursive filter lookup table used at the time. However, this has since been resolved and all of the variables had the same application of the background errors at all levels, based solely on the length scales derived from the NMC method.

The GSI solution to the analysis problem is solved iteratively using a conjugate-gradient minimization algorithm, where the solution is preconditioned by the background error covariance. The version that was implemented uses two outer loops consisting of 100 iterations each, whereas the operational SSI, a less efficient code than the GSI (which has been made more efficient through both improved software engineering as well as the removal of spectral transforms in the inner loop), was using two outer loops of 150 iterations at the time it was replaced. The second outer loop consists of an update using the solution of the minimization after 100 iterations, and then reprocessing all of the forward, nonlinear observation operators (including the radiative transfer model) and quality control. If using the same number of processors and using the same number of iterations (in practice, they were reduced from 150 to 100), the GSI results in a reduction of 20% in wall-clock time and 14% memory usage relative to the SSI.

Implementation of the GSI as part of the GDAS cycle occurred along with changes to the Global Forecast System (GFS) model. The version of the GFS model implemented with the package is a T382 spectral model with 64 levels, defined using a hybrid sigma-pressure coordinate, which replaced the terrain-following sigma coordinate. This version of the model also included a change to modularize the radiation scheme (the performance of which was shown to be extremely close to the operational version). Independent tests showed that the new model resulted in changes to the upper troposphere and lower stratosphere in the extratropics where the new vertical coordinate was designed to improve forecasts of the midlatitude jet stream. The analysis is created on the model vertical levels, but on a horizontal Gaussian–linear grid that corresponds to the T382 spectral definition, or 768 × 384 grid points. There is an interface that is part of the GSI system that takes care of the necessary variable transforms, including the spectral to grid transforms (or vice versa), upon reading and writing the forecast and analysis files, respectively. The same horizontal grid was used in computing the latitude-dependent background error statistics. The assimilation of global positioning system (GPS) radio occultation data (Cucurull et al. 2007) and Microwave Humidity Sounder (MHS) brightness temperatures from the National Oceanic and Atmospheric Administration-18 (NOAA-18) polar-orbiting satellite was also turned on with this implementation.

3. Analysis comparison

a. Single-observation tests

As shown in previous work (WPP02; Derber et al. 1991), single-observation tests can be used to demonstrate the impacts of the background error, including the structure functions and multivariate correlations. To evaluate the impacts of going from the SSI definition of the background error to the GSI, many single-observation tests were carried out. Figure 2 shows a comparison of the zonal wind analysis solution for a single zonal wind observation at 300 hPa with a difference from the background as well as an observation error of 1 m s−1, for the SSI and GSI, respectively. Despite the two systems having different variable definitions for the wind, balance operators, and background error definitions, the resultant analysis increments are similar. The GSI wind increment (Fig. 2b) exhibits a smaller spatial scale than that of its SSI counterpart (Fig. 2a).

Fig. 2.

Zonal wind increment (analysis minus background) at 300 hPa (CI is 0.1 m s−1, negative values are dashed, zero contour is omitted) resulting from a single 300-hPa zonal wind observation at 45°N, 90°W with a 1 m s−1 residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 2.

Zonal wind increment (analysis minus background) at 300 hPa (CI is 0.1 m s−1, negative values are dashed, zero contour is omitted) resulting from a single 300-hPa zonal wind observation at 45°N, 90°W with a 1 m s−1 residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Through their respective multivariate variable definition and balance operators, both the GSI and SSI induce a temperature increment from the same zonal wind observation (Fig. 3). A cross section shows that both systems induce a north–south temperature gradient from the zonal wind observation, consistent with a mass–wind balance relationship in the extratropical troposphere. The horizontal spatial scale in the gridpoint system (Fig. 3b) is much smaller, consistent with the smaller spatial scale of the corresponding zonal wind increment. The vertical distribution of the resultant temperature increment is different between the two systems, likely because of the different vertical covariance operators and the different multivariate coupling between the mass and wind, as well as a dynamic constraint applied in the GSI. The SSI tends to project increments on a much deeper vertical scale, a consequence of using empirical orthogonal functions to define the vertical covariances. The application of a vertical recursive filter in the gridpoint system allows for the increment to remain more localized. Additionally, the magnitude of the induced temperature increment is less for the GSI, a consequence of the multivariate coupling between the wind and temperature and the different balance operators applied in the two systems.

Fig. 3.

Cross section taken at 90°W of the temperature increment (CI is 0.01 K, negative values are dashed, zero contour is omitted) resulting from a single 300-hPa zonal wind observation at 45°N, 90°W with a 1 m s−1 residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 3.

Cross section taken at 90°W of the temperature increment (CI is 0.01 K, negative values are dashed, zero contour is omitted) resulting from a single 300-hPa zonal wind observation at 45°N, 90°W with a 1 m s−1 residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

The resultant analysis comparison for single temperature observation experiment (with a 1-K difference from the background and observation error) is also revealing. The smaller horizontal and vertical scales in the gridpoint system for the temperature observation and subsequent increment can be seen in Figs. 4 and 5. The magnitude of the temperature increment is also reduced in the GSI, an indication of the smaller values for the background error variance for temperature. Additionally, the choice of normalized relative humidity as the moisture variable allows the analysis to induce a flow-dependent moisture increment from the single temperature observation (Fig. 1), without making large changes to the relative humidity.

Fig. 4.

Temperature increment at 850 hPa (CI is 0.05 K, negative values are dashed, zero contour is omitted) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 4.

Temperature increment at 850 hPa (CI is 0.05 K, negative values are dashed, zero contour is omitted) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 5.

Cross section taken at 90°W of the temperature increment (CI is 0.05 K, negative values are dashed, zero contour is omitted) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 5.

Cross section taken at 90°W of the temperature increment (CI is 0.05 K, negative values are dashed, zero contour is omitted) resulting from a single 850-hPa temperature observation at 45°N, 90°W with a 1-K residual and observation error for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

These single-observation tests were performed using the new hybrid (sigma pressure) vertical coordinate definition in the GSI, whereas the SSI experiments were performed using the sigma-coordinate definition. The use of new generalized coordinates is not easily done using the SSI system without significant coding changes and careful recalculation and evaluation of the background error statistics (particularly because the SSI uses vertical empirical orthogonal functions for the vertical correlations). Despite the many obvious differences in the formulation and application processes, the GSI system produces analysis increments that are qualitatively similar to those of the SSI, with the biggest differences being the effective horizontal and vertical scales. The GSI increments have been made to be as well balanced as the SSI increments, through the addition of the tangent linear normal mode constraint to complement the multivariate balance definition (Kleist et al. 2009).

b. Single-analysis impact

To more fully visualize the impacts of the new gridpoint system, the GSI and SSI were run for a single analysis using the same background, observation set, and satellite bias correction. All of the files that went into the creation of the operational analysis at 1200 UTC on 20 April 2007 were retrieved and formatted to complete a direct comparison of the old versus new systems. This included a reformatting of the forecast files into the new vertical coordinate (interpolation from sigma to hybrid sigma pressure). Just as in the single-observation tests, differences were expected due to the background error and balance operator definitions. However, many other factors contributed to differences when including the entire observational dataset because of enhancements within the GSI, including improved quality control and observational data handling, and an improved minimization algorithm, as well as better forward models such as an updated radiative transfer model (Han et al. 2006).

Despite the many differences between the two systems, the resultant 250-hPa zonal wind analysis increment for the aforementioned case contain many features that are quite similar (e.g., the large dipole over the central Atlantic; see Fig. 6). Though only a snapshot, this result was typical of other test cases. As was shown in the single-observation tests, the spatial scale of the increment is smaller in the gridpoint system, though the increment structure and magnitude are comparable between the two. However, this is not necessarily the case for the analysis increments of temperature (not shown) and surface pressure (Fig. 7) in the tropics. This is consistent with the new multivariate variable definition and balance operator in the gridpoint system, which does much less in forcing an incremental mass–wind relationship in the tropics than the spectral system. Additionally, the estimates for the standard deviation and horizontal length scale are substantially smaller in the tropics for unbalances surface pressure (relative to the spectral-based representation used in the SSI). It is worth noting that for these individual cases, the resultant fit of the observations to the analysis for each system is comparable (Table 1). The largest discrepancies in the fit of the observations for this case can be seen in the surface pressure [likely the result of the tangent linear normal mode constraint; Kleist et al. (2009)] and relative humidity observations.

Fig. 6.

Zonal wind increment at 250 hPa (CI is 2 m s−1, negative values are dashed, zero contour is omitted) resulting from the assimilation of all operationally available observations within a 6-h time window for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 6.

Zonal wind increment at 250 hPa (CI is 2 m s−1, negative values are dashed, zero contour is omitted) resulting from the assimilation of all operationally available observations within a 6-h time window for (a) SSI and (b) GSI analyses valid at 1200 UTC 20 Apr 2007.

Fig. 7.

As in Fig. 6, but for the surface pressure increment (CI is 0.5 hPa, negative values are dashed, zero contour is omitted).

Fig. 7.

As in Fig. 6, but for the surface pressure increment (CI is 0.5 hPa, negative values are dashed, zero contour is omitted).

Table 1.

The RMS of the fit of conventional observations to the SSI and GSI analyses (OA) after the assimilation of all operationally available observations within the 6-h time window using the same 6-h forecast as the background valid at 1200 UTC 20 Apr 2007.

The RMS of the fit of conventional observations to the SSI and GSI analyses (O − A) after the assimilation of all operationally available observations within the 6-h time window using the same 6-h forecast as the background valid at 1200 UTC 20 Apr 2007.
The RMS of the fit of conventional observations to the SSI and GSI analyses (O − A) after the assimilation of all operationally available observations within the 6-h time window using the same 6-h forecast as the background valid at 1200 UTC 20 Apr 2007.

4. Preimplementation testing

Although the new system appeared to generate reasonable analyses, it was important to test the impacts of the new gridpoint system, the assimilation of new observations, and the changes to the forecast model on the quality of the model forecasts relative to the operational system prior to implementation. Due to the nature of working with an operational system, it was necessary to perform the evaluation by running the suite to mimic operations as closely as possible, at full resolution, with limited available resources. As is often the case, this required developers to bundle together many changes into a single package for preoperational testing. For example, the change in the model vertical coordinate could not have made done without the GSI development (as no resources were put forth to generalize the SSI for other vertical coordinates). Rather than performing two sets of preimplementation experiments, the model changes were bundled with the analysis package to cut the resource requirements in half. Finally, there were many improvements that were developed specifically for the GSI (such as an improved minimization algorithm, enhancements to the data selection and quality control processes, and the aforementioned dynamic balance operator), making it difficult to design truly clean control and experimental runs.

Because the analysis system is such a critical part of the GDAS, and the nature of the changes was deemed significant, it was decided that extensive retrospective and real-time testing would be necessary. Three retrospective test periods were evaluated including the 2005 hurricane season (15 June–5 November), the 2006 hurricane season (30 July–5 November), and the 2006–07 Northern Hemisphere winter (24 October–5 February). All of these retrospective tests included full cycling of the late data cutoff (GDAS) analysis, which also includes the 9-h forecast. For the winter retrospective, the early analysis (GFS) and subsequent 15-day free forecast were only created for the 0000 UTC cycle to minimize resource usage (for all free forecasts, a reduction in resolution to T190 occurs 7.5 days into the simulation, just as in operations). For the hurricane retrospectives, it was necessary to rerun the Geophysical Fluid Dynamics Laboratory (GFDL) hurricane model from the new GFS analyses and boundary conditions for 0000 and 1200 UTC. As such, an additional early analysis and 84-h free forecasts were created at 1200 UTC, with the full 15-day free forecast being run at 0000 UTC. In total, the retrospective testing required creating analyses at four analysis cycles per day (at 0000, 0600, 1200, and 1800 UTC) for 330 simulated days, each of which had a 0000 UTC 15-day free forecast at full operational resolution. Additionally, a real-time parallel that mimicked operations as closely as possible was performed and evaluated prior to implementation (i.e., creating all early and late analyses and forecasts for all four cycles). This parallel was started on 1 March 2007 and was run continuously in real time until after the actual implementation occurred on 1 May 2007. For each of these experiments, a new set of satellite bias correction coefficients were calculated from short, cycled experiments that were completed prior to each retrospective evaluation period. Due to the availability of the data, only those experiments covering dates after 1 November 2006 assimilated GPS radio occultation observations.

A thorough evaluation of all retrospective and parallel experiments was completed relative to the operational system. It is worth mentioning again that due to constraints on the available resources, the package that was tested and implemented included changes to the analysis (including data, as well as options that could have been implemented into the SSI) and forecast model. This of course makes it difficult to interpret the results, as there is not a clean control to compare the results against. Although each component of the new package was tested individually for shorter time periods, it is difficult to attribute improvements or degradations found to any individual part of the new system.

All of the grid-based verification was done by comparing the forecasts from each system against their own respective analyses. The new system performed better in the Northern Hemisphere and comparably in the Southern Hemisphere when comparing the average 500-hPa anomaly correlation scores for the retrospective period covering 15 November 2006–1 February 2007 (Fig. 8). The improvements in the Northern Hemisphere scores were just on the edge of being significant at the 95% confidence level. Though the details were often different for each of the retrospective periods, the new system was always comparable in terms of the anomaly correlation scores, often performing better in one hemisphere and slightly worse or equally well in the other, relative to the operational SSI-based system. Similar behavior was noted for extratropical geopotential height root-mean-square errors at both 1000 and 500 hPa (not shown).

Fig. 8.

The average 500-hPa geopotential height anomaly correlation scores by forecast day for SSI (black, operational) and GSI (gray, experimental) forecasts verifying daily from 15 Nov 2006 through 1 Feb 2007 in the (a) Northern and (b) Southern Hemispheres. The error bars represent the significance of the difference between the two runs at the 95% confidence level. Forecasts for each run were initialized at 0000 UTC and verified against their own analyses, respectively.

Fig. 8.

The average 500-hPa geopotential height anomaly correlation scores by forecast day for SSI (black, operational) and GSI (gray, experimental) forecasts verifying daily from 15 Nov 2006 through 1 Feb 2007 in the (a) Northern and (b) Southern Hemispheres. The error bars represent the significance of the difference between the two runs at the 95% confidence level. Forecasts for each run were initialized at 0000 UTC and verified against their own analyses, respectively.

In the tropics, the new system showed marginal improvements in the RMS vector wind error at 200 hPa (Fig. 9a) and substantial improvements at 850 hPa (Fig. 9b) for the 15 November 2006–1 February 2007 retrospective experiment. This improvement was consistent for all test periods, with the reduction in the RMS vector wind error at 850 hPa being approximately 10%–15%. Based on these experimental results as well as what was shown in WPP02, it would appear that the upgraded package of changes including the GSI as part of GDAS has yielded an improved tropical circulation, both in the analysis and forecasts. As pointed out in WPP02, a smoother solution has an unfair advantage when vector RMS is used as the verification metric, especially for tropical winds where there is sizable analysis uncertainty. However, the tropical analysis increments for wind are of similar magnitude in the new system (Fig. 6), and although these increments have smaller spatial scales, they are not notably smoother. A further evaluation of the fit to the observations showed that the new system fit as closely (for analyses) or closer (for the forecasts) to the tropical observations than did the operational package (not shown).

Fig. 9.

As in Fig. 8, but for the wind vector RMSE (m s−1) in the tropics at (a) 200 and (b) 850 hPa.

Fig. 9.

As in Fig. 8, but for the wind vector RMSE (m s−1) in the tropics at (a) 200 and (b) 850 hPa.

The new system consistently reduced the short-term forecast error and bias, specifically for 6-h forecasts. This was especially true for the 500-hPa geopotential height bias, where the operational system exhibited larger biases over the Southern Ocean, Central America and the Gulf of Mexico, and parts of Asia and Europe (Fig. 10a), with some of the average biases exceeding 15 m. This substantial reduction in the short-term forecast bias (Fig. 10b) suggests that the new analyses may be much more compatible with the forecast model, possibly through an improvement of the mass–wind relationship. It is also possible that the SSI-based system was projecting some of the observational information on much too large of a scale, through the spectral representation of the background error. However, it is nearly impossible to attribute this significant improvement to any single component of the package. The improvement to the forecast at short lead times was also evident for 24-h forecasts (Fig. 11).

Fig. 10.

The average 500-hPa geopotential 6-h forecast error (CI is 2.5 m) from fully cycled systems using the (a) SSI (operations) and (b) GSI (parallel) for the period spanning 0000 UTC 1 Aug–1200 UTC 30 Aug 2007. The mean is computed for all four cycles per day with all forecasts verified against their own analyses, respectively.

Fig. 10.

The average 500-hPa geopotential 6-h forecast error (CI is 2.5 m) from fully cycled systems using the (a) SSI (operations) and (b) GSI (parallel) for the period spanning 0000 UTC 1 Aug–1200 UTC 30 Aug 2007. The mean is computed for all four cycles per day with all forecasts verified against their own analyses, respectively.

Fig. 11.

As in Fig. 10, but for the 24-h forecast error for the period spanning 1 Aug–15 Oct 2007. The mean is computed only from the 0000 UTC forecasts with all forecasts verified against their own analyses, respectively.

Fig. 11.

As in Fig. 10, but for the 24-h forecast error for the period spanning 1 Aug–15 Oct 2007. The mean is computed only from the 0000 UTC forecasts with all forecasts verified against their own analyses, respectively.

Another notable improvement was a significant reduction in the low bias in the short-term forecasts of the East Asian jet (Fig. 12), a longstanding problem in the old system. This finding was independently verified by looking at the fits of the forecasts to conventional observations such as radiosondes. This difference may be primarily attributable to the change in the model vertical coordinates rather than the GSI. There was also improvement found in the fit of the 1- and 2-day stratospheric wind and temperature forecasts to observations in all regions, showing the benefits of the new hybrid vertical coordinate in the forecast model as well as the inclusion of the GSI. A further comparison of the analyses to observations showed that the new system was often quite similar to the operational system, with consistent improvement seen in the fit of the global analyzed moisture to the radiosonde observations (not shown).

Fig. 12.

As in Fig. 11, but for the 200-hPa wind speed 48-h forecast error (CI is 0.5 m s−1).

Fig. 12.

As in Fig. 11, but for the 200-hPa wind speed 48-h forecast error (CI is 0.5 m s−1).

In terms of precipitation forecasts, it was found that the new system exhibited superior skill for nearly all lead times when compared to the operational system over the continental United States (Figs. 13 and 14) for the period covering 15 November 2006–1 February 2007. For convective periods such as in the Northern Hemisphere summer, the new system noticeably reduced the high forecast bias while improving the equitable threat scores (not shown), though, in general, precipitation forecasts from the new system appeared to be qualitatively very similar when compared to operational forecasts. Perhaps some of the improvement may be attributed to the reduced spinup or spindown times in the new system, in addition to the better fits of the analyzed moisture to the observations.

Fig. 13.

Verification of precipitation relative to observations over North America as a function of threshold amounts for 1-day forecasts (average of 12–36-h forecasts) verifying between 15 Nov 2006 and 1 Feb 2007 for systems using the SSI (black, operations) and GSI (gray, parallel). Verification includes the (a) equitable threat score and (b) bias skill score. The gray numbers above (b) indicate the numbers of cases that went into the verification of each threshold category.

Fig. 13.

Verification of precipitation relative to observations over North America as a function of threshold amounts for 1-day forecasts (average of 12–36-h forecasts) verifying between 15 Nov 2006 and 1 Feb 2007 for systems using the SSI (black, operations) and GSI (gray, parallel). Verification includes the (a) equitable threat score and (b) bias skill score. The gray numbers above (b) indicate the numbers of cases that went into the verification of each threshold category.

Fig. 14.

As in Fig. 13, but for the 2-day forecasts (average of 36–60-h forecasts).

Fig. 14.

As in Fig. 13, but for the 2-day forecasts (average of 36–60-h forecasts).

Significant effort was put into assessing the impacts of the new system on tropical cyclone forecasting. The GFS model is itself an important piece of forecast guidance for tropical cyclones at the National Hurricane Center, but it is also used to create initial and boundary conditions for the GFDL nested hurricane model (Bender et al. 2007), another important forecast guidance tool. As such, tropical cyclone forecasts from both the GFS model, as well as the GFDL reruns initialized from those very same GFS forecasts were evaluated. It was found that for the 2005 and 2006 Atlantic hurricane seasons, the new system led to a reduction in the average track error for both the GFS (Fig. 15) and GFDL (not shown) models compared to their operational counterparts. The reduction in track error seemed to be especially significant for lead times beyond 48 h. This reduction in the track error is consistent with the improvement seen in the tropical circulation and reduced vector wind RMS error in the new system. However, the results for the eastern Pacific storms over this same period were more mixed, with the new system having similar track errors to the operational system.

Fig. 15.

The average GFS tropical cyclone track error (n mi) out to 120 h from systems using the SSI (black, operations) and GSI (gray, parallel) for the combined 2005 and 2006 hurricane seasons in both the (a) Atlantic basin and (b) east Pacific basin. The average track error was computed from a homogeneous sample of cases, with the number of cases for each lead time identified below each forecast hour.

Fig. 15.

The average GFS tropical cyclone track error (n mi) out to 120 h from systems using the SSI (black, operations) and GSI (gray, parallel) for the combined 2005 and 2006 hurricane seasons in both the (a) Atlantic basin and (b) east Pacific basin. The average track error was computed from a homogeneous sample of cases, with the number of cases for each lead time identified below each forecast hour.

5. Summary and discussion

After extensive development and testing, the GSI was implemented as part of the NCEP GDAS cycle on 1 May 2007, along with changes to the GFS forecast model. Retrospective and real-time testing showed that the new system was competitive with the SSI-based system and often exhibited superior performance. The most noteworthy improvements included a reduction in the bias of forecasts of the East Asian jet, smaller short-range 500-hPa geopotential height and sea level pressure biases, better precipitation forecasts over the continental United States, and more skillful tropical cyclone track forecasts.

Shortly after the implementation described in this manuscript, a minor data upgrade occurred to the GDAS on 26 May 2007 that included improvements to the quality control and observation errors for satellite-derived winds. Additionally, High Resolution Infrared Radiation Sounder (HIRS), MHS, and Advanced Microwave Sounding Unit-A (AMSU-A) brightness temperatures from the European Organization for the Exploitation of Meteorological Satellites’ METOP-A, which flew in an orbit similar to that of NOAA-17, began to be assimilated. Finally, the assimilation of the Geostationary Operational Environmental Satellite (GOES) averaged 5 × 5 brightness temperatures were replaced with the assimilation of detector-specific nonaveraged (i.e., every field of view) observations. The impacts of this data upgrade was a small, but positive improvement to the analyses and forecasts (not shown).

All of the improvements gained with the new system were achieved without exercising many of the options that the GSI was originally developed to carry out, a very encouraging development as this implementation was only a first step aimed at substantially improving global and regional analyses at NCEP. The background error structure used in the implementation was designed to be similar to those defined for the SSI, homogeneous in the zonal direction, and isotropic. With the analysis now formulated in grid space, more realistic background error structures (i.e., anisotropic, inhomogeneous, flow dependent) can be developed and tested.

The modular characteristics and inclusion of a generalized vertical coordinate structure within the GSI has allowed for a streamlining of the development, improving the data assimilation between various groups within NCEP. The GSI has already been implemented for many different applications at NCEP, including within the North American Mesoscale Model Data Assimilation system (NDAS, implemented 20 June 2006) and the Real-Time Mesoscale Analysis (RTMA), and to compute sea surface temperature retrievals from the Marine Modeling and Analysis branch at NCEP’s Environmental Modeling Center (EMC). Great benefit has also been achieved through collaboration with groups outside of NCEP, especially from the Global Modeling and Assimilation Office at the National Aeronautics and Space Administration (NASA), where work is on going to extend the GSI for use in 4DVAR assimilation.

Now that the GSI has been implemented as part of the GDAS at NCEP, further improvements can be expected in the near future through the development and testing of more advanced data assimilation techniques and the assimilation of new observations. There are a few options that have already been added and recently implemented as part of the GSI component of GDAS. One of these options is the inclusion of an efficient method of employing flow-dependent background error variance structure. Additional changes have allowed for the better use of conventional observation data, through the addition of variational quality control. Finally, through the Joint Center for Satellite Data Assimilation, the assimilation of the next generation of satellite observations using GSI can be rapidly tested and implemented.

Acknowledgments

The authors thank our many collaborators at EMC, in particular R. James Purser, whose work on recursive filters made this development possible. Additionally, we wish to thank our many partners within the Joint Center for Satellite Data Assimilation, especially James Jung and Paul van Delst for their many contributions. Many useful comments and suggestions were provided by David Groff, George Gayno, and three anonymous reviewers, which greatly improved the manuscript and are much appreciated.

REFERENCES

REFERENCES
Bannister
,
R. N.
,
2008a
:
A review of forecast error covariance statistics in atmospheric variational data assimilation. I: Characteristics and measurements of forecast error covariances.
Quart. J. Roy. Meteor. Soc.
,
134
,
1951
1970
.
Bannister
,
R. N.
,
2008b
:
A review of forecast error covariance statistics in atmospheric variational data assimilation. II: Modelling the forecast error covariance statistics.
Quart. J. Roy. Meteor. Soc.
,
134
,
1971
1996
.
Bender
,
M. A.
,
I.
Ginnis
,
R.
Tuleya
,
B.
Thomas
, and
T.
Marchok
,
2007
:
The operational GFDL coupled hurricane–ocean prediction system and a summary of its performance.
Mon. Wea. Rev.
,
135
,
3965
3989
.
Cucurull
,
L.
,
J. C.
Derber
,
R.
Treadon
, and
R. J.
Purser
,
2007
:
Assimilation of global positioning system radio occultation observations into NCEP’s Global Data Assimilation System.
Mon. Wea. Rev.
,
135
,
3174
3193
.
Daley
,
R.
, and
E.
Barker
,
2001
:
NAVDAS: Formulation and diagnostics.
Mon. Wea. Rev.
,
129
,
869
883
.
Dee
,
D. P.
, and
A. M.
da Silva
,
2003
:
The choice of variable for atmospheric moisture analysis.
Mon. Wea. Rev.
,
131
,
155
171
.
De Pondeca
,
M.
, and
Coauthors
,
2007
:
The status of the real time mesoscale analysis at NCEP.
Preprints, 22nd Conf. on Weather Analysis and Forecasting/18th Conf. on Numerical Weather Prediction, Park City, UT, Amer. Meteor. Soc., 4A.5. [Available online at http://ams.confex.com/ams/pdfpapers/124364.pdf]
.
Derber
,
J. C.
, and
F.
Bouttier
,
1999
:
A reformulation of the background error covariance in the ECMWF Global Data Assimilation System.
Tellus
,
51A
,
195
221
.
Derber
,
J. C.
,
D. F.
Parrish
, and
S. J.
Lord
,
1991
:
The new global operational analysis system at the National Meteorological Center.
Wea. Forecasting
,
6
,
538
547
.
Devenyi
,
D.
,
S. G.
Benjamin
,
J. M.
Middlecoff
,
T. W.
Schlatter
, and
S. S.
Weygandt
,
2005
:
Gridpoint Statistical Interpolation for Rapid Refresh.
Preprints, 21st Conf. on Weather Analysis and Forecasting/17th Conf. on Numerical Weather Prediction, Washington, DC, Amer. Meteor. Soc., P1.56. [Available online at http://ams.confex.com/ams/pdfpapers/95149.pdf]
.
Devenyi
,
D.
,
S. S.
Weygandt
,
T. W.
Schlatter
,
S. G.
Benjamin
, and
M.
Hu
,
2007
:
Hourly data assimilation with the Gridpoint Statistical Interpolation for Rapid Refresh.
Preprints, 22nd Conf. on Weather Analysis and Forecasting/18th Conf. on Numerical Weather Prediction, Park City, UT, Amer. Meteor. Soc., 4A.2. [Available online at http://ams.confex.com/ams/pdfpapers/124535.pdf]
.
Fisher
,
M.
, and
P.
Courtier
,
1995
:
Estimating the covariance matrices of analysis and forecast error in variational data assimilation.
ECMWF Research Dept. Tech. Memo. 220, 26 pp. [Available from European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG2 9AX, United Kingdom]
.
Han
,
Y.
,
P.
van Delst
,
Q.
Liu
,
F.
Weng
,
B.
Yan
,
R.
Treadon
, and
J.
Derber
,
2006
:
JCSDA Community Radiative Transfer Model (CRTM)—Version 1.
NOAA Tech. Rep. NESDIS 122, 33 pp
.
Holm
,
E.
,
2003
:
Revision of the ECMWF humidity analysis: Construction of a Gaussian control variable.
Proc. Workshop on Humidity Analysis, Reading, United Kingdom, ECMWF/GEWEX, 1–6. [Available from European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG29AX, United Kingdom]
.
Kleist
,
D. T.
,
D. F.
Parrish
,
J. C.
Derber
,
R.
Treadon
,
R. M.
Errico
, and
R.
Yang
,
2009
:
Improving incremental balance in the GSI 3DVAR analysis system.
Mon. Wea. Rev.
,
137
,
1046
1060
.
Lorenc
,
A. C.
, and
Coauthors
,
2000
:
The Met. Office global 3-dimensional variational data assimilation scheme.
Quart. J. Roy. Meteor. Soc.
,
126
,
2991
3012
.
Parrish
,
D. F.
, and
J. C.
Derber
,
1992
:
The National Meteorological Center’s spectral statistical-interpolation system.
Mon. Wea. Rev.
,
120
,
1747
1763
.
Purser
,
R. J.
,
W-S.
Wu
,
D. F.
Parrish
, and
N. M.
Roberts
,
2003a
:
Numerical aspects of the application of recursive filters to variational statistical analysis. Part I: Spatially homogeneous and isotropic Gaussian covariances.
Mon. Wea. Rev.
,
131
,
1524
1535
.
Purser
,
R. J.
,
W-S.
Wu
,
D. F.
Parrish
, and
N. M.
Roberts
,
2003b
:
Numerical aspects of the application of recursive filters to variational statistical analysis. Part II: Spatially inhomogeneous and anisotropic general covariances.
Mon. Wea. Rev.
,
131
,
1536
1548
.
Rabier
,
F.
,
A.
McNally
,
E.
Anderson
,
P.
Courtier
,
P.
Unden
,
J.
Eyre
,
A.
Hollingsworth
, and
F.
Bouttier
,
1998
:
The ECMWF implementation of three-dimensional variational assimilation (3D-Var). II: Structure functions.
Quart. J. Roy. Meteor. Soc.
,
124
,
1809
1829
.
Rabier
,
F.
,
H.
Järvinen
,
E.
Klinker
,
J-F.
Mahfouf
, and
A.
Simmons
,
2000
:
The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplified physics.
Quart. J. Roy. Meteor. Soc.
,
126
,
1143
1170
.
Rawlins
,
F.
,
S. P.
Ballard
,
K. J.
Bovis
,
A. M.
Clayton
,
D.
Li
,
G. W.
Inverarity
,
A. C.
Lorenc
, and
T. J.
Payne
,
2007
:
The Met Office global four-dimensional variational data assimilation scheme.
Quart. J. Roy. Meteor. Soc.
,
133
,
347
362
.
Wu
,
W-S.
,
2005
:
Background error for NCEP’s GSI analysis in regional mode.
Proc. Fourth Int. Symp. on Analysis of Observations in Meteorology and Oceanography, Prague, Czech Republic, WMO, 3A.27
.
Wu
,
W-S.
,
D. F.
Parrish
, and
R. J.
Purser
,
2002
:
Three-dimensional variational analysis with spatially inhomogeneous covariances.
Mon. Wea. Rev.
,
130
,
2905
2916
.

Footnotes

* Additional affiliation: Science Applications International Corporation, Camp Springs, Maryland.

Corresponding author address: Daryl T. Kleist, NOAA Science Center 207, 5200 Auth Rd., Camp Springs, MD 20746-4304. Email: daryl.kleist@noaa.gov