Part III:   How to build an ensemble prediction
system (EPS)?  

1-Dimentional EPS

    Only IC uncertainty is considered in 1-D EPS by perturbing initial conditions. Three basic properties need to be followed in designing a perturbation: Realism, Divergence and Orthogonality. Realism is that the magnitude of perturbation needs to be within the size of realistic analysis error and should exhibit a realistic spectral distribution over spatial scales such as larger uncertainty in smaller-scale waves (difficult to observe) and smaller uncertainty in larger-scale waves (easier to observe). Divergence is that the perturbation needs to have dynamically growing structure leading members to diverge as much as possible during model integration to cover all possible solutions in a model space. Orthogonality is that the perturbation needs to be orthogonal to each other among members to maximize information content contained in an ensemble, which is especially important for a small size ensemble. Initial conditions or states needed to be perturbed include interior, lower-, upper- and lateral-boundary (if limited-area model) conditions. There are currently five or so different categories proposed for perturbing ICs especially for interior states.

    (1) Random Perturbation (Monte Carlo approach): Perturbation is randomly generated based on some kind of statistics (usually a normal distribution) representing typical uncertainty in analysis (such as the average difference between two analyses over a long period of time) (Errico and Baumhefner, 1987; Mullen and Baumhefner, 1994; Du et. al., 1997). Although this type of perturbation represents well the average magnitude of analysis uncertainty (Realism), it is lack of dynamically growing spatial structure and not reflecting “error of the day”. As a result, the perturbation growth rate is low and, therefore, the Divergence in solution among members is usually not ideal. Random Perturbation is usually applied in places where one is unsure which other methods work better.

    (2) Time-Lagged approach: There are “Direct Time-Lagged” and “Scaled Time-Lagged” two kinds. Direct Time-Lagged approach (Hoffman and Kalnay, 1983) directly pulls multiple forecasts which are initiated from different past times but verified at a same time together as an ensemble (a mixture of old and young forecasts). This method views the error of a past forecast at t=0 (current initial time) directly as IC perturbation which should reflect “error of the day” and has dynamically growing structure leading to larger ensemble spread than random perturbation. The advantage of the method is that the generation of perturbation is absolutely free and there is no need to purposely generate IC perturbation fields for the ensemble, which implies that all operational NWP centers have this type of ensemble automatically. However, a main concern is that the quality (magnitude) of perturbation depends on the age of a forecast since forecast quality usually decreases with lead time. To avoid this weakness, past forecast errors are first scaled down by their “ages” (assuming error growth is quasi-linear) at t=0 to have similar magnitude in all perturbations and then added to or subtracted from the current control analysis to create multiple analyses to initiate an ensemble of forecasts (Ebisusaki and Kalnay, 1983; Kalnay, 2003). This modified version is called Scaled Time-Lagged method and can be simply described by the equation (4):

               Initial perturbation = scaling x (time-lagged forecast – current analysis)     (4)

The latter is not only able to control perturbation size but also doubles the ensemble membership by using both addition and subtraction procedures with very little extra computing cost (which is very important in real-time production) in generating perturbation. Note that this Scaled Time-Lagged method is actually already the same in idea and similar in technical procedure to the Breeding method (see the next paragraph). The Time-Lagged approach has been used in many ensemble research and operation such as NCEP operational seasonal ensemble forecast system (Saha et. al., 2006; Hou et. al., 2001; Lu et. al., 2006; Brankovic et. al., 2006; Mittermaier, 2007). A limitation of the Time-Lagged method is that it cannot create an ensemble with large enough membership size since the number of “good” old forecasts available is limited in reality. Otherwise, the ensemble quality will be severely contaminated if too old forecasts are included to have a large size ensemble. By the way, Time-Lagged forecasts are often used as “initial seeds” to cold start an ensemble which uses other perturbation methods such as Breeding.

    (3) Breeding: Different from Scaled Time-Lagged method, Breeding uses two concurrent forecasts (at a recent past time t=-T) rather than a time-lagged forecast and a current analysis to calculate raw perturbation at t=0. The difference is then be scaled down and added to or subtracted from the current control IC (Toth and Kalnay, 1993 and 1997). In this way, one can create as many members as he wants to have a large size ensemble as long as he can create enough initial random seeds or use other approaches (such as borrowing from other existing forecasts like time-lagged ones) in the beginning to initiate or have enough different forecasts to start with (cold start). So, Breeding method can overcome the membership limitation of Scaled Time-Lagged method. Since all the past forecasts used are now concurrent, the scaling factor has no need to be forecast-age based but can be anything related to analysis uncertainty to control perturbation magnitude and/or spatial structure. The perturbation can now be simply described by the equation (5):

               Initial perturbation = scaling x (forecast 1 – forecast 2)                                (5)

