Clim Dyn DOI 10.1007/s00382-012-1302-7

Flow-dependent empirical singular vector with an ensemble Kalman filter data assimilation for El Nino prediction Yoo-Geun Ham • Michele M. Rienecker

Received: 13 June 2011 / Accepted: 18 January 2012 Ó Springer-Verlag 2012

Abstract In this study, a new approach for extracting flow-dependent empirical singular vectors (FESVs) for seasonal prediction using ensemble perturbations obtained from an ensemble Kalman filter (EnKF) assimilation is presented. Due to the short interval between analyses, EnKF perturbations primarily contain instabilities related to fast weather variability. To isolate slower, coupled instabilities that would be more suitable for seasonal prediction, an empirical linear operator for seasonal timescales (i.e. several months) is formulated using a causality hypothesis; then, the most unstable mode from the linear operator is extracted for seasonal time-scales. It is shown that the flow-dependent operator represents nonlinear integration results better than a conventional empirical linear operator static in time. Through 20 years of retrospective seasonal predictions, it is shown that the skill of forecasting equatorial SST anomalies using the FESV is systematically improved over that using Conventional ESV (CESV). For example, the correlation skill of the NINO3 SST index using FESV is higher, by about 0.1, than that of CESV at 8-month leads. In addition, the forecast skill improvement is significant over the locations where the correlation skill of conventional methods is relatively low, indicating that the FESV is effective where the initial uncertainty is large.

Y.-G. Ham (&)  M. M. Rienecker Global Modeling and Assimilation Office, NASA/GSFC Code 610.1, Greenbelt, MD, USA e-mail: [email protected] Y.-G. Ham Goddard Earth Sciences Technology and Research Studies and Investigations, Universities Space Research Association, Baltimore, MD, USA

Keywords Singular vector  Seasonal prediction  Ensemble Kalman filter  El Nino  ENSO

1 Introduction It has been widely accepted that ensemble forecasts using slightly perturbed initial conditions can have a beneficial impact on the forecast skill by means of ensemble averaging and probabilistic forecasts (Molteni and Palmer 1993; Buizza et al. 1998). Hence, one of the aspects crucial for maximizing forecast skill is the generation of optimal initial perturbations (Leith 1974; Toth and Kalnay 1993, 1997; Molteni et al. 1996; Zhu et al. 2002; Buizza et al. 2005; Leutbecher and Palmer 2008). Most existing studies tried to generate fast-growing perturbations that can effectively represent/reduce forecast uncertainty with limited numbers of ensemble members. One well known approach for generating fast-growing perturbations is the singular vector method (Farrell 1989; Palmer et al. 1994; Molteni et al. 1996; Xue et al. 1997a, b; Chen et al. 1997; Fan et al. 2000; Tang et al. 2006). This method was first proposed for weather forecasting at the European Centre for Medium-Range Weather Forecasts (ECMWF). Unstable perturbation patterns were defined from the singular vectors of the linear propagator of the primitive equations. This method was also adopted for intermediate coupled models (Moore and Kleeman 1996, 1997a, b, 1998). Although the singular vector method is a powerful tool and used operationally in several ensemble prediction systems, it has limitations for Coupled General Circulation Models (CGCMs) because a linearized version of the model is required. To overcome these deficiencies, an alternative method for extracting singular vectors of an empirical linear

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

operator has been developed (Tziperman et al. 2008; Hawkins and Sutton 2009; Kug et al. 2010). They extracted singular vectors empirically using an historical integration, obviating the need for a linearized version of the model or additional computations. Kug et al. (2010) demonstrated that the first leading singular vector from the empirical linear operator represents, to a large extent, the fastgrowing mode in the nonlinear integration. With the aid of the fast-growing property of the empirical singular vector, they also found that with the optimal perturbations the forecast skill for equatorial SST anomalies improved at most lead times and regions. Similarly, Tziperman et al. (2008) and Hawkins and Sutton (2009) examined the unstable perturbations associated with decadal predictability of Atlantic Ocean anomalies using a CGCM. With this approach, they diagnosed the structure of the optimal initial perturbation related to the Atlantic Meridional Overturning Circulation (AMOC). In addition, they demonstrated that the evolution of temperature and salinity fields in the Atlantic can be effectively described by an empirically-derived linear system. However, these studies capture a time-fixed fast growing mode, so their empirical singular vectors will sometimes be different from the actual fast-growing error mode, which could be highly flow-dependent (Toth and Kalnay 1997; Hamill et al. 2003; Buehner 2005; Yang et al. 2006). Since the time-varying component of error statistics are significant, a new method is needed to generate optimal flowdependent perturbations. To formulate the linear operator for extracting flowdependent empirical singular vectors, multiple realizations of the state are needed. For example, to extract empirical singular vectors for winter 1997, ensemble integrations of the forecast model with perturbed initial conditions for winter 1997 are needed. This strategy has been applied in an examination of the dependency of the optimal growth pattern on the seasonal cycle and the ENSO phase and the implication for ENSO predictability (e.g. Lorenz 1965; Xue et al. 1994, 1997a, b; Chen et al. 1997; Fan et al. 2000; Kleeman et al. 2003; Tang et al. 2006). Another technique is the so-called breeding method (Toth and Kalnay 1993, 1997; Corrazza et al. 2003), which performs a sequence of perturbed free integrations, rescaling perturbations to extract those that are fast-growing. Toth and Kalnay (1993) first suggested the breeding method for the fastest-growing perturbation in a weather forecast system. Later studies have applied bred vectors for seasonal to interannual time scales using an intermediate coupled model (Cai et al. 2003) and CGCMs (Yang et al. 2006, 2008) for seasonal forecast ensembles. However, these methods require additional computations to obtain ensemble perturbations. These concepts can be merged into an ensemble data assimilation technique (Houtekamer and Mitchell 1998;

123

Evensen and Van Leeuwen 2000; Keppenne and Rienecker 2002; Evensen 2003; Zhang et al. 2005; Leeuwenburgh 2007; Keppenne et al. 2008). The main assumption in this technique is that error statistics can be represented by an ensemble of possible model states (Evensen 1994). One of the well known approaches of the ensemble data assimilation approach is the ensemble Kalman filter (EnKF), which estimates the flow-dependent error statistics from an evolving ensemble of the forecast model. Fast-growing perturbations that would be optimal can be extracted from EnKF system when the ensemble size is large enough to span all possible directions. However, in practical applications of the EnKF, the ensemble size is limited. In addition, due to the short interval between analyses, EnKF perturbations contain instabilities related to the faster weather or intra-seasonal variability, and do not successfully isolate the slowly varying coupled instability relevant for seasonal forecasts. For example, the analysis time interval for most operational ocean data assimilation systems for seasonal forecasts is a few days (for example, 2 days for the NASA/GMAO EnKF system; Keppenne et al. 2008). On the other hand, a physically meaningful choice of rescaling period for bred vectors to capture the slowly varying coupled instability and saturate the growth of weather signals would be longer than 2 weeks (Pen˜a and Kalnay 2004; Yang et al. 2006). This implies that both ensemble perturbations from the ensemble data assimilation technique and empirical singular vectors extracted from these ensemble perturbations are not optimal perturbations for the seasonal forecasts. In this study, a method to extract fast-growing perturbations optimized for the seasonal forecasts from an EnKF system is developed. This paper is organized as follows. In Sect. 2, the forecast model used in this study and the data assimilation experiments are described. In Sect. 3, a new algorithm to extract flow-dependent empirical singular vectors for seasonal prediction within the EnKF framework, and its application to the forecast model, is described. Section 4 describes the ensemble seasonal prediction experiments. In Sect. 5, seasonal prediction results using the new algorithm are presented and compared with those using the other ensemble generation methods. A summary of the study and discussions are included in Sect. 6.

