A- A+
Alt. Display

# 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.

Keywords:
How to Cite: Žagar, N., 2017. A global perspective of the limits of prediction skill of NWP models. Tellus A: Dynamic Meteorology and Oceanography, 69(1), p.1317573. DOI: http://doi.org/10.1080/16000870.2017.1317573
Published on 01 Jan 2017
Accepted on 23 Jan 2017            Submitted on 2 May 2016
1.

## 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.

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. Althoughthe 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 small-scale 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 state-of-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 (Buizza et al., 2008), 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 mid-latitudes. 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.

2.

## 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 forecast-error statistics. In addition, it provides statistics on the initial uncertainties.

2.1.

### 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 terrain-following levels.

The horizontal motions for a given vertical mode m ($m=1,\dots ,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 ${\mathit{\chi }}_{n}^{k}\left(m\right)$, which describes the two wind components (uv) and the geopotential height h, is obtained by the following formula:

((1) )
$\begin{array}{cc}\hfill {\mathit{\chi }}_{n}^{k}\left(m\right)& =\frac{1}{2\mathit{\pi }}{\int }_{0}^{2\mathit{\pi }}{\int }_{-1}^{1}\left({\mathbf{S}}_{m}^{-1}\sum _{j=1}^{J}{\left(u,v,h\right)}^{T}{G}_{m}\left(j\right)\right)\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}·{\left[{\mathbf{H}}_{n}^{k}\left(m\right)\right]}^{\ast }\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\mathit{\mu }\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\mathit{\lambda }.\hfill \end{array}$

Equation (1) consists of the vertical projection (within the parenthesis) followed by the horizontal step. The vertical structure functions ${G}_{m}\left(j\right)$ 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, $\left(m>1\right)$, are referred to as baroclinic modes. Each vertical mode (m) is associated with $\left(m-1\right)$ 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 ${\mathbf{H}}_{n}^{k}\left(m\right)$; 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 ${\mathbf{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 inertio-gravity modes which are denoted WIG; $R={N}_{R}+{N}_{E}+{N}_{W}$. Parameters $\mathit{\lambda }$, $\mathit{\phi }$ stand for the geographical longitude and latitude, respectively, where $\mathit{\mu }=sin\left(\mathit{\phi }\right)$. 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’.

2.2.

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 $\overline{{\mathit{\chi }}_{n}^{k}\left(m\right)}$ is computed for every mode from $P=50$ samples as $\overline{{\mathit{\chi }}_{n}^{k}\left(m\right)}=1/P{\sum }_{p=1}^{P}{\mathit{\chi }}_{n}^{k}\left(m;p\right)$. This is followed by the computation of the ensemble spread ${\mathrm{\Sigma }}_{n}^{k}\left(m\right)$ according to the following formula:

((2) )
$\begin{array}{cc}\hfill {\left[{\mathrm{\Sigma }}_{n}^{k}\left(m\right)\right]}^{2}& =\frac{1}{P-1}\sum _{p=1}^{P}g{D}_{m}\left({\mathit{\chi }}_{n}^{k}\left(m;p\right)-\overline{{\mathit{\chi }}_{n}^{k}\left(m;p\right)}\right)\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}×\left({\mathit{\chi }}_{n}^{k}\left(m;p\right)-\overline{{\mathit{\chi }}_{n}^{k}\left(m;p\right)}{\right)}^{\ast }.\hfill \end{array}$

The unit of ${\mathrm{\Sigma }}_{n}^{k}\left(m\right)$ is ${\mathrm{ms}}^{-1}$. As shown in ŽBT, the ensemble spread computed by (2) is equivalent to the spread computed at the point $\left({\mathit{\lambda }}_{i},{\mathit{\phi }}_{j},m\right)$ 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 shallow-water 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).

Figure 1.

Growth of the global forecast uncertainties. Ensemble spread in May 2015 ECMWF ENS for zonal wind at (a–b) model level 45 (around 150 hPa) and (c–d) model level 55 close to 300 hPa in (a,c) 24-h forecasts and (b,d) 120-h forecasts. Colourbar is in m/s.

Figure 2.

