A global perspective of the limits of prediction skill of NWP models

Abstract This paper discusses the scale-dependent growth of the global forecast uncertainties simulated by the operational ensemble prediction system of the European Centre for Medium-Range Weather Forecasts. It is shown that the initial uncertainties are largest in the tropics and have biggest amplitudes at the large scales. The growth of forecast uncertainties (ensemble spread) takes place at all scales from the beginning of forecasts. The growth is nearly uniform in the zonal wavenumbers 1–5 and strongly scale-dependent in the larger wavenumbers. Moreover, the growth from initial uncertainties at large scales appears dominant over the impact of errors cascading up from small scales. A decomposition of the ensemble spread in components associated with the balanced and unbalanced dynamics shows that the initial uncertainties are primarily in the unbalanced modes, especially at the subsynoptic scales. The growth of uncertainties is found to be faster in the balanced than in the unbalanced modes and after 0.5–1 day of forecasts the balanced errors become dominant except at the subsynoptic scales.


Introduction
In spite of the progress in global numerical weather prediction (NWP) (e.g. Magnusson and Källén, 2013) and the progress in the simulation of tropical variability (e.g. Bechtold et al., 2008;Vitart and Molteni, 2010), the tropics remain a region with the largest uncertainties in the initial state (analysis) for NWP (e.g. Park et al., 2008) as well as the largest forecast errors in the short range (e.g. Žagar et al., 2013).
Properties of the global forecast errors in the short range (3-12 h) are especially important as their spatial correlations are used in data assimilation to spread the impact of observations in model space (e.g. Kalnay, 2003). In the mid-latitudes, the largest forecast errors are known to develop in the storm track regions and are understood with the help of quasi-geostrophic dynamics. In the tropics, the equatorially trapped inertio-gravity waves play a role in dynamics at all scales, and they are characterized by small phase speeds that are comparable with those of the equatorial Rossby waves. Žagar et al. (2005, 2007, 2013) showed that the equatorial large-scale inertio-gravity waves contribute a large portion of the short-range forecast-error variance statistics. However, little is known about the impact of these waves on the growth of forecast errors in the tropics and their effect on the mid-latitudes. * Corresponding author. e-mail: nedjeljka.zagar@fmf. uni-lj.si The loss of predictability and limits of prediction skill with relevance for NWP have been estimated using the models and variables which describe mid-latitude dynamics. Limits of prediction skill have been estimated to vary from about two weeks in the 1980s (Lorenz, 1982;Dalcher and Kalnay, 1987) to around one month (Shukla, 1981) and beyond (i.e. 45 days in mid-latitudes, Kleeman, 2008). Operational NWP centres such as the European Centre for Medium-Range Weather Forecasts (ECMWF) are now running the combined medium-range and monthly forecasting systems that provide estimates of forecast uncertainty in planetary scales for one month ahead and up to day 15 for smaller scales (Vitart et al., 2008).
Most of studies of the prediction skill in NWP models have relied on the comparison of forecasts with analyses. Global forecast errors estimated in this way have usually been based on the average properties of the geopotential height of 500 hPa level in the Northern Hemisphere extra-tropics (e.g. Lorenz, 1982;Tribbia and Baumhefner, 2004). This is the middle troposphere level where divergence is minimal and the level where the barotropic vorticity equation, which has been used for the studies of the error growth in turbulent flows with the kinetic energy spectrum proportional to the −3 power of the horizontal wavenumber (Lorenz, 1969), best applies. Although the idealized studies such as the barotropic vorticity equation are of little practical importance they illustrate the fundamental reasons for the predictability loss as a consequence of the upscale growth of initially small errors in small scales. Consequently, much research related to the growth of forecast errors in NWP models focused on the inverse cascade picture of the error growth from small-scale uncertainties in the initial state for NWP. On the other hand, Rotunno and Snyder (2008) showed that quasi-geostrophic dynamics with a kinetic energy spectrum proportional to the −5/3 power of the wavenumber produces a rapid saturation of forecast errors in small scales. Tribbia and Baumhefner (2004) discussed that the role of cascading smallscale errors is primarily to perturb the baroclinically unstable scales in the extra-tropics that is followed by a rapid error growth in these scales and the loss of large-scale predictability.
The present study discusses the role of large scales in the scale-dependent growth of the global forecast errors in the stateof-the-art NWP systems. The large-scale errors in NWP models were previously presented by Dalcher and Kalnay (1987) who compared the ECMWF forecasts with their verifying analyses for various global wavenumbers and showed that the forecast errors grow from the start of the forecasts in all scales. Their results have been reaffirmed by the recent paper of Žagar et al. (2015) who analyzed 7-day long ensemble forecasts from the operational ensemble prediction system of ECMWF (ENS). They showed that analysis uncertainties, as simulated by the combination of the ensemble of analyses produced by 4D-Var data assimilation (EDA, Isaksen et al., 2010) and singular vectors , exist in all spatial scales and grow from the start of the forecasts in all horizontal and vertical scales. In other words, initial uncertainties are larger in large scales than in small scales, and grow to larger forecast uncertainties in large scales compared to small scales. This result is in accord with studies by Rotunno and Snyder (2008) and by Durran and Gingrich (2014) who stressed that the impact of initial errors in small scales in the simplified model of quasi-geostrophic dynamics may be relatively unimportant with respect to relatively small initial errors in large scales.
This paper extends the results of Žagar et al. (2015, hereafter ŽBT) to the 15-day forecasts of the ECMWF ENS and contrasts analysis and forecast uncertainties in the tropics and the midlatitudes. Results provide a scale-dependent distribution of analysis and forecast uncertainties associated with quasi-geostrophic and inertio-gravity dynamics. The latter is of special interest as it is little understood. Presented results support the idea that the upscale error cascades are not the key to the loss of predictability in the current global ensemble prediction systems.
The paper consists of four sections. Section 2 presents the scale decomposition of the model-level data from ECMWF ENS. Section 3 presents the scale-dependent properties of the ensemble spread of analyses and 15-day long forecasts in relation to dynamics. Discussion and conclusions are given in Section 4.