2 Model and data assimilation experiments In this study, a hybrid coupled model is used to investigate the ensemble perturbations for ENSO predictions. The oceanic component of the hybrid coupled model is based on an intermediate ocean model similar to the Cane-Zebiak (CZ) model (Zebiak and Cane 1987). It differs from the original in that it uses a new method for parameterization of the

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

subsurface temperature (Kang and Kug 2000). The horizontal resolution of the modified CZ model is 5.625° (2°) in the longitudinal (latitudinal) direction. The atmospheric component of the hybrid model is the SPEEDY (Simplified Parameterizations, primitivE-Equation DYnamics, Molteni 2003) AGCM (Atmospheric Global Circulation Model). According to Molteni (2003), the SPEEDY model simulates the general structure of the global atmospheric circulation fairly well, and some aspects of the systematic errors are similar to many AGCMs. The resolution of the model is T42L10 (horizontal spectral truncation of 42 wavenumbers, or approximately 2.8°, and 10 vertical levels). The air–sea coupling interval is 10 days. The oceanic (atmospheric) model provides its anomalous SST (anomalous wind stress), and receives anomalous zonal and meridional wind stresses (anomalous SST), whose value is 10-day averaged. Note that the oceanic part of the hybrid coupled model calculates only an anomaly. Details on the coupled model and its performance can be found in Ham et al. (2009), or An et al. (2010). The CZ-SPEEDY coupled model has been used for data assimilation to generate initial conditions for seasonal forecasts. The experiments are performed in a perfect model context, i.e., one realization of the coupled model is regarded as the true state. To obtain the analysis values (initial conditions for prediction), the EnKF is configured with 16 ensemble members. The analysis interval is 10 days, and both SST and thermocline depth anomalies are assimilated. In this study, both SST and thermocline depth observations are extracted from the coupled model simulation every 11.25° (4°) in the longitudinal (latitudinal) direction. To mimic the errors in real observations, the prescribed observational errors in the SST and thermocline depth anomaly are 0.35°C and 4 m, respectively, and commensurate random perturbations are added to the model output to mimic the measurement and representative errors in observations. This perfect model approach is often used to evaluate the quality of data assimilation system under the perfect model assumption (Houtekamer and Mitchell 1998; Zhang et al. 2006). No covariance localization or inflation is applied. The assimilation period is 25 years, from year 15 to 40 of the coupled model simulation.

3 Flow-dependent empirical singular vector (FESV) within an EnKF framework 3.1 Flow-dependent empirical singular vector (FESV) In this section, details of the method to extract flowdependent empirical singular vectors (FESV) within the EnKF framework are provided. Assume that the data

assimilation interval is T, and assimilation is performed n times from time 0. Then ensemble perturbations which are freely integrated from time 0 to T, T to 2T, 2T to 3T, and so on, until the last analysis time are available. By using these ensemble members, an empirical linear operator can be obtained for each analysis time (Blumenthal 1991; Moore and Kleeman 2001; Penland 1996; Kug et al. 2010; Ham and Kang 2010). The procedure to calculate the empirical linear operator is as follows. First, assume that the nonlinear integration can be approximated by a linear operator (L) for the evolution of the state vector from time 0 to time T as follows: 0

0

YT ¼ LðT; 0ÞX0 þ 2;

ð1Þ

0

Where X0 , Y0 T, and L(T, 0) are, respectively, ensemble perturbations at time 0, T, and an empirical linear operator 0 which approximately transforms the initial state X0 into a 0 forecast state Y0 T. The dimensions of X0 , and Y0 T are m 9 k where m (k) is the size of the state vector (number 0 of ensemble members). Note that X0 is the ensemble perturbation (the deviation of the state vector from the  0 ) after the EnKF analysis at ensemble mean X00 ¼ X0  X 0 time 0, and Y T is the ensemble perturbation after the short term forecast from time 0 to T from the initial condition X0. The variable 2 denotes the error in the linear approximation. It is possible to estimate the operator empirically from ensemble samples, i.e., the linear operator can be calculated as: 0

0

0

0

LðT; 0Þ ¼ YT X0T ðX0 X0T Þ1 :

ð2Þ