Comparing Eqs. (4) and (5), one can also see that the perturbation in Breeding (bred vector) is no longer pure forecast error but the difference between two past forecasts which is a nonlinear extension of Lyapunov vector (Kalnay, 2003). Experience (e.g., NCEP SREF) shows that bred vector becomes mature in structure and leads to large ensemble spread growth after cycling for about two to three days after the cold start. Toth and Kalnay pointed out that the spatial structure of a mature bred vector is not sensitive to scaling period (T) and norm selected and that bred vector reflects analysis error (error of the day) introduced during data assimilation cycle well (Realism). Although the difference between two past forecasts should, theoretically, reflect the error growing structure of the immediate past cycle but not of the forecast period in future (i.e., a looking-backward approach), experiments show that bred vector’s growth rate (Divergence) is quite satisfactory in practice and higher than that by using either Monte Carlo or Scaled Time-Lagged methods (Toth and Kalnay, 1993 and 1997). Since this method is simple with no mathematical simplification or assumption (but using full nonlinear primitive-equation based model) and easy to implement, has little cost in computing power and gives good ensemble spread, it’s widely used and tested by many such as NCEP ensemble systems (Du and Tracton, 2001; Tracton and Kalnay, 1993) and CMA (China Meteorological Administration). However, it’s reported that bred vectors among ensemble members are not orthogonal enough but highly correlated to each other, resulting less optimal in information content contained in an ensemble (Wang and Bishop, 2003; Martin et. al., 2007). As a result of this, the ensemble spread growth (mainly in magnitude but not structure) is closely related to initial amplitude of bred vector. To orthogonize bred vectors, Ensemble Transform (ET) technique is used to make bred vectors more orthogonal to each other by applying a simplex transformation matrix to transform forecast-based perturbations to analysis perturbations (Wei et. al., 2007). Experiment shows that ET technique can improve ensemble performance over classical Breeding method. Therefore, ET has been implemented in the NCEP global ensemble system to enhance the Breeding method (Wei et. al., 2007). Another approach proposed to improve classical Breeding method is called Geometric Breeding which controls spatial correlation of bred vectors among members to make them less correlated to each other (Martin et. al., 2007). Geometric Breeding shows better spread-skill relation than classic Breeding. Since bred vector mainly depicts synoptic scale baroclinic instability but not smaller-scale convective instability (Toth and Kalnay 1993), it is, however, desired to have smaller-scale instabilities included in perturbation for a mesoscale EPS focusing on predicting, say, convective system related heavy precipitation events. Chen et. al. (2003) suggested that differencing two forecasts from a same model but with different versions of a convective scheme (instead of one same version as in classical Breeding) will help to depict convective instability in perturbation and, therefore, enhances ensemble performance in predicting heavy precipitation. On other hand, due to the characteristic that fast, small scale modes quickly saturated, and only slow, large scale modes left during breeding cycle, bred vector is a good candidate to be used in ocean-atmosphere coupled ensemble prediction system which is mainly associated with slow modes (Cai et. al., 2002, Yang et. al. 2006).  Recently, Prof. Eugenia Kalnay (2007) reported that bred vector has capability of predicting weather regime transition.

    (4) Singular Vector (SV): This method first needs to develop a linear version of a nonlinear model (called Tangent Linear Model, TLM) as well as an Adjoint (Errico, 1997) of the TLM. Then over a desired optimal future time period such as 0-48h, TLM is integrated forward in time and then the Adjoint integrated backward in time to find initial sensitive areas to the forecast, say, at t=48h. This “forward and backward” cycling process needs to be iterated many times to obtain leading singular vectors. Then, a linear combination including rescaling and orthogonal rotating is applied to the vectors to construct a desired number of perturbations. With addition and subtraction of the perturbations to and from a control analysis, an ensemble of forecasts can be formed. Unlike the bred vector, the structure of SV is sensitive to the norm used and the cycling period selected (Errico and Vukiceric, 1992; Palmer et. al., 1998). ECMWF chooses total energy as norm and 0-48 h cycling window to calculate SVs in their global ensemble system (Buizza, 1994; Palmer et. al., 1998). Obviously, SV is a looking-forward approach rather than looking-backward in time like Breeding does and mathematically optimizes large perturbation growth and orthogonarity to have large ensemble spread and contain more information at a pre-selected targeted forecast time. SV method is widely used and tested in both research and operation such as ECMWF and Canadian Meteorological Service’s regional ensemble (Li et. al., 2007). A disadvantage of the method is that it is costly in computing because the number of the iteration of this “forward and backward” integration required is usually about 3 times the number of SVs you want to create (e.g., it needs to integrate about 3x50x2=300 times of 48h-forecast to obtain 50 SVs optimizing at 48h lead time. Therefore, SVs have to be calculated in a much reduced model resolution (comparing to actual forecast model resolution) to save computing time in operation. Another disadvantage is that a finite forecast lead time has to be pre-defined at which SVs are targeted to be optimal in growth. So that SV-based ensemble might not be optimal in performance crossing multiple time ranges such as short and median range at the same time. The linear assumption that perturbation is small enough so that its evolution can be governed by the linear versions (TLM and Adjoint) of a nonlinear model is also a concern in calculating traditional SVs although it retains nonlinear property in some degree by calculating and combining multiple SVs. To overcome the linear assumption, some efforts have been made such as modifying the iteration procedure (Oortwin and Barkmeijer, 1995; Barkmeijer, 1996) and introducing the nonlinear singular vector concept (Mu, 2000) and the conditional nonlinear optimal perturbation (CNOP) method (Mu et. al., 2003; Mu and Duan, 2003; Mu and Zhang, 2006). With simple models, CNOP method showed improved capability of depicting nonlinear features in perturbation comparing to SV method although further studies are needed with full NMP models. Adding moist physics in TLM and Adjoint (Ehrendorfer et. al., 1999) is another step to be closer to reality and showed improved performance. Since SV is mathematically looking forward and focusing on perturbation growing structure at a future time but not related to any immediate past, a reasonable question to ask is that does SV-based perturbation actually reflect error of the day which is introduced by the most recent data assimilation cycle? The following efforts are addressing this kind of concern and yield improved results (Barkmeijer et. al., 1998; Fischer et. al., 1998): e.g, evolving-SV approach by adding the final or evolved singular vector from an immediate previous cycle, say, 48h-period prior to the current model initiation time (analysis time) to the current cycle’s SV perturbation; using analysis error covariance in replacing total energy as norm to calculate SVs; and the application of Kalman filter and so on. It’s interesting to notice that the resulting SV perturbation is closer to Lyapunov vector or Bred vector by either switching to evolving-SV or using analysis error covariance as norm (Kalnay, 2003; Reynolds and Errico, 1999).

    (5) Coupling with Data Assimilation: The simplest version of this method is directly using multiple analyses available to initiate an ensemble of forecasts (Tracton et. al., 1998; Grimit and Mass, 2002). However, the number of available analyses is quite limited, which will restrict ensemble size. By perturbing observations, Houtekamer et. al. (1996) and Houtekamer and Mitchell (1998) purposely generated multiple analyses to initiate their global ensemble system, which shed lights to a new IC perturbation generating approach for ensemble so called Ensemble Transform Kalman Filter (ETKF) approach (Anderson, 1996). The ETKF approach was further explored in detail by Wang and Bishop (2003), Wang et. al. (2004) and Wei et. al. (2006) for generating ensemble. In their studies, ETKF transforms forecast perturbations into analysis perturbations by multiplying a transformation matrix. Using observational information, the magnitude of the analysis perturbations was adjusted before the perturbations are added to a control analysis to initiate an ensemble of forecasts. The transformation matrix used can also guarantee all perturbations to be orthogonal to each other, a desired property for ensemble forecasting. Although ETKF was not used in the data assimilation procedure to directly output multiple analyses in their studies, ETKF itself is an ensemble-based data assimilation technique (Tippett et. al., 2003; Anderson, 2001; Whitaker and Hamill, 2002; Ott et. al., 2004; Szunyogh et. al., 2004; Hamill, 2006; Zhang et. al., 2004; Wang et. al., 2007). Therefore, ETKF technique has a potential to directly link ensemble prediction and data assimilation (DA) into one unified procedure in an NWP system: ensemble forecast variance provides background error covariance information for DA, while DA provides an ensemble of analyses to initiate an ensemble of forecasts. In such a coupled system, not only EPS can be improved by having more realistic IC perturbations reflecting true error of the day in the analysis, but also the quality of analysis improved by using flow-dependent background error information (Hamill, 2006; Zhang, 2005). Therefore, it’s believed that ETKF method has a great potential. This method has been explored to be applied to operation such as regional ensemble systems at UK Met-Office (Mylne, personal communication) and U.S. Navy and Air Force (McLay et. al., 2007).  Active research is, however, still going on in this area. For example, it might need a large ensemble size to provide adequate error covariance information for DA, which is very costly in production. A dual-resolution idea has, therefore, been proposed to run such an ETKF-based system to significantly reduce computing cost (Gao and Xue, 2007; Gao et. al., 2007).

    It is worth pointing out that all the currently existing IC perturbation generating methods do not work well and yield very little ensemble spread in tropics. This is because the current methods are mainly dealing with slower and larger-scale baroclinic instability which dominates middle- and high-latitude atmospheric motion but not faster and smaller-scale barotropic and convective instabilities which dominate tropical atmospheric motion. Therefore, special consideration needs to be researched for tropics. A simple comparison among various methods can be found in Bowler (2006).

    Besides interior initial states, the content of initial condition should also include lower-, upper- and lateral-boundary forcing. The lower boundary forcing is introduced by land and water surface initial parameters such as sea surface temperature, heat and moisture flux, ice and snow cover, soil properties including moisture, temperature and type, surface albedo, roughness and greenness etc.. Out of those, sensitivity of initial soil moisture uncertainty to ensemble prediction in short-range (0-3 days) has been paid particular attention so far. It is found that soil moisture uncertainty plays an important role in convective precipitation during warm season but less important in large-scale precipitation during cool season (Sutton et. al., 2006; Aligo et. al., 2007; Du et. al., 2007a). The sensitivity of some surface variables such as 2-meter temperature to initial soil moisture remains high in both warm and cool seasons although such a sensitivity demonstrates strong diurnal variation related to radiation (much stronger during daytime than nighttime) and geographically preferred regions in particular forecasts (Du et. al., 2007a). Du et. al., 2007a) also reported that the effectiveness of soil moisture perturbation to ensemble forecast depends on perturbation’s spatial structure and magnitude: spatially uniform and larger magnitude perturbations produce larger ensemble spread than spatially random and smaller magnitude ones do. In general, comparing with the sensitivity to atmospheric initial condition uncertainty, the sensitivity to land-surface initial parameter uncertainty seems to be secondary for short-range forecasts, which, of course, could be just opposite for longer-range seasonal forecasts which could be more sensitive to lower-boundary forcing uncertainty. Therefore, it’s recommended that perturbation in surface forcing needs to be combined with other perturbations (such as IC and physics) in order to build a robust ensemble system. As for upper boundary forcing from space such as solar activity or how the upper boundary being treated in a model, not much attention is yet paid to how sensitive a weather forecast will be or how much ensemble spread could be attributed to due to the uncertainty in its initial description. It needs to be studied quantitatively. For a regional ensemble system, lateral boundary condition (LBC) could play a dominant role in defining ensemble spread of many variables (except for precipitation) if model domain is small (Du and Tracton, 1999; Warner et. al., 1997). Therefore, it is recommended that (a) a large enough model domain should be used in limited-area model (LAM) based ensemble system to avoid negative impact from LBCs; and (b) LBCs need to be perturbed too to ensure diverse ensemble solutions. Currently, a common approach in practice is to use different ensemble members from an available global ensemble system as LBCs for different members in a LAM-EPS such as NCEP SREF (Du et. al., 2004). Nutter et. al. (2004a and 2004b) suggested an approach to compensate spread loss due to LBC. How does the inconsistency or consistency in structure between LBC and internal IC perturbation affect ensemble performance is an issue yet needs to be investigated. It’s reported that Spanish Weather Service (INM) has built a LAM-EPS purely based on diversity in LBCs (Garcia-Moya, 2006, personal communication).

Contact  Jun Du