Modal decomposition of the ensemble prediction system of ECMWF
The ECMWF ensemble system is described in numerous papers (e.g. Buizza et al., 2008;Buizza and Leutbecher, 2015). The ENS data from 7-day long operational ensemble forecasts in December 2014 were used by ŽBT to present the scale-dependent growth of the ensemble spread. They showed that the growth of the ensemble spread (standard deviation) follows fairly well the growth of the root-mean-square-differences between the ensemble mean and the control analysis (forecast errors). This supports the use of the ensemble spread as a proxy for the forecasterror statistics. In addition, it provides statistics on the initial uncertainties.

Modal decomposition method
The methodology for the scale decomposition of the 3D global ensemble spread relies on the representation of the global baroclinic atmosphere in terms of M global shallow-water systems, each characterized by its own fluid depth for the horizontal motions that is known as the equivalent depth. This parameter, denoted D m , is a constant which couples the non-dimensional oscillations of the horizontal wind and geopotential height fields with the vertical structure equation discretized for J terrainfollowing levels.
The horizontal motions for a given vertical mode m (m = 1, . . . , M) are represented by a series of Hough harmonics which consist of the Hough vector functions in the meridional direction and waves in the longitudinal direction. The complex expansion coefficient χ k n (m), which describes the two wind components (u, v) and the geopotential height h, is obtained by the following formula: Equation (1) consists of the vertical projection (within the parenthesis) followed by the horizontal step. The vertical structure functions G m ( j) are orthogonal and solved using the finite difference method for the model's temperature and stability profiles. The vertical structure of the first vertical mode, m = 1, has no zero-crossing point in the vertical and it is referred to as barotropic mode. The higher vertical modes, (m > 1), are referred to as baroclinic modes. Each vertical mode (m) is associated with (m−1) zero crossings in the vertical profile. The largest amplitude of the corresponding G m functions moves downward towards the bottom level as m increases. Examples of the vertical structure functions are presented in ŽBT.
The Hough harmonics are denoted by H k n (m); for every given vertical mode the Hough harmonics are characterized by the two indices for the zonal wavenumber k and meridional mode n. The scaling matrix S m is a 3 × 3 diagonal matrix which makes the input data vector after the vertical projection dimensionless. A discrete solution of (1) is obtained by replacing the integration by a finite series of the Hough harmonic functions including the zonally averaged state, K zonal waves and R meridional modes. The maximal number of meridional modes R combines N R balanced modes, denoted BAL, N E eastward-propagating inertio-gravity modes, denoted EIG and N W westward-propagating inertiogravity modes which are denoted WIG; R = N R + N E + N W . Parameters λ, ϕ stand for the geographical longitude and latitude, respectively, where μ = sin(ϕ). The details of the projection procedure are outlined in ŽBT and references therein. Note that ŽBT had used 'ROT' to denote balanced modes due to their predominantly rotational nature. However, as there is a rotational component in the inertio-gravity circulation, especially in the tropics, and a small component of divergence in the balanced flow, the balanced modes are here denoted 'BAL'.