Note that the dimension of L is m 9 m. After calculating the empirical linear operator for each analysis time, a linear operator for a longer time-scale is calculated using the hypothesis of causality as follows (Terwisscha and Dijkstra 2005): LðnT; 0Þ  LF ðnT; 0Þ ¼ LðnT; ðn  1ÞTÞ. . .Lð2T; TÞLðT; 0Þ;

ð3Þ

where LF (nT, 0) is the final linear operator which can propagate ensemble perturbations from time 0 to nT. Note that because of the analysis procedure between time 0 and 0 Tn, LF (nT, 0) cannot be calculated directly from X0 and Y0 nT, i.e., the relationship of ensemble perturbations (ensemble covariance) between time 0 and nT with assimilation is different from the relationship during a freerunning integration without assimilation. Also, the error term in (1) has been dropped because for ENSO forecasts the linear assumption is quite valid for a two or three season forecast (e.g. Xu and von Storch 1990; Tang 1995). Using the longer time-scale linear operator in Eq. 3, singular vectors for longer time-scales can be extracted by

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

applying a singular value decomposition (SVD) of the linear operator LF (nT, 0): ui YnT ¼ si vi X0 ;

ð4Þ

where si, ui, and vi are the ith singular value of LF (nT, 0) and its right and left singular vectors, respectively. If the singular value si, is greater than one, this implies that the singular mode grows in the linear operator. 3.2 Application of FESV to the CZ-SPEEDY model In order to apply the FESV to prediction with the CZSPEEDY model, the linear operator is obtained for a reduced space defined through EOF analysis. For both initial and final state vectors X and Y in Eq. 1, the EOF analysis is applied to instantaneous thermocline depth data using a 100-year free running simulation. Because the memory for ENSO prediction resides primarily in the thermocline depth, only this variable is adopted to represent perturbations of FESV. Then, magnitudes of 16 ensemble members, which is the number of ensemble members used in the EnKF assimilation, projected onto the first seven EOFs are used to construct the X and Y matrices. Note that the seven dominant EOFs explain about 85% of the total themocline depth variability, and the minor EOF modes, which are hard to explain physically, are excluded. Therefore, the dimension of the X and Y matrices is 7 9 16 in this study. Note that the major results of this study are not sensitive to the number of the EOF modes. After generating FESV under EOF space, the spatial patterns of the FESV are reconstructed by multiplying the FESV amplitude by the corresponding EOF eigen-vectors. The optimal time interval for the FESV is 3-months using the 3-month assimilation data prior to forecast initialization. For example, to obtain the empirical operator for September 1 of any year, assimilation data from June 1 to September 1 of that year are used. Note that this approach clearly does not use any observations during forecast periods.

4 The seasonal forecast experiments For the seasonal forecast experiments, the first unstable mode from the FESV is added and subtracted to ensemble mean value of initial conditions generated by the EnKF algorithm. For comparison, seasonal predictions were also conducted from initial conditions generated by adding and subtracting the first unstable mode of the conventional ESV (CESV) which is static in time, with an optimal time (i.e. time lags between initial and final variables) of 3 months, to the same ensemble mean initial condition used for the

123

seasonal prediction experiments with the FESV. Note the 3-month optimal time is selected to make a fair comparison with the FESV. The CESV for each calendar month was calculated using the single 100-year long simulation of Kug et al. (2010). Therefore, CESVs only have seasonal cycles. In this study, we use eight ensemble members (i.e. four sets of CESV and FESV perturbations). To generate additional sets of ESV perturbations, the difference between randomly selected ensemble member among 16 ensemble members and ensemble mean is added to ESV perturbations. Note that the difference is added for all oceanic and atmospheric initial condition to reduce the inconsistency between variables and minimize the initial shock. With these initial conditions using 8 ensemble members, predictions of 12-month duration are performed over a 20-year period from year 21 to 40, with initialization on the first day of every month. Thus, there are total of 240 prediction experiments in all. Hereafter, forecasts made using the conventional ESV are denoted as ‘‘CESV prediction’’, and those made with the Flow-dependent ESV are denoted as ‘‘FESV prediction’’.

5 Results 5.1 Validation of the flow-dependent linear operator Prior to validating the quality of the flow-dependent linear operator, it is worthwhile to mention the importance of the time-varying component of the ensemble perturbation. To measure the magnitude of the ensemble perturbation at each grid point, the ensemble spread at each grid point is defined as follows. sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rens 1 nX  i;j Þ2 ; Ensemble Spreadi;j ¼ ðXk  X ð5Þ nrens k¼1 i;j  i;j are the number of ensemble Where, nrens, and Xki;j  X members, and the deviation of the state vector of kth ensemble member from the ensemble mean at ith (jth) grid point in zonal (meridional) direction, respectively. Then, the magnitude of climatological (time-mean) ensemble spread is compared to the SD of ensemble spread to measure the relative magnitude of the time-mean perturbation magnitude to that of time-varying perturbation magnitude. Figure 1 shows the climatological ensemble spread and SD of the ensemble spread of the thermocline using the CZ-SPEEDY model. The climatological ensemble spread is largest over the equatorial regions. Similarly, the SD of the ensemble spread is largest over the equatorial regions, about 40–50% of the climatological spread. This means that the

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

perturbation growth is time-dependent, justifying the need for a flow-dependent linear operator. To investigate how well the flow-dependent linear operator represents the nonlinear operator, a statistical prediction test in conducted using the flow-dependent operator. To compare the quality of the statistical prediction, statistical prediction using a static (time-fixed) linear operator is also performed. For example, statistical prediction using a flow-dependent linear operator for September 1, year 21 can be obtained YFESV ¼ LF ðJun 1 21; Mar 1 21ÞXJune 1 21 . Note that the flow-dependent operator is constructed using results from a previous data assimilation cycle, so that no observations during the forecast are used. Similarly, statistical prediction for September 1 year 21 is conducted using the static linear operator from June to September. Figure 2 shows the skill for the thermcline depth anomaly measured as the root-mean-square error (RMSE) between the statistical predictions and the nonlinear integration. The total number of cases is 240 (12 months 9 20 years) from year 21 to 40. In statistical prediction using the static operator, the RMSE is relatively large over the eastern Pacific, and off-equatorial Pacific. The relative large RMSE over the eastern Pacific is due to the large