As in Fig. 1 but for the zonally averaged zonal wind spread in (a) 12-h forecast, (b) 24-h forecast and (c) 120-h forecast.

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\le k\le 3$ describe the planetary scales (denoted PS), wavenumbers $4\le k\le 14$ define synoptic scales (denoted SS) and all $k\ge 15$ belong to the subsynoptic scales (denoted MS 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

((3) )
$\begin{array}{cc}\hfill {E}_{bal}^{PS}& =\sum _{m=1}^{M}\sum _{n=0}^{{N}_{R}}\sum _{k=0}^{3}{\mathrm{\Sigma }}_{n}^{k}\left(m\right),\hfill \\ \hfill {E}_{bal}^{SS}& =\sum _{m=1}^{M}\sum _{n=0}^{{N}_{R}}\sum _{k=4}^{14}{\mathrm{\Sigma }}_{n}^{k}\left(m\right),\hfill \\ \hfill {E}_{bal}^{MS}& =\sum _{m=1}^{M}\sum _{n=0}^{{N}_{R}}\sum _{k=15}^{K}{\mathrm{\Sigma }}_{n}^{k}\left(m\right),\hfill \end{array}$

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}_{bal}^{PS}+{E}_{bal}^{SS}+{E}_{bal}^{MS}$. 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}^{PS}$, ${E}^{SS}$, and ${E}^{MS}$, respectively. The inertio-gravity component is ${E}_{ig}={E}_{eig}+{E}_{wig}$.

2.3.

### 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).

Figure 3.

Ensemble spread as a function of the zonal wavenumber and meridional mode for the balanced vertical mode $m=1$ at (a) initial time, (b) 24-h forecast, and (c) 7-day forecast. Colourbar is in meters per second.

Figure 4.

As in Fig. 3 but for the balanced vertical mode $m=13$ with the equivalent depth ${D}_{m}$ around 25 m.

Figure 5.

As in Fig. 3 but the sum of the balanced spread in all vertical modes.

Figure 6.

As in Fig. 3 but for the EIG spread.

Figure 7.

As in Fig. 4 but for the EIG spread.

Figure 8.

As in Fig. 3 but for the WIG spread.

Figure 9.

As in Fig. 5 but for the WIG spread.

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.

3.

## 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 extra-tropics 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.

3.1.

### Two-dimensional structure

Scale-dependent properties are presented in Figs. 39 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).

Figure 10.

The growth of estimated forecast errors in selected meridional modes and zonal scales as a function of the zonal wavenumber and vertical mode (as specified in the legend). (a) meridional mode $n=1$ and vertical mode $m=1$, (b) meridional mode $n=1$ and all vertical modes integrated and (c) meridional mode $n=10$ and zonal wavenumber $k=120$. The presented quantity is $\mathrm{l}og\left[E\left(t\right)\right]$ integrated over individual modes as specified in the legends.

Figure 11.

The growth of ensemble spread in selected balanced modes as a function of the zonal wavenumber. The bottom curve in each panel shows spread at the initial time whereas every next curve applies to 12-h longer forecasts. (a) meridional mode $n=0$ and vertical mode $m=1$, (b) meridional mode $n=1$ and vertical mode $m=1$, and (c) meridional mode $n=1$ and vertical mode $m=4$.

Figure 12.

As in Fig. 11 but for the IG modes. (a) EIG mode with meridional index $n=0$ and the first vertical mode $m=1$, (b) EIG mode with meridional index $n=0$ and vertical mode $m=4$, and (c) WIG meridional mode $n=0$ and vertical mode $m=4$.

Figure 13.

Ensemble spread as a function of the zonal wavenumber normalized by the spread at the forecast range 360 h in the same zonal wavenumber. (a) Total spread, $E\left(t\right)/E\left(t=360$ h), (b) balanced spread, ${E}_{bal}\left(t\right)/{E}_{bal}\left(t=360$ h), (c) eastward-propagating inertio-gravity spread, ${E}_{eig}\left(t\right)/{E}_{eig}\left(t=360$ h) and (d) westward-propagating inertio-gravity spread, ${E}_{wig}\left(t\right)/{E}_{wig}\left(t=360$ h). The axes are logarithmic.

Figure 14.