Scale-dependent ensemble spread
The computation of the ensemble spread is carried out as follows. Each member of the operational 50-member ensemble, 00 UTC run, is projected following Equation (1) every 12 h during the 15-day period. For every time instant, the ensemble spread is computed in two steps. First, the ensemble mean χ k n (m) is computed for every mode from P = 50 samples as χ k n (m) = 1/P P p=1 χ k n (m; p). This is followed by the computation of the ensemble spread k n (m) according to the following formula: The unit of k n (m) is ms −1 . As shown in ŽBT, the ensemble spread computed by (2) is equivalent to the spread computed at the point (λ i , ϕ j , m) of physical space for the two wind components and geopotential normalized by D m for the 3D baroclinic model atmosphere represented in terms of M shallowwater layers by the vertical transform procedure. The applied spectral truncations are K = 120 zonal wavenumbers and N R = N E = N W = 60 meridional modes for each of the three species of oscillations. In the vertical direction, we use M = 50 modes. The data-set consists of 31 samples including initial states and forecasts every 12 h during 15 days (the length of operational ENS).
The 3D orthogonality of the NMF expansion allows us to consider the spread in every scale and mode independently.
Once Equation (2) has been evaluated, the ensemble spread can be integrated independently for the three motion types (BAL, EIG and WIG) along any of the three spatial directions. Three characteristic regimes, the planetary, synoptic and subsynoptic scales, are defined following Jung and Leutbecher (2008): the zonal wavenumbers 0 ≤ k ≤ 3 describe the planetary scales (denoted P S), wavenumbers 4 ≤ k ≤ 14 define synoptic scales (denoted SS) and all k ≥ 15 belong to the subsynoptic scales (denoted M S for mesoscale). In the mid-latitudes the zonal scales under 1000 km thus belong in the subsynoptic range whereas in the tropics the subsynoptic range starts at around 1300 km. We denote the integrated spread by E, and compute its components associated with the balanced modes as for the planetary, synoptic and subsynoptic (or mesoscale) range, respectively. Equivalent expressions apply for the spread associated with the eastward-propagating and westward-propagating inertio-gravity modes, E eig and E wig , respectively. The summation of the three spread components in (3) provides the total spread associated with the balanced modes, denoted E bal = E P S bal + E SS bal + E M S bal . The globally integrated total spread in the ENS system is then E = E bal + E eig + E wig and it can be split into the global planetary, synoptic and subsynoptic components, E P S , E SS , and E M S , respectively. The inertiogravity component is E ig = E eig + E wig .