variability over the eastern Pacific related to the ENSO. The large RMSE over the off-equatorial Pacific is due to the fact that the static operator, which is fixed in time, fails to predict the evolution of propagating oceanic Rossby waves over the off-equatorial regions. On the other hand, RMSE using the flow-dependent operator is systematically reduced over the off-equatorial regions. Thus, the flowdependent linear operator better represents the time-varying oceanic wave propagation over the off-equatorial regions. In turn, the propagation of the oceanic Rossby waves reflected from off-equatorial regions to equator affects the SST forecast skill over the eastern Pacific (Kang and An 1998). 5.2 The comparison of the FESV with CESV Figure 3 shows the singular value of FESV and CESV starting from April which is the month when the time-

(a)

(a) (b)

(b)

Fig. 2 The root-mean-square errors (RMSE) for the thermcline depth from the nonlinear integration and the statistical predictions using a the static operator, and b the flow-dependent operator

(c)

Fig. 1 Spatial pattern of the 20-year averaged a climatological ensemble spread of thermocline depth, b SD of the ensemble spread, and c ratio of the SD of the ensemble spread to the climatological ensemble spread multiplied by 100. Sixteen ensemble members are used, and the forecast lead time is 10 days. The unit is m

Fig. 3 The singular value of FESV (red line), and CESV (black line) starting at April 1st

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

averaged singular value of the CESV is the largest. It is obvious that there is systematic year-to-year variability of singular value in FESV. For example, the singular value at year 24 is smallest as 1.6, and that at year 36 is 4.08. It implies that the unstable mode captured as singular vector is different in time, therefore, it emphasize the positive impact of time-dependent perturbations like FESV. On the other hand, the singular value of CESV is fixed in time, because it captures state-independent singular vectors. The time-averaged singular value for FESV is 2.86, and the singular value for CESV is 2.79, which is similar to each other. Figure 4 shows the first EOF of FESV and CESV starting from April. For comparison, dominant EOF of the Bred Vector (BV) starting from April is also shown. For the breeding experiment, the norm is SST averaged over the equatorial central-eastern Pacific (180–80oW, 10oS– 10oN). The breeding cycle is 3-monthly, which is the same as the optimal time of the FESV and CESV. Breeding is performed from year 15 to 40. The control run is started from the ensemble mean initial condition from the EnKF assimilation, and the initial perturbation for the perturbed run is defined as the difference of the analyses on January 15 and February 15. Note that each breeding cycle is

started from ensemble mean initial condition from the EnKF assimilation to obtain the BV under the initialized states. In addition, to investigate the sensitivity of optimal time to extract the ESV, ESV are also calculated using the linear operator from the last analysis cycle whose time interval is 10 days, [i.e., L (nT, (n - 1)) in Eq. 3]. This is done to extract the fast-growing direction from all the directions spanned by the ensemble perturbations of the EnKF, therefore, it will be denoted as ‘EnKF perturbation’ in Fig. 4c, g. The initial state of the thermocline of the CESV (Fig. 4a) shows a deepened thermocline over the eastern and western Pacific, and shoaling in the off-equatorial region. Even though it is somewhat reminiscent of the recharge pattern of the equatorial heat content in the Recharge Oscillator (e.g. Jin 1996, 1997a, b), it also reflects the thermocline pattern of the mature phase of ENSO, with a maximum over the eastern Pacific. In contrast, the spatial pattern of the first EOF of the FESV shows a clear ENSO recharge pattern. The initial state of the thermocline of the FESV shows a deepened thermocline over the entire equatorial Pacific with a peak value near 150oW. This spatial pattern is more consistent with that of the observed/theoretical thermocline

(a)

(e)

(b)

(f)

(c)

(g)

(d)

(h)

Fig. 4 The first EOF of FESV, CESV, EnKF perturbations, and BV starting from April. Left panels show the initial thermocline depth anomaly of a CESV, b FESV, c EnKF perturbations, and d BV with a

123

3-month breeding interval. Right-hand panels show the final SST anomalies of e CESV, f FESV, g EnKF perturbations, and h BV

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

anomaly pattern during the developing phase of ENSO when the anomalous SST growth is maximum (Jin 1997a, b; Zelle et al. 2004). Different again, the dominant unstable mode from the EnKF perturbations shows an east–west contrast pattern with positive (negative) peak over the far eastern (western) Pacific, which is similar to the pattern of the mature phase of ENSO. Using a simple two-strip version of a Cane-Zebiaktype coupled model, Jin (1997a, b) identified two essential components for the instability related to ENSO. He emphasized the role of equatorial waves in the establishment of the quasi-equilibrium equatorial Sverdrup balance and the slow recharge/discharge of equatorial heat content. The first process is related to Bjerknes feedback. When a westerly wind anomaly is induced, it generates anomalous thermocline deepening (shoaling) over the eastern (western) Pacific with the related positive SST anomalies over the eastern Pacific. Because the westerly wind is maximum over the central Pacific, the east–west contrast in thermocline depth anomaly results. The positive SST anomaly over the eastern Pacific then reinforces the westerly wind anomaly due to the weakening of the Walker Circulation. Thus, positive air–sea feedback reinforces the anomalous westerly wind, the east– west contrast in the thermocline depth anomaly, and the resultant SST anomaly over the eastern Pacific. This process is relatively fast compared to the time-scale of ENSO. On the other hand, the second process is relatively slow and has a comparable time-scale to ENSO. The shallow (deep) zonal mean thermocline depth is controlled by weakening (strengthening) of the trade wind system, and is maintained during the whole El Nino (La Nina) event. Based on the concepts above, the differences between FESV (or CESV), and EnKF perturbations may be caused by capturing different ENSO-related instabilities with different time scales for optimal growth. It is likely that both FESV and CESV, whose optimal time is 3 months, captures the instability related to the slow recharge/discharge process. This process is characterized by a zonally-elongated anomaly which is similar to the FESV. On the other hand, the EnKF perturbation, whose assimilation interval is 10 days, captures the instability related to the fast Bjerknes feedback. This feedback is characterized by the thermocline depth with an east–west contrast pattern, which is similar to the dominant mode of EnKF perturbations. It is also interesting to look at the dominant pattern that emerges from breeding. The BV pattern shows a clear east– west contrast pattern even though the breeding cycle is 3 months. This is the results of the breeding methodology in which the BV pattern would be similar to the pattern during the mature phase of the target phenomenon because the BV is the perturbation ‘‘grown’’ over the breeding period. Note that, the BV is used as the initial perturbation of the next breeding cycle after the rescaling process.