(a) Balanced ensemble spread as a function of the zonal wavenumber normalized by the total spread at the same forecast range and zonal wavenumber, ${E}_{bal}\left(t\right)/E\left(t\right)$, (b) Portion of the global spread at initial time and 15-day range in balanced and IG components in each zonal wavenumber. Different colour lines are specified in legend. Blue lines are (dashed) ${E}_{bal}\left(t=0\right)/E\left(t=0\right)$ and (full) ${E}_{bal}\left(t=360\right)/E\left(t=360\right)$, red lines represent (dashed) ${E}_{ig}\left(t=0\right)/E\left(t=0\right)$ and (full) ${E}_{ig}\left(t=360\right)/E\left(t=360\right)$. Black line corresponds to the portion of the total spread in each zonal wavenumber at 2-week range multiplied by 10.

Figure 15.

The scale-dependent growth of estimated forecast errors in the balanced and inertio-gravity regimes. The presented quality is the global error integrated over (a) all scales ($\mathrm{l}og\left[E\right]$, $\mathrm{l}og\left[{E}_{bal}\right]$, $\mathrm{l}og\left[{E}_{ig}\right]$), (b) planetary scales ($\mathrm{l}og\left[{E}^{PS}\right]$, $\mathrm{l}og\left[{E}_{bal}^{PS}\right]$, $\mathrm{l}og\left[{E}_{ig}^{PS}\right]$), (c) synoptic scale ($\mathrm{l}og\left[{E}^{SS}\right]$, $\mathrm{l}og\left[{E}_{bal}^{SS}\right]$, $\mathrm{l}og\left[{E}_{ig}^{SS}\right]$) and (d) subsynoptic scales ($\mathrm{l}og\left[{E}^{MS}\right]$, $\mathrm{l}og\left[{E}_{bal}^{MS}\right]$, $\mathrm{l}og\left[{E}_{ig}^{MS}\right]$).

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-$\mathit{\beta }$ 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. 4(c)).

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 mid-latitudes). 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${}^{\circ }$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. 5(c)).

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. 69). 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. 7(c)). 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. 9(c) in synoptic scales in 7-day forecasts extended towards $n=5$ that is coupled to the balanced spread in the same modes in Fig. 3(c).

3.2.

### 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. 10(b)). 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. 10(c), 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 scale-dependent 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. 11(b)). Figure 11(b) 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. 11(c)). 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. 12(c)) 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 1012 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 scale-dependent 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. 13(c) 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. 13(b)).

3.3.

### 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. 14(b)). 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. 14(b)). 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. 15(d) 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 large-scale 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.

4.

## 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 large-scale 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.

## Disclosure statement

1No potential conflict of interest was reported by the author.

## Acknowledgements

The input data were provided by the European Centre for Medium Range Weather Forecasts (ECMWF), http://www.ecmwf.int. The author is very grateful to Prof. Eugenia Kalnay and Dr Peter Bechtold for their constructive criticism and to Blaž Jesenko for drawing Figures 1-2.

## References

1. Bechtold , P. , Köhler , M. , Jung , T. , Doblas-Reyes , F. , Leutbecher , M. , et al. . 2008 . Advances in simulating atmospheric variability with the ECMWF model: from synoptic to decadal time-scales . Q. J. R. Meteorol. Soc. 134 , 1337 – 1351 .

2. Buizza , R. and Leutbecher , M. 2015 . The forecast skill horizon . Q. J. R. Meteorol. Soc 141 , 3366 – 3382 .

3. Buizza , R. , Leutbecher , M. and Isaksen , L. 2008 . Potential use of an ensemble of analyses in the ECMWF Ensemble Prediction System . Q. J. R. Meteorol. Soc. 134 , 2051 – 2066 .

4. Dalcher , A. and Kalnay , E. 1987 . Error growth and predictability in operational ECMWF forecasts . Tellus 39A , 474 – 491 .

5. Durran , D. R. and Gingrich , M. 2014 . Atmospheric predictability: why butterflies are not of practical importance . J. Atmos. Sci. 71 , 2476 – 2488 .

6. Durran , D. R. , Reinecke , P. A. and Doyle , J. D. 2013 . Large-scale errors and mesoscale predictability in Pacific Northwest snowstorms . J. Atmos. Sci. 70 , 1470 – 1487 .

7. Hayashi , J. 1981 . Space-time spectral analysis and its applications to atmospheric waves . J. Meteorol. Soc. Japan 60 , 156 – 171 .

8. Isaksen , L. , Bonavita , M. , Buizza , R. , Fisher , M. , Haseler , J. , et al. 2010 . Ensemble of data assimilations at ECMWF . ECMWF Technical Memorandum No. 636. (Available from ECMWF, Shinfield Park, Reading RG2-9AX, UK) . Online at: http://old.ecmwf.int/publications .