ECMWF ensemble data
The present study is based on the 15-day operational ensemble forecasts in May 2015, Cycle 41r1. The ENS data are produced with a variable resolution; horizontal resolution T L 639 is applied during the first 10 days and T L 319 (i.e. about 65 km spacing) thereafter. The data are used on the N64 regular Gaussian grid (256 × 128 points). The top model levels which are known to show some systematic errors are excluded (see ŽBT) so that the decomposition (1) is carried out for 73 model levels under 10 hPa (i.e. model levels 19 to 91).
Some important properties of the ECMWF ENS ensemble spread in May 2015 are illustrated in Figs. 1 and 2 for the zonal wind spread. First, Fig. 1 presents the average spread at two model levels in 24-h and 5-day forecasts. The presented levels are near 150 and 300 hPa; the latter is the level with the maximal spread in the mid-latitudes whereas the former is close to the tropical tropopause. Figure 1 illustrates that the largest spread in 24-h forecasts is found in the upper troposphere at the latitude of the inter-tropical convergence zone north of the equator. The zonally averaged spread at three forecast steps, 12, 24 and 120 h is shown in Fig. 2. It confirms that the previously reported structure of the simulated 12-h forecast errors in the ECMWF system (Žagar et al., 2013) applies also to the present data; that is, a larger component of the ensemble spread is in the tropics than in the mid-latitudes at all levels and especially in the upper troposphere. The structure of the initial spread is the same as for 12-h forecasts, but the amplitude is smaller (not shown). The spread is maximal at levels around 140 hPa in the tropics. After 36 h of forecasts, the maximum of the zonally averaged spread is located in the mid-latitudes in relation to the growth of spread associated with developing baroclinic waves on the westerly flow. The amplitude of the mid-latitude spread continues growing to exceed the tropical spread in the 5-day range by a factor of two and more (note different scales of the colourbars in Figs. 1 and 2). The rest of the paper deals with the scale-dependent properties of the presented statistics in relation to the dynamics. The 'errors of the day', which may not be well predicted in individual ensemble runs, are not studied here. Furthermore, the statistics can be considered less representative for the tropics than for the mid-latitudes as discussed by ŽBT. The presented discussion relies on the ensemble spread being an indicator of the predictability loss in operational NWP ensembles.

Scale-dependent distribution of the global analysis and forecast uncertainties
An important property of the modal decomposition is that it is global and three-dimensional. This prevents us from a direct comparison with results valid for a single level in the extratropics that is a more usual approach. For example, if the spread   along a mid-latitude circle in Fig. 1 is decomposed using 1D Fourier decomposition at different forecast lengths, it would show a shift of the initial spread from the synoptic towards planetary scales (e.g. Satterfield and Szunyogh, 2011). As our decomposition includes both the tropical and the mid-latitude belts, we cannot differentiate between them in modal space. Discussion instead focuses on the balanced and IG components in relation to the tropical and mid-latitude large-scale dynamics. For a comparison with physical space, various components of the spread can be filtered from modal to physical space as illustrated in ŽBT.