Therefore, the BV pattern is likely to be similar to the pattern during the mature phase of ENSO. The spatial pattern for the final SST anomaly is calculated by linear regression, although only thermocline depth data for the final state are used in the calculation of the FESV and CESV. From these initial states, unstable modes show large SST anomalies over the eastern Pacific. This feature is similar to the singular mode of Xue et al. (1997a, b). Note that the final SST state from the FESV includes a negative SST anomaly over the off-equatorial western Pacific, while final SST state from other methods does not. Thus, the dominant mode of the FESV leads to growth in SST anomaly over the western as well as the eastern Pacific. The negative SST anomaly over the western Pacific is well observed and simulated in CZ-type models during the mature phase of ENSO (Cane 1983; Zebiak and Cane 1987). These empirical singular vectors are also compared to the classical SV analysis based on tangent-linear model (TLM) (Chen et al. 1997; Moore et al. 2003; Tang et al. 2006; Hawkins and Sutton 2009; Hawkins et al. 2011). To calculate the classical SV using a tangent linear model (TLM), we selected one target year when the FESV is most similar to the dominant EOF of FESV. Then, the perturbation is added to the ensemble mean data assimilation analysis to formulate the TLM. In this study, we introduced 15 dominant EOF modes of thermocline depth as perturbations. The perturbations for other variables are generated as linearly regressed field onto 15 dominant EOF PCs. Note that the 15 dominant EOFs explain 92% of total thermocline depth variability, therefore, 15 perturbations based on EOF modes would extract the SV effectively with reduced computational cost (Moore et al. 2003; Tang et al. 2006; Hawkins and Sutton 2009; Hawkins et al. 2011). Figure 5 shows the initial SV of thermocline depth using the TLM, the conventional empirical singular vector method, and the flowdependent empirical singular vector method at the target year. Note that the final SST patterns using all methods are El Nino-like, as already shown in Fig. 4. The initial classical SV shows the maximum positive values over the central Pacific between 180 and 150oW, which is consistent results with Xue et al. (1997a). The positive maximum of the initial FESV is also located over the central Pacific similar to the initial classical SV, while the maximum of initial CESV is located over the eastern Pacific. Interestingly, the offequatorial spatial pattern of the initial FESV is similar to that of initial classical SV, but in the opposite hemisphere. However, the FESV captures the dominant features of the classical time-dependent singular vector. 5.3 Impacts of FESV on seasonal forecasts The impact on seasonal forecasts of the optimal initial perturbations from FESV and CESV is shown in Fig. 6.

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ðN  1Þð1 þ r23 Þ t ¼ ðr12  r13 Þu  u  t N1 ðr12 þ r13 Þ2 ð1  r23 Þ2 2 jRj þ N3 4

(a)

(b)

(c)

Fig. 5 The initial singular vector of thermocline depth using a Tangent Linear Model (TLM), b conventional empirical singular vector (CESV) method, and c flow-dependent empirical singular vector (FESV) method at the target year. Note that the target year is selected when the FESV is most similar to the dominant FESV mode

The FESV predictions have a better skill than the CESV predictions over all forecast lead times, particularly for forecasts at 6-month or longer leads. The correlation improvement increases with lead time, and the difference saturates after about 7 months. These results support the view that the optimal initial perturbation can improve the prediction skill on seasonal time scales. The RMSE statistics, as shown in Fig. 6b, show similar results. To confirm that the improvement of correlation skill using the FESV is significant, we performed a Hotelling’s t test in Meng et al. (1992). The t test takes the form

Fig. 6 a Correlation, and b RMS errors for NINO3 forecast SST anomalies based on the FESV prediction (red line), CESV (black line). Note that the dots on the anomaly correlation denote that the correlation improvement in FESV prediction is significant with 90% confidence level

123

(a)

ð6Þ   where jRj ¼ 1  r212  r213  r223 þ ð2r12 r13 r23 Þ; and r12, r13, r23 is the correlation between ESV and true, ESV and true, and ESV and CNTL, respectively. N is 240 in this study. Note that there is possibility that the significance is exaggerated because there is autocorrelation between the samples to reduce the effective number of degree of freedom. The significant test shows that the correlation improvement using FESV is significant with 90% confidence level after 6-month lead (i.e. the dot on the correlation denotes the correlation improvement in FESV prediction is above the 90% confidence level). Because the fast-growing perturbation is extracted, the FESV should grow very rapidly as the prediction starts. Therefore, it is expected that the ensemble spread of the FESV will be larger than other initial perturbations. In order to check the growth of the spread using the ensemble spread of the NINO3 index (SST anomaly averaged over 150–90oW, 5oS–5oN). By calculating the ratio of the ensemble spread of the initial condition and the monthly mean forecast at each forecast lead time, we can estimate the growth rate of the SST perturbation. Figure 7 shows this estimated growth rate for the CESV and FESV predictions. In the following, comparisons are made only between FESV and CESV to focus on the positive impact of the flowdependent operator. The small initial spread grows rapidly as the prediction starts in both cases. However, the spread of the FESV is larger than that of the CESV from month 3 and is largest in 6–10-month lead forecasts. In particular, for forecasts with a 9-month lead, the ensemble spread of the FESV is about twice that of the CESV. This supports the view that the FESV grows faster than the CESV.