9. Jung , T. and Leutbecher , M. 2008 . Scale-dependent verification of ensemble forecasts . Q. J. R. Meteorol. Soc 134 , 973 – 984 .

10. Kalnay , E. 2003 . Atmospheric Modeling, Data Assimilation and Predictability Cambridge University Press , 341 p .

11. Kleeman , R. 2008 . Limits, variability, and general behaviour of statistical predictability of the midlatitude atmosphere . J. Atmos. Sci. 65 , 263 – 275 .

12. Lorenz , E. 1969 . The predictability of a flow which possess many scales of motion . Tellus XXI , 289 – 307 .

13. Lorenz , E. 1982 . Atmospheric predictability experiments with a large numerical model . Tellus 34 , 505 – 513 .

14. Mapes , B. , Tulich , S. , Nasuno , T. and Satoh , M. 2008 . Predictability aspects of global aqua-planet simulations with explicit convection . J. Meteorol. Soc. Japan 86A , 175 – 185 .

15. Magnusson , L. and Källén , E. 2013 . Factors influencing skill improvements in the ECMWF forecasting system . Mon. Wea. Rev. 141 , 3142 – 3153 .

16. Matsuno , T. 1966 . Quasi-geostrophic motions in the equatorial area . J. Meteorol. Soc. Japan 44 , 25 – 42 .

17. Park , Y.-Y. , Buizza , R. and Leutbecher , M. 2008 . TIGGE: preliminary results on comparing and combining ensembles . Q. J. R. Meteorol. Soc 134 , 2029 – 2050 .

18. Rotunno , R. and Snyder , C. 2008 . A generalization of Lorenz’s model for the predictability of flows with many scales of motion . J. Atmos. Sci. 65 , 1063 – 1076 .

19. Reynolds , C. A. , Webster , P. J. and Kalnay , E. 1994 . Random error growth in NMC’s global forecasts . Mon. Wea. Rev. 122 , 1281 – 1305 .

20. Satterfield , E. and Szunyogh , I. 2011 . Assessing the performance of an ensemble forecast system in predicting the magnitude and the spectrum of analysis and forecast uncertainties . Mon. Wea. Rev. 139 , 1207 – 1223 .

21. Shukla , J. 1981 . Dynamical predictability of monthly means . J. Atmos. Sci. 38 , 2547 – 2572 .

22. Tribbia , J. and Baumhefner , D. P. 2004 . Scale interactions and atmospheric predictability: an updated perspective . Mon. Wea. Rev. 132 , 703 – 713 .

23. Vitart , F. , Buizza , R. , Alonso Balmaseda , M. , Balsamo , G. , Bidlot , J.-R. , et al. . 2008 . The new VAREPS-monthly forecasting system: a first step towards seamless prediction . Q. J. R. Meteorol. Soc 134 , 1789 – 1799 .

24. Vitart , F. and Molteni , F. 2010 . Simulation of the Madden-Julian Oscillation and its teleconnections in the ECMWF forecast system . Q. J. R. Meteorol. Soc 136 , 842 – 855 .

25. Žagar , N. , Andersson , E. and Fisher , M. 2005 . Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors . Q. J. R. Meteorol. Soc. 131 , 987 – 1011 .

26. Žagar , N. , Andersson , E. , Fisher , M. and Untch , A. 2007 . Influence of the quasi-biennial oscillation on the ECMWF model short-range forecast errors in the tropical stratosphere . Q. J. R. Meteorol. Soc. 133 , 1843 – 1853 .

27. Žagar , N. , Blaauw , M. , Jesenko , B. and Magnusson , L. 2016 . Diagnosing model performance in the tropics . ECMWF Newsletter 147 , 26 – 33 .

28. Žagar , N. , Buizza , R. and Tribbia , J. 2015 . A three-dimensional multivariate modal analysis of atmospheric predictability with application to the ECMWF ensemble . J. Atmos. Sci. 72 , 4423 – 4444 .

29. Žagar , N. , Isaksen , L. , Tan , D. and Tribbia , J. 2013 . Balance and flow-dependency of background-error variances in the ECMWF 4D-Var ensemble . Q. J. R. Meteorol. Soc. 139 , 1229 – 1238 .