Two-dimensional structure
Scale-dependent properties are presented in Figs. 3-9 by showing separately the balanced and IG components of the spread for the barotropic mode (m = 1, Fig. 3), for vertical mode m = 13 (Fig. 4) and for the vertically integrated total balanced spread (Fig. 5).
Forecast uncertainties in the balanced modes are largest in the synoptic horizontal scales with deepest vertical structure (smallest m). This is shown in Fig. 3 for m = 1 that can be considered a modal equivalent of the global mid-troposphere features including winds and geopotential height coupled through the    linear balance equation. An important difference to the statistics based on the mid-latitude data-sets is that a large part of the spread in Fig. 3a (initial state) is contributed by analysis uncertainties in the large-scale tropical winds that project to the lowest meridional modes. After 24 h, the maximum of the spread is shifted to larger n (Fig. 3b), in agreement with a larger growth of the spread in the mid-latitudes than in the tropics in Fig. 2b. In 7-day forecasts, the spread is largest in n = 5 (Fig. 3c) that corresponds to the large forecast uncertainties in the storm-track regions of the mid-latitudes (Fig. 2c). There is also a shift from zonal wavenumber k = 6−8 in initial states (Fig. 3a) to k = 3−5 in 7-day forecasts (Fig. 3c).
The amplitude of analysis and forecast uncertainties reduces as m increases since the vertical depth reduces and the meridional structure of the associated Hough functions becomes stronger bounded to the tropics. For example, vertical mode m = 13 shown in Fig. 4 has D 13 = 25 m which can not represent input data from high latitudes. Such small equivalent depths have been extensively used to characterize equatorial waves based on linear wave theory on the equatorial-β plane (e.g. Matsuno, 1966;Hayashi, 1981). Žagar et al. (2005, 2007) showed that such equivalent depths provide a successful application of the equatorial wave theory for the representation of the tropical forecast errors. Compared to m = 1, analysis uncertainties in m = 13 in Fig. 4 are located in larger scales (Fig. 4a) in agreement with Fig. 1. In longer forecasts, the spread in small m shifts to larger n that corresponds to latitudes further away from the equator (Fig. 4c).
With 50 vertical modes summed up, the global initial uncertainties in the balanced modes are maximal in the leading meridional mode n = 1 and zonal wavenumbers k = 1 − 6 ( Fig. 5a). Notice also a nearly isotropic nature of the initial balanced spread in Figs. 3a, 4a and 5a beyond zonal wavenumber k = 15 (about 1300 km at the equator, and about 900 km in midlatitudes). The tropical balanced spread is significantly smaller than the extra-tropical spread in 7-day forecasts, but the sum of the spread in the tropics (the region ±30 • off the equator represents the half of the atmosphere) and the mid-latitudes results in the maximal value of the globally integrated balanced spread in the zonal wavenumber k = 1 (Fig. 5c).
Analysis uncertainties in the tropics reside in large horizontal scales as illustrated in Fig. 1. These uncertainties largely project to the IG modes with the lowest meridional modes of different m (Figs. 6-9). For the eastward IG modes (EIG), the mode with the largest spread is the n = 0 mode that represents the Kelvin wave. Its analysis uncertainty is larger in the vertical mode m = 13 than m = 1 since the former corresponds to the Kelvin waves with smaller phase speeds that are coupled to convection (Fig. 7a versus Fig. 6a). During forecasts, the spread in the equatorially trapped Kelvin waves maximizes in the synoptic scales related to the organized convection generating the Kelvin waves (Fig.  7c). These are also the scales and the mode identified in ŽBT as ones that lack a variability in the ensemble prediction system. The vertically integrated EIG spread is dominated by the zonal wavenumber k = 1 Kelvin mode associated with the spread in the tropical stratosphere (not shown).
The spread associated with the westward IG modes (WIG) has on average smaller amplitudes than in the case of the EIG modes (Figs. 8 and 9). In physical space, the WIG spread is found also in the mid-latitudes that is associated with the forecast errors in baroclinic waves projecting partly to the IG modes that propagate eastward. Therefore, the horizontal distribution of the WIG spread is characterized by larger amplitudes in higher n and larger k than in the case of EIG modes in both barotropic mode (Fig. 8) and in the vertically integrated spread (Fig. 9). In particular, notice a region of larger spread in Fig. 9c in synoptic scales in 7-day forecasts extended towards n = 5 that is coupled to the balanced spread in the same modes in Fig. 3c.