(b)

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

Fig. 7 The ratio of ensemble spread of monthly-mean NINO3 (150–90oW, 5oS–5oN) SST at each forecast month to that of initial conditions in the CESV (black) and FESV (red) prediction over the duration of the forecast. The total number of cases is 240

Figure 8 shows the spatial patterns of SST anomaly correlation skill from the CESV prediction (contour). It also shows the difference in skill between the FESV and CESV predictions with 90% confidence level (shading). A positive sign denotes that the forecast skill of the FESV prediction is higher than that of the CESV prediction. The correlation skill of the CESV prediction is relatively low in the western equatorial Pacific compared to that in the

off-equatorial eastern Pacific at 6-month lead forecasts. This deficiency over the equatorial eastern Pacific is rectified somewhat in FESV predictions. For example, the improvement in the FESV prediction is concentrated over the equatorial regions at 6- and 9-month lead forecasts. It is also noteworthy that the correlation skill is systematically improved over regions where the skill is relatively low. For example, at 6- and 9-month forecast leads, improvements in forecast skill are significant over the far eastern Pacific, where the forecast skill of the control prediction is relatively low. According to Lorenz, the forecast skill of the atmospheric model depends not only on the accuracy of the initial conditions but also on the instability of the flow itself in the perfect model context (Lorenz 1963a, b, 1965). The forecast skill is degraded when the initial errors grow fast even if the initial errors are relatively small. Hence, the effectiveness of the FESV is related to the predictability, and the positive impact of the FESV optimal ensemble perturbations is emphasized over the regions where there is strongly unstable flow.

6 Summary and discussions In this study, we present a new approach for extracting flowdependent empirical singular vectors (FESVs) for seasonal prediction using ensemble perturbations obtained from an

(a)

(b)

(c)

(d)

Fig. 8 The spatial pattern of correlation for forecast SST anomalies in the CESV prediction (contour), and forecast skill difference between the FESV and CESV predictions (shading). Note that a

positive difference (red shading) denotes that the FESV forecast skill improves upon that of the CESV with 90% confidence level. The shaded contour interval of shading is 0.02

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

ensemble Kalman filter (EnKF) assimilation. Due to the short interval between analyses, EnKF perturbations contain the instability related to the fast time scale variability. To isolate slower coupled instabilities which would be more relevant to seasonal prediction, an empirical linear operator for seasonal time-scales (e.g. several months) is formulated using a hypothesis of causality; then ESVs from a linear operator appropriate for seasonal time-scales is extracted. In this study, nine sets of 10-day interval ensemble perturbations are gathered to formulate the linear operator with a 3-month optimal time. Statistical prediction tests using the linear operators show that the flow-dependent linear operator constructed using a causality assumption predicts the nonlinear integration results better than a static operator. Because the ensemble perturbation extracted using the flow-dependent linear operator represents the actual fastgrowing perturbations well, it grows faster than the other initial perturbations tested in this study. For SST over the central-eastern Pacific, the FESV grows about twice as fast as that of the CESV at a 9-month lead forecast. Through a 20-year series of 12-month seasonal predictions initialized every month, it is shown that the average forecast skill using FESV significantly improves upon that using CESV. For example, the correlation skill of the NINO3 index using FESV is higher by about 0.1 than that from CESV at 8-month leads with 90% confidence level. In addition, the forecast skill using the FESV is superior to that using the other methods at most lead times and regions. In particular, the improvement of the forecast skill is significant where the correlation skill of CESV is relatively low. One can be curious that 16 ensemble members to derive FESV are enough to extract stable optimal perturbations, because the size of 16 samples for statistical analysis seems to be small. Therefore, a sensitivity experiment of FESV modes to ensemble size is checked. Figure 9 shows the dominant EOF of thermocline depth FESV with 16 and 32 ensemble members. It is found that that the spatial pattern of dominant EOF of FESV is not sensitive to the number of ensemble members. As well as the spatial pattern of the dominant EOF, the time-averaged pattern correlation between FESV with 16 ensemble members and that with 32 ensemble members for each year is 0.58, even though the pattern correlation is relative small in some years (ex. year 27). This insensitivity of FESV to ensemble size implies that the 16 ensemble members are somewhat enough to span possible directions of error growth. This is probably related to the fact that the degrees of freedom in the CZ-SPEEDY model are smaller than that in a complex coupled model. The present method has several benefits compared to the breeding method, which is a well-known method for detecting growing perturbations under time-varying flows. Firstly, no additional computational costs are required to obtain the FESV. For example, FESV are easily obtained

123

using ensemble members from the EnKF which have already been generated for the data assimilation system. In contrast, additional runs are needed to obtain a BV. Secondly, the FESV may be used to generate several ensemble perturbations. For example, when n EOF modes are used to formulate the linear operator, n ensemble perturbations are extracted which are orthogonal each other. This means that ensemble predictions using FESVs would be more effective than the single perturbation from breeding because each ensemble perturbation represents an independent direction of possible forecast errors. In contrast, the BV is designed to capture a single fastest growing perturbation among many unstable modes in nature. Even though there have been several attempts to detect independent bred vectors using masked breeding, or different norms, these have not yet proven successful for seasonal prediction. For example, Yang et al. (2009) tested the sensitivity to the norm over the tropics for coupled BVs using a CGCM, but they concluded that equatorial BVs generally have similar dominant structures, indicating that different norms are not effective in capturing multiple unstable modes for seasonal prediction. The optimal ensemble perturbations for seasonal prediction should contain the large-scale dynamics associated with the ENSO. For example, Yang et al. (2006, 2008) show the dominant BV patterns are similar to the anomalies during the mature phase of ENSO. Such patterns are not extracted with a short-time (less than a week) free integration interval because the accompanying growth

(a)

(b)

(c)

Fig. 9 The dominant EOF of thermocline depth FESV with a 16 ensemble members, and b 32 ensemble members starting at April 1st. c The pattern correlation between FESV with 16 ensemble members and that with 32 ensemble members for each year

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector

related to the coupled instability is no larger than that of weather instability. Therefore, it is unlikely that an EnKF perturbation whose analysis interval is less than a few days would contain the proper magnitude of coupled instability. On the other hand, the FESV which are extracted from a flow-dependent linear operator with a time scale of several months can represent the time-varying coupled instability related to ENSO. Therefore, if the method presented here translates to a CGCM, it offers the possibility of improving current state-of-art seasonal prediction. This will be the subject of a future investigation. Acknowledgements We appreciate the helpful suggestions and comments from two anonymous reviewers. This research was supported by the NASA Modeling, Analysis and Prediction program under WBS 802678.02.17.01.25. Computational resources for this study were provided by the NASA Center for Climate Simulation.

References An S-I, Ham Y-G, Kug J-S, Timmermann A, Choi J, Kang I-S (2010) The inverse effect of annual mean state and annual cycle changes on ENSO. J Clim 23:1095–1110 Blumenthal MB (1991) Predictability of a coupled ocean-atmosphere model. J Clim 4:766–784 Buehner M (2005) Ensemble-derived stationary and flow-dependent background error covariances: evaluation in a quasi-operational NWP setting. Q J R Meteorol Soc 131:1013–1044 Buizza R, Petroliagis T, Palmer TN, Hamrud M, Hollingsworth A, Simmons A, Wedi N (1998) Impact of model resolution and ensemble size on the performance of an ensemble prediction system. Q J R Meteorol Soc 124:1935–1960 Buizza R, Houtekamer PL, Toth Z, Pellerin G, Wei M, Zhu Y (2005) A comparison of the ECMWF, MSC, and NCEP global ensemble prediction systems. Mon Weather Rev 133:1076–1097 Cai M, Kalnay E, Toth Z (2003) Bred vectors of the Zebiak-Cane model and their potential application to ENSO predictions. J Clim 16:40–56 Cane M (1983) Oceanographic events during El Nino. Science 222:1189–1195 Chen Y-Q, Battisti DS, Palmer TN, Barsugli J, Sarachik ES (1997) A study of the predictability of tropical Pacific SST in a coupled atmosphere-ocean model using singular vector analysis: the role of the annual cycle and the ENSO cycle. Mon Weather Rev 125:831–845 Corrazza M, Kalnay E, Patil DJ, Yang S-C, Morss R, Cai M, Szunyogh I, Hunt BR, Yorke JA (2003) Use of the breeding technique to estimate the structure of the analysis ‘‘errors of the day’’. Nonlinear Process Geophys 10:233–243 Evensen G (1994) Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J Geophys Res 99:10 143–10 162 Evensen G (2003) The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dyn 53:343–367 Evensen G, van Leeuwen PJ (2000) An ensemble Kalman smoother for nonlinear dynamics. Mon Weather Rev 128:1852–1867 Fan Y, Allen MR, Anderson DLT, Balmaseda MA (2000) How predictability depends on the nature of uncertainty in initial conditions in a coupled model of ENSO. J Clim 13:3298–3313 Farrell BF (1989) Optimal excitation of baroclinic waves. J Atmos Sci 46:1193–1206

Ham Y-G, Kang I-S (2010) Growing-error correction of ensemble Kalman filter using empirical singular vectors. Q J R Meteorol Soc. doi:10.1002/qj.711 Ham Y-G, Kug J-S, Kang I-S (2009) Optimal initial perturbations for El Nino ensemble prediction with ensemble Kalman filter. Clim Dyn 33:959–973. doi:10.1007/s00382-009-0582-z Hamill TM, Snyder C, Whitaker JS (2003) Ensemble forecasts and the properties of flow-dependent analysis-error covariance singular vectors. Mon Weather Rev 131:1741–1758 Hawkins E, Sutton R (2009) Decadal predictability of the Atlantic Ocean in a coupled GCM: forecast skill and optimal perturbations using linear inverse modeling. J Clim 22:3960–3978 Hawkins E, Robson J, Sutton R, Smith D, Keenlyside N (2011) Evaluating the potential for statistical decadal predictions of surface temperatures with a perfect model approach. Clim Dyn. doi:10.1007/s00382-011-1023-3 Houtekamer PL, Mitchell HL (1998) Data assimilation using an ensemble Kalman filter technique. Mon Weather Rev 126:796– 811 Jin F-F (1996) Tropical ocean-atmosphere interaction, the Pacific cold tongue, and the El Nin˜o-Southern Oscillation. Science 274:76–78 Jin F-F (1997a) An equatorial ocean recharge paradigm for ENSO. Part I: conceptual model. J Atmos Sci 54:811–829 Jin F-F (1997b) An equatorial ocean recharge paradigm for ENSO. Part II: a stripped-down coupled model. J Atmos Sci 54:830–847 Kang I-S, An S-I (1998) Kelvin and Rossby wave contributions to the SST oscillation of ENSO. J Clim 11:2461–2469 Kang I-S, Kug J-S (2000) El-Nino prediction system using an intermediate ocean and a statistical atmosphere. Geophys Res Lett 27:1167–1170 Keppenne CL, Rienecker MM (2002) Initial testing of a massively parallel ensemble Kalman filter with the Poseidon isopycnal ocean general circulation model. Mon Weather Rev 130:2951–2965 Keppenne CL, Rienecker MM, Jacob JP, Kovach R (2008) Error covariance modeling in the GMAO ocean ensemble Kalman filter. Mon Weather Rev 136:2964–2982 Kleeman R, Tang Y, Moore AM (2003) The calculation of climatically relevant singular vectors in the presence of weather noise as applied to the ENSO problem. J Atmos Sci 60:2856– 2868 Kug J-S, Ham Y-G, Kimoto M, Jin F-F, Kang I-S (2010) New approach for optimal perturbation method in ensemble climate prediction with empirical singular vector. Clim Dyn 35:331–340. doi:10.1007/s00382-009-0664-y Leeuwenburgh O (2007) Validation of an EnKF system for OGCM initialization assimilating temperature, salinity, and surface height measurements. Mon Weather Rev 135:125–139 Leith CE (1974) Theoretical skill of Monte Carlo forecasts. Mon Weather Rev 102:409–418 Leutbecher M, Palmer TN (2008) Ensemble forecasting. J Comput Phys 227:3515–3539 Lorenz EN (1963a) Deterministic non-periodic flow. J Atmos Sci 20:130–141 Lorenz EN (1963b) The predictability of hydrodynamic flow. Trans NY Acad Sci Ser II 25:409–432 Lorenz EN (1965) A study of the predictability of a 28-variable atmospheric model. Tellus 17:321–333 Meng XL, Rosenthal R, Rubin DB (1992) Comparing correlated correction coefficients. Psychol Bull 111:172–175 Molteni F (2003) Atmospheric simulations using a GCM with simplified physical parameterizations. I: model climatology and variability in multi-decadal experiments. Clim Dyn 20:175–191 Molteni F, Palmer TN (1993) Predictability and finite time instability of the northern winter circulation. Q J R Meteorol Soc 119:269–298