One-dimensional growth
In this subsection, we look in detail at the average spread growth in individual modes and the growth of integrated spread as a function of the zonal wavenumber. First we show in Fig. 10 the growth of spread in several zonal wavenumbers during the 360-h forecast. Presented n = 1 mode is one of the most studied modes of the global atmosphere. The spread in n = 1 is shown for a single vertical mode (Fig. 10a) and summed up for all vertical modes representing the model depth up to 10 hPa (Fig. 10b). Figure 10a and b shows that the spread in small scales saturates quickly while the spread in planetary and synoptic scales steadily increases. In particular, the spread in synoptic scales increasingly projects to the barotropic mode (m = 1). A quick saturation of the spread in small scales is associated with the distribution of initial perturbations in the leading vertical mode that contains most of the variance in large zonal scales. In larger meridional modes, which contain larger initial perturbations in the subsynoptic scales such as mode n = 10 shown in Fig. 10c, the spread in the smallest scales is quickly growing in forecasts, especially in deep modes. In Fig. 10, one can also notice a bump at day 10 that is related to the resolution upscale of the ENS system. Figure 10a and b is an effective illustration of the scaledependent amplitude of forecast uncertainties. It shows that the larger the zonal scale, the greater the amplitude of the ensemble spread. At initial time, this applies to all zonal wavenumbers and vertical modes. Later in the forecasts, the growth of spread in synoptic scales in the upper troposphere, which largely projects on the barotropic vertical mode, dominates over the growth in planetary scales (Fig. 10a). Overall Fig. 10 suggests that in the global NWP framework, the upscale error cascades are not the most dominant process for the error growth. In the following figures, we shall elaborate this result in more detail focusing on the role of the tropics and the relative importance of the balanced and IG components.     Figures 11 and 12 show the evolution of spread in several meridional and vertical modes for the balanced and IG components. In almost all cases, the large-scale uncertainties are initially dominant. After the first 12 h of forecasts, spread in synoptic scales in the barotropic mode significantly grows in comparison to the planetary scales, especially in n = 1 (Fig.  11b). Figure 11b can be compared to Fig. 3 and to Fig. 2 in physical space and it is similar to the results obtained by other analysis methods usually applied to the mid-latitude data (e.g. Satterfield and Szunyogh, 2011). At the same time, Fig. 11 shows that the spread grows also in planetary scales from the start of the forecasts. In higher vertical modes (e.g. m = 4) that are more associated with the tropics the spread in planetary and synoptic scale grow equally fast (Fig. 11c). An exception is the n = 0 balanced mode, the mixed Rossby-gravity mode shown in Fig. 11a, that is dominated by the spread in zonal wavenumbers k = 4 − 6 in the tropics.
The growth of spread in two vertical modes of the Kelvin wave is shown in Fig. 12a and b. The barotropic mode m = 1 contains a component of the spread from each model level whereas the mode m = 4 has the largest amplitude in the stratosphere. The growth of spread greatly increases in m = 4 after the second day of forecasts, in agreement with the picture of the tropical forecast errors that propagate vertically to the stratosphere while amplifying (Žagar et al., 2016). A recent study by Žagar et al. (2016) showed that after the first two days of forecasts the tropical forecast errors in the ECMWF model project primarily to the balanced modes, especially in the stratosphere. This reaffirms conclusions that the planetary-scale spread is largely contributed by the tropical region, whereas the synoptic-scale maximum of the spread in both balanced (Fig. 11) and WIG modes (Fig. 12c) is more associated with the mid-latitudes processes. Notice also that the spread 'bump' in synoptic scales in Fig. 11 would be larger if the data were from the winter season as typically used.
Figures 10-12 also show that the spread in the subsynoptic range appears nearly saturated after a couple of days. It saturates at much smaller value than estimated from climatological variance in ŽBT suggestive of a lack of mesoscale variability in the ensemble system. In contrast, the spread in the synoptic and planetary scale approaches its asymptotic values more slowly following the fast growth in the first two days, in agreement with the parametric models of the forecast error growth (e.g. Dalcher and Kalnay, 1987). In order to illustrate the relative growth across many scales, we present in Fig. 13 the scaledependent spread normalized by its values at the final step (360 h) in the same wavenumber for the balanced and IG modes. First we observe in Fig. 13a that the global growth of the spread appears relatively flat for the wavenumbers 1 to 5 and strongly scale-dependent for greater wavenumbers; in this range the exponential spread growth applies. The growth in the balanced modes appears similar to the total growth. The large-scale IG modes contain 40 to 50% of their 15-day spread in the initial states ( Fig. 13c and d) that is due to the already discussed large uncertainties in the initial state in the tropics. In contrast, initial uncertainties in the planetary and synoptic scales in the balanced modes have magnitudes only 10% of their values in the 15-day forecast range (Fig. 13b).

Spread partitioning between the balanced and inertio-gravity components
During the forecasts, the percentage of the balanced spread increases as shown in Fig. 14a. The balanced spread portion grows from 38% initially to about 67% of the total spread at day 15. In 15-day range, the IG component exceeds the balanced spread only in scales smaller than 500 km (k < 30). Initially, around 25% of the total spread is in the subsynoptic scales (k < 15) and about 50% in scales k < 30. As the forecasts evolve, a portion of the spread in large scales increases and after 15 days there is 55 and 75% of the spread in scales k < 15 and k < 30, respectively (integration of the black line in Fig. 14b). At this time, a crossing point between the balanced and IG spread can be seen close to zonal wavenumber 35 that corresponds to 400 km in the mid-latitudes and about 600 km in the tropics (Fig. 14b). Such crossing point was presented in ŽBT for the climatological variance at somewhat larger scale since at the 15-day forecast range the spread has not yet reached its asymptotic value in planetary scales. This can be seen more clearly in the next figure that shows the growth of the total global spread in the planetary, synoptic and subsynoptic scales split between the balanced component and IG component. Figure 15 confirms that the IG component is initially greater than the balanced spread in all scales, but especially in the subsynoptic range. However, the global growth of errors associated with quasi-geostrophic dynamics makes the balanced spread dominant after 12 h in planetary and after 24 h in the synoptic scales. In the subsynoptic range, the two components converge until forecast day 10 when the model horizontal resolution reduces. Equal levels of the balanced and IG spread in the subsynoptic range in Fig. 15d may be associated also with the lack of mesoscale variability in the ensemble as seen in Figs. 11 and 12 and discussed in ŽBT. The total global spread associated with the balanced modes becomes greater than the spread projecting to IG modes after two days of forecasts (Fig. 15a). In the 15-day forecast range used for the operational medium-range ensemble  prediction, the global spread, now mainly in the balanced largescale modes, appears not yet saturated in agreement with the recent estimate of the limits of global prediction skill by Buizza and Leutbecher (2015). As discussed in ŽBT, it is possible that the partitioning of spread between the balanced and IG parts in the initial state is associated with perturbed observations used in EDA, i.e. that initial perturbations for ensemble prediction are not optimal. It would be interesting to check if other ensemble prediction systems show the same properties of the global spread, especially the distribution of analysis uncertainties in the tropics.