123

Y.-G. Ham, M. M. Rienecker: Flow-dependent empirical singular vector Molteni F, Buizza R, Palmer TN, Petroliagis T (1996) The ECMWF ensemble prediction system: methodology and validation. Q J R Meteorol Soc 122:73–119 Moore AM, Kleeman R (1996) The dynamics of error growth and predictability in a coupled model of ENSO. Q J R Meteorol Soc 122:1405–1446 Moore AM, Kleeman R (1997a) The singular vectors of a coupled ocean-atmosphere model of ENSO I: thermodynamics, energetic and error growth. Q J R Meteorol Soc 123:953–981 Moore AM, Kleeman R (1997b) The singular vectors of a coupled ocean-atmosphere model of ENSO II: thermodynamics, energetic and error growth. Q J R Meteorol Soc 123:953–981 Moore AM, Kleeman R (1998) Skill assessment for ENSO using ensemble prediction. Q J R Meteorol Soc 124:557–584 Moore AM, Kleeman R (2001) The differences between the optimal perturbations of coupled models of ENSO. J Clim 14:138–163 Moore AM, Tang Y, Moore A (2003) The calculation of climatically relevant singular vectors in the presence of weather noise. J Atmos Sci 60:2856–2867 Palmer TN, Buizza R, Molteni E, Chen Y-Q, Corti S (1994) Singular vectors and the predictability of weather and climate. Philos Trans R Soc Lond 348:459–475 Pen˜a M, Kalnay E (2004) Separating fast and slow modes in coupled chaotic systems. Nonlinear Process Geophys 11:319–327 Penland C (1996) A stochastic model of IndoPacific sea surface temperature anomalies. Physica D 96:534–558 Tang B (1995) Periods of linear development of the ENSO cycle and POP Forecast experiments. J Clim 8:682–691 Tang Y, Kleeman R, Miller S (2006) ENSO predictability of a fully coupled GCM model using singular vector analysis. J Clim 19:3361–3377 Terwisscha van Scheltinga A, Dijkstra H (2005) Nonlinear data assimilation using implicit models. Nonlinear Process Geophys 12:515–525 Toth Z, Kalnay E (1993) Ensemble forecasting and NMC: the generation of perturbations. Bull Am Meteorol Soc 74:2317–2330 Toth Z, Kalnay E (1997) Ensemble forecasting at NCEP and the breeding method. Mon Weather Rev 125:3297–3318 Tziperman E, Zanna L, Penland C (2008) Nonnormal thermohaline circulation dynamics in a coupled ocean–atmosphere GCM. J Phys Oceanogr 38:588–604

123

Xu J-S, von Storch H (1990) Predicting the state of the Southern Oscillation using principal oscillation pattern analysis. J Clim 3:1316–1329 Xue Y, Cane MA, Zebiak SE, Blumenthal B (1994) On the prediction of ENSO: a study with a low-order Markov model. Tellus 46A:512–528 Xue Y, Cane MA, Zebiak SE (1997a) Predictability of a coupled model of ENSO using singular vector analysis. Part I: optimal growth in seasonal background and ENSO cycles. Mon Weather Rev 125:2043–2056 Xue Y, Cane MA, Zebiak SE (1997b) Predictability of a coupled model of ENSO using singular vector analysis. PartII: optimal growth and forecast skill. Mon Weather Rev 125:2057–2073 Yang S-C, Cai M, Kalnay E, Rienecker M, Yuan G, Toth Z (2006) ENSO bred vectors in coupled ocean-atmosphere general circulation models. J Clim 19:1422–1436 Yang SC, Kalnay E, Cai M, Rienecker MM (2008) Bred vectors and tropical pacific forecast errors in the NASA coupled general circulation model. Mon Weather Rev 136:1305–1326 Yang SC, Keppenne C, Rienecker M, Kalnay E (2009) Application of coupled bred vectors to seasonal-to-interannual forecasting and ocean data assimilation. J Clim 22:2850–2870 Zebiak SE, Cane MA (1987) A model for El Nino—southern oscillation. Mon Weather Rev 115:2262–2278 Zelle H, Appeldoorn G, Burgers G, van Oldenborgh GJ (2004) On the relationship between sea surface temperature and thermocline depth in the eastern equatorial Pacific. J Phys Oceanogr 34:643–655 Zhang S, Harrision MJ, Wittenberg AT, Rosati A (2005) Initialization of an ENSO forecast system using a parallelized ensemble filter. Mon Weather Rev 133:3176–3201 Zhang F, Meng Z, Aksoy A (2006) Tests of an ensemble Kalman filter for Mesoscale and regional-scale data assimilation. Part I: perfect model experiments. Mon Weather Rev 134:722–736. doi: 10.1175/MWR3101.1 Zhu Y, Toth Z, Wobus R, Richardson D, Mylne K (2002) The economic value of ensemble-based weather forecasts. Bull Am Meteorol Soc 83:73–83