Discussion and conclusions
We assessed the scale-dependent growth of forecast uncertainties in a state-of-the-art global forecasting system, the operational ECMWF ensemble prediction system. Presented results complement previous studies of the forecast error growth that focused on the extra-tropical quasi-geostrophic dynamics and often considered the error-free large-scale initial state. In contrast, the operational NWP systems are characterized by uncertainties in the initial state at all scales. The largest initial uncertainties and the largest initial growth of forecast errors is found to be in the tropics. The tropical error growth is associated with the model deficiencies and parametrization of deep convection (e.g. Reynolds et al., 1994) in addition to its great dependence on the simulated initial perturbations (Žagar et al., 2013, 2015). It is only after about 1.5-2 days into the forecasts when the forecast uncertainties in the mid-latitudes become a dominant feature of the global prediction.
Large-amplitude large-scale forecast errors early in the forecasts have been presented by Dalcher and Kalnay (1987) for the 500 hPa height. However, this important property of the forecast errors has not received much attention while the community focused on the error cascades from the smaller scales. The historical paper by Lorenz (1969) included the impact of the initial errors at large scales but it seems to have been overlooked in the conclusions and in much of the subsequent research (Durran and Gingrich, 2014). Today we have increasing evidence of significant uncertainties in the initial state for NWP at the large scales in the tropics (Park et al., 2008;Žagar et al., 2015). The importance of the initial large-scale errors in the tropics was recognized also by Mapes et al. (2008) in an aqua-planet simulation whereas Durran et al. (2013) made a similar point for a mesoscale NWP model in the mid-latitudes. Several studies described the structure of the short-term tropical forecast errors in relation to unbalanced dynamics, especially the largescale equatorial waves (Žagar et al., 2005, 2007, 2013). These properties are further discussed in the present paper to make evident the need to include tropical processes in the discussion of the limits of prediction skill in global NWP models. Tropical analysis uncertainties are associated with a lack of observations, especially wind observations and deficiencies in the tropical data assimilation and NWP modelling. The reduction of these uncertainties is needed for the practical predictability to converge towards the theoretical limit.
The presented results provide a scale-dependent quantification of the global forecast uncertainties associated with the balanced and inertio-gravity (or unbalanced) modes across many scales. In particular, it is shown that the larger the zonal scale, the greater the amplitude of analysis uncertainties in selected individual modes. Nevertheless, globally integrated analysis uncertainties reside mainly (around 75% of the global value) at the subsynoptic scales (zonal wavenumber 15 and greater). At the start of forecasts, the growth of forecast uncertainties takes place in all wavenumbers. Moreover, the growth of uncertainty at large scales appears dominant over the impact of errors cascading up from small scales; in other words, the upscale error cascade is not the key process for the growth of the global forecast uncertainties in the medium range.
The applied decomposition method shows that the analysis uncertainties, as represented by the combination of the ensemble data assimilation and singular vectors in the ECMWF ENS data, project to a greater part to the inertio-gravity modes than to the balanced modes in all zonal wavenumbers. A fast growth of the forecast uncertainties in the balanced modes exceeds initially dominant unbalanced uncertainties after 0.5-1 day except at the subsynoptic scales. At this scale range, the initial growth of the balanced forecast uncertainties is largest but the unbalanced component remains dominant during the two-week forecasts.