Estimating trajectory uncertainties due to flow dependent errors in the atmospheric analysis

Atmospheric trajectory calculations is a widely used method to determine the origin or fate of air parcels in the atmosphere. For example, the use of back trajectories is common in atmospheric chemistry to determine source regions of sampled air at a particular site. The uncertainty of a calculated trajectory is dependent on the uncertainty in the atmospheric analysis. The errors in meteorological analyses are partly due to observation errors and the fact that observations can include local phenomena that are not representative for an entire grid box. In order to produce a best estimate of the atmospheric state, numerical weather prediction centers use advanced data assimilation systems. In data assimilation, the information from observations combined with a background (usually a short forecast valid at the time of the analysis) is used to obtain an analysis. However, both of these components contain uncertainties. Thus the resulting analysis will not perfectly match the current atmosphere and will therefore impact on the calculated trajectory path.


Introduction
Atmospheric trajectory calculations is a widely used method to determine the origin or fate of air parcels in the atmosphere.For example, the use of back trajectories is common in atmospheric chemistry to determine source regions Correspondence to: A. Engström (anderse@misu.su.se) of sampled air at a particular site (see e.g.Verver et al., 1999).Furthermore, it is an efficient way to couple individual measurements, on e.g.aerosol composition, with the general characteristics of a source region.Although trajectory calculations is a powerful tool in atmospheric research, the calculations contain errors that must be accounted for.To this end, numerous studies have investigated the impact from e.g.spatial and temporal interpolation errors, calculation method errors and wind field analysis differences (see Stohl, 1998, and references therein).One of the largest errors in trajectory calculations originate in differences in meteorological data sets.Harris et al. (2005) found, when comparing different contributions to trajectory uncertainty, that using different reanalysis data sets resulted in an uncertainty of 30-40% of the average distance traveled.However, it was also concluded that no data set could be considered superior to another despite that differences existed in the resulting trajectories (Harris et al., 2005).
To calculate trajectories information on the atmospheric motion is needed.In the case of forecast trajectories, information from a numerical weather prediction model is used.Errors in such trajectory calculations will to a large extent be subject to uncertainties in the predicted wind field.To calculate backward trajectories, one can instead use atmospheric analysis data, for instance the two widely used re-analysis data sets: the ECMWF ERA-40 Reanalysis Data Set (Uppala et al., 2005) or the NCEP/NCAR Reanalysis Data Set (Kalnay et al., 1996).Even if an analysis is the best estimate of the true state of the atmosphere it still contains errors.One widely used method to account for such errors is to calculate many trajectories, each with a perturbed receptor/emission point to obtain an estimate of the uncertainty in the area of origin (Merrill et al., 1985).However, this method may not accurately sample the uncertainty in the analysis since it has not been investigated whether an uncertainty in the atmospheric analysis can be represented using this approach.
A. Engström and L. Magnusson: Trajectory uncertainties due to analysis errors The errors in meteorological analyses are partly due to observation errors and the fact that observations can include local phenomena that are not representative for an entire grid box.In order to produce a best estimate of the atmospheric state, numerical weather prediction centers use advanced data assimilation systems.In data assimilation, the information from observations combined with a background (usually a short forecast valid at the time of the analysis) is used to obtain an analysis.However, both of these components contain uncertainties.Thus the resulting analysis will not perfectly match the current atmosphere.Even if one knows that errors are present, it is difficult to estimate the error amplitude and to simulate the corresponding error structures.For instance, one known property of the analysis uncertainty is the existence of spatial correlations of, and correlations between, meteorological variables.The error in one grid point is related to the error in adjacent grid points.Furthermore, the structure of the correlations and the amplitude of the errors vary on a day-to-day basis due to the current flow situation.In data assimilation, this error correlation is obtained using differences between forecasts of different length or with flow dependent covariance matrices obtained from an ensemble of short forecasts (Bannister, 2008).
The effect of the analysis error on the forecast uncertainty is simulated in ensemble weather forecasting.Due to the chaotic behavior of the atmosphere, small uncertainties in the analysis will grow with increasing forecast length.In order to simulate these analysis uncertainties several methods are proposed in the literature.These are e.g.singular vectors (Palmer, 1993), breeding vectors (Toth and Kalnay, 1997), Ensemble Transform method (Wei et al., 2008) or the Ensemble Transform Kalman filer (Wang and Bishop, 2003).The common purpose of these methods is to add perturbations to the estimated analysis.Another viable alternative is to perform multiple analyses, ensemble data assimilation (Houtekamer and Mitchell, 1998;Buizza et al., 2008).In Magnusson et al. (2009) the properties of the singular vector method and the Ensemble Transform method were compared.It was shown that the Ensemble Transform perturbations are closer to the expected properties of the analysis error compared to singular vectors.The Ensemble Transform method is a development of the breeding method (Toth and Kalnay, 1997) and uses the previous ensemble to create new perturbations with the purpose to sample perturbation structures that have grown in the past due to a sensitive flow situation.The perturbations are expected to have a spatial and multivariate structure that is similar to correlations in the analysis error and have uncertainties that are dependent on the current flow situation.
In this study we will compare two different methods to generate a set of trajectories that can be used to estimate the spatial uncertainty, due to analysis errors, in a calculated trajectory.The first method is the Ensemble analysis method (EA) where all trajectories are initiated at the same point, but each trajectory path is calculated from a perturbed anal-ysis used in weather prediction to sample the analysis uncertainty.The ensemble of analyses is created with the Ensemble Transform method that is used to perturb the analysis from the ECMWF data assimilation system.The second method is the Initial spread method (IS) where a set of trajectories, each with an individual receptor point, is used to sample uncertainties in the atmospheric flow.This is the most common method used in atmospheric chemistry applications since it allows for a large number of trajectories to be computed from a high-resolution analysis and still accounting for errors by perturbing the receptor point.
The use of the ECMWF ensemble prediction system to estimate forecast trajectory uncertainty was studied in Scheele and Siegmund (2001) and Straume (2001) where a trajectory was started from each member of the ensemble.Straume (2001) found that the ensemble spread is mainly important for long trajectory calculations (more than 24 h) and may not be suitable for shorter estimates of trajectory uncertainties.However, both Scheele and Siegmund (2001) and Straume (2001) studied forecast trajectories based on the ensembles created with the singular vector method.The singular vector method is aimed at finding structures that will grow most rapidly during a period of optimization.However, as mentioned above, these structures might not necessarily resemble the expected properties of the analysis error.Therefore, in the case of backward trajectories where analysis data is used, other methods might be more appropriate.
This paper is organized as follows.In Sect. 2 we describe the trajectory calculations and the two methods that we have used in this study.In Sect. 3 the results from two selected case studies are presented together with statistics for several different regions.Finally, in Sect. 4 the conclusion are presented.

Methods
To investigate the difference between the two methods and to study how the analysis error impacts on the calculated back trajectory path, daily five-day back trajectories are calculated for for 15 different locations around the globe.The period selected in this study spans 29 subsequent days between July and August 2005 and 25 subsequent days between December 2005 and January 2006.The start-point locations are 60 • S, 30 • S, 0 • N, 30 • N and 60 • N for 120 • W, 0 • E and 120 • E. By selecting these regions we can compare the two methods subject to different meteorological conditions.
We use a kinematic three-dimensional trajectory model available at the ECMWF where the trajectory path is determined by the full three-dimensional wind field (P.Kållberg, personal communication, 2009).An analysis field (resolution 1.125 • ×1.125 • and 40 vertical levels) is provided every 6 h and the wind field is interpolated linearly between each analysis.The provided analysis field is interpolated from the operational analysis which uses a resolution of 0.5 • ×0.5 • and 60 vertical levels.The trajectory calculation time step is 30 min.All trajectories are released at either the 850 hPa or the 300 hPa level.

Ensemble analysis method (EA)
In order to simulate the day-to-day variations in the analysis uncertainty, we have adopted the Ensemble Transform method from ensemble forecasting (Bishop and Toth, 1999;Wei et al., 2008).In ensemble forecasting, the aim is to sample analysis uncertainty in order to simulate their impact on the forecast.In this study, the aim is to simulate the impact from the analysis uncertainties on backward trajectories.One difference between the two purposes is that in ensemble forecasting, the dynamical model can simulate the development of the perturbation properties (e.g. the amplitude relative to the forecast error and spatial correlations) during the first hours of integration which is the case, for example, with singular vectors starting with a low initial perturbation amplitude (Magnusson et al., 2009).Since no forward integration of the model takes place for the trajectory calculations, the above described difference is the rationale for using Ensemble Transform perturbations in the present study since the Ensemble Transform initial perturbations are closer to the expected properties of the analysis error (Magnusson et al., 2009).The Ensemble Transform perturbations are calculated by orthonormalizing the ensemble of perturbations every six hours (for more details on how the perturbations are obtained see Wei et al., 2008;McLay et al., 2008).This procedure is undertaken in order to sample fast growing error structures.The Ensemble Transform method is a further development of the breeding technique (Toth and Kalnay, 1997) and yields perturbations orthonormal to the inverse analysis error norm.The procedure generates perturbation structures that have spatial correlations dependent on the current flow situation and are normalized by the estimated analysis error.The ability for breeding vectors (and also applicable for Ensemble Transform perturbations) to sample the flow dependent error is further discussed in Corazza et al. (2003).
Our data set consists of 20 perturbed analyses, obtained every 6 h.The ensemble of analyses is generated by adding perturbations to the original analysis, the same methodology as in ensemble weather forecasting.By using a simplex transformation (described in e.g.Wei et al., 2008), the mean of all perturbed analyses equal to the unperturbed analysis is obtained.Each perturbed analysis will be a realization of the analysis error.The differences between the perturbed analyses is for instance small changes in the magnitude and direction of the wind and differences in the temperature field (i.e.geopotential).The unperturbed analysis is used to calculate a control trajectory.One additional trajectory is calculated for each perturbed analysis.This means that for each receptor point 21 trajectories are calculated.Note that when referring to the method used to create an ensemble of analyses we use the terminology Ensemble Transform method.The terminology Ensemble analysis method (EA) is used when referring to the trajectory calculations based on this ensemble of analyses.
The amplitude of the analysis perturbations created with the Ensemble Transform method is a measure of the expected uncertainty in the analysis and forecast.If the ensemble standard deviation is large then it is likely that the unperturbed analysis is an unreliable estimate of the true state of the atmosphere, or if it is small it is likely a reliable estimate.However, what is the desired ensemble standard deviation that correspond to the true uncertainties?To verify whether the perturbation amplitude (in the mean) corresponds to the true uncertainties in an analysis or forecast, one can compare the ensemble standard deviation and the RMS error of the ensemble-mean forecast, averaged over time.These two quantities should be equal if the ensemble samples the uncertainties correctly in the mean (see Palmer et al., 2006).In Fig. 1 the standard deviation of the analysis perturbations for the zonal wind as a mean for the two included periods is shown (the meridional wind component behave very similar and is therefore not shown).The ensemble standard deviation is lower in the tropical band compared with the extratropics, mainly due to a lower variability in the wind field over the tropics compared to the mid-latitudes.In order to estimate if the perturbation amplitude corresponds to the true uncertainties for the winter period, we compare the ensemble standard deviation (Fig. 1a, red, solid) and the RMS error of the ensemble-mean forecast (Fig. 1a, red, dashed) after two days.The reason for not comparing the two quantities at the initial time is because it is difficult to obtain independent estimates of the RMS error.The ensemble is somewhat   over-dispersive (i.e. the analysis error is somewhat overestimated) in the northern extra-tropics and under-dispersive (the error is somewhat under-estimated) in the tropics.In the southern extra-tropics the two quantities agree well.In Fig. 1b the corresponding information is shown for the summer period included in this study.We see that the ensemble is still under-dispersive in the tropics and over dispersive in the northern-extra tropics.Independent of season the analysis error appears slightly over-estimated in the northern extra tropics.The lower forecast error in the northern extra-tropics could be explained by the fact the the Northern Hemisphere is well observed in terms of meteorological observations, reducing the uncertainty in the 48-hour forecast.The general characteristics of the analysis error is however similar during both the winter and summer period.We can also see that the initial perturbations (Fig. 1a,b, black lines) are slightly higher in the winter hemisphere of each period, reflecting the more dynamically active winter season.

Initial spread method (IS)
A widely used method to obtain an estimate of the uncertainty in the calculated trajectory path is to slightly change the receptor point of the trajectory.The logical basis behind this method is to see how small uncertainties in the trajectory start point grow with time, thus giving an estimate of the uncertainty in the region of origin.This method to generate a set of trajectories has been used in previous studies to estimate such uncertainties (Merrill et al., 1985;Draxler, 2002).Using this approach we create a set of 26 trajectories with a receptor point displaced by at most ±1 • and 10 hPa relative to the control trajectory start point.This initialization does introduce a latitude dependent bias.We have however accounted for this bias in an extra set of simulations and found no significant difference in the presented result.The trajectories are evenly distributed within the ±1 • box and the distribution is the same for all geographical locations.For reference, we also create a set of trajectories only perturbed in the horizontal plane.This set is not used in the case studies in section 3.However, when considering the statistical behavior of the EA and IS methods we show results from both variations of the IS method.

Case studies
To qualitatively illustrate the evolution of the uncertainties in the trajectory calculations, due to a high analysis uncertainty, we have selected two case studies.We compare trajectories calculated with the EA method and the IS method.For the first case the trajectories are released at 60 • N, 0 • E (North Atlantic) and for the other case at 0 • N, 120 • E (tropics).Five-day backward trajectories are calculated using the two different methods.In total 21 trajectories are calculated for the EA method and 26 trajectories for the IS method for each case study.

North Atlantic case study
The North Atlantic is characterized by low-pressure systems and corresponding areas with strong wind speeds originating from baroclinic instability.When calculating backward trajectories in such regions the location of, for instance, high and low pressure systems could significantly change the calculated trajectory path.Furthermore, due to saddle points in the flow between troughs and ridges, small differences between trajectory locations can change which flow regime they follow.If a trajectory enters the centre of a low-pressure system, fast changes in the vertical level can occur due to the strong vertical wind speeds.
In Fig. 2 the 850 hPa level wind and wind speed standard deviation between the ensemble members (obtained with the Ensemble Transform method) are shown for four subsequent days.In this case we examine trajectories released on the 25th of December, 2005 at 60 • N, 0 • E just north of a highpressure centre located over the British Isles.Two low pressure systems are located over Greenland and Scandinavia.The wind in the region is governed by a relatively weak westerly flow and the location of the receptor point represents a saddle point in the flow.East of the receptor point the wind direction shifts to more northerly winds and to the west the winds are more southerly.The trajectory locations (excluding the first day of the five day period) are shown in Fig. 2 for the EA and IS method, respectively.The region shown in Fig. 2 is characterized by a relatively high standard deviation in wind speed between the ensemble members.This is related to uncertainty in the position of two developing low pressure systems, one over the Atlantic and one over the Norwegian Sea.The main difference found between the EA and IS method is the more irregular behavior of the EA trajectories.For the first day, both methods place the trajectories in the same area.However, after about two days the trajectories enter the circulation of the low-pressure systems.This results in a fast-growing difference between the EA trajectories since the simulated analysis uncertainty determines if the trajectories follow the circulation of the Norwegian Sea low pressure system or the Atlantic low pressure system or if they end up in the high pressure system over western Europe.Thus, the result is three significantly different, but all plausible, regions of origin that are not sampled by the IS method.The more chaotic behavior of the EA trajectories can also be studied in the vertical displacement (not shown) where the IS method trajectories are significantly more coherent in their behavior.The trajectories originating from the same pressure level all follow the same general vertical path.
Both methods sample the error growth due to differences in the atmospheric flow.Once an error in the trajectory location has been introduced, this error will grow due to differences in the atmospheric flow.However, the EA method adds an additional error due to the sampling of the analysis error.This in turn results in a faster error growth.This illustrates that with the IS method one initially mainly sample a slowgrowing error related to the position of the trajectories.The error related to the uncertainty in the atmospheric analysis is not sampled with this method.

Tropical case study
The Ensemble Transform method is mainly aimed at capturing dynamically unstable parts of the analysis error, not observation errors.In tropical regions the analysis error may not be dominated by dynamically growing error structures and, as discussed in Section 2.1, the Ensemble Transform perturbations generally underestimate the analysis error in the tropical region.Figure 3 shows the 850 hPa level wind and wind speed standard deviation between the ensemble members for four subsequent days over the Indo-Asian region.Compared to the North Atlantic case, there are smaller differences between the ensemble members when considering the wind speed.This is due to the smaller initial perturbations in the tropical region.However, if the trajectories were to pass through a region where the analysis uncertainty is increased due to dynamically growing error structures, one would expect differences between the two methods used, similar to the previous case study.Such regions are for instance visible within, and near, the circulation patterns found around 10 • N and 10 • S, which is probably related to equatorial wave activity.In this case we examine trajectories released on the 19th of December, 2005 at 0 • N, 120 • E. The trajectory locations (excluding the first day of the five day period) are shown in Fig. 3.As expected, both methods initially display similar characteristics.However, in Fig. 3b the current flow situation increases the analysis uncertainty at the location of the trajectories.This results in a split between the EA trajectories dividing them into two different source regions, one where the air originated from the east over the Indian Ocean and one where the air originated from the north, south of Taiwan.This is to some extent captured with the IS method as two IS trajectories also indicate the Indian Ocean source region.
Since analysis perturbations using the Ensemble Transform method are generally small in the tropics (cf.using other methods, aimed at capturing the non-growing part of the analysis error (i.e. the effects of observational errors and the lack of observations) could improve the sampling of the analysis uncertainties in tropical regions.Additional errors related to convective transport or turbulence, which is not included in the trajectory model, could also add to the trajectory uncertainty in this region.However, the EA method does show that the two suggested source regions are equally plausible and that the uncertainty in the calculated area of origin is higher than using the IS method, albeit not as high as in the mid-latitudes.

Statistics
Even if the Ensemble Transform method is mainly aimed at sampling flow-situations, where high analysis uncertainty results from dynamically growing error structures, statistics of the main behavior of the trajectory spread is of interest.In order to obtain statistics of the deviation between the

30N
EA method IS method IS method (only horizontal spread) Fig. 4. Trajectory mean deviation for the EA and the IS method as a function of the trajectory length for 60 • N, 30 • N, 0 • N, 30 • S and 60 • S. Calculated as the difference between one perturbed trajectory (either using the EA method or the IS method) and the unperturbed (control) trajectory.The mean deviation is calculated as the mean of both the summer and winter period.24 Fig. 4. Trajectory mean deviation for the EA and the IS method as a function of the trajectory length for 60 • N, 30 • N, 0 • N, 30 • S and 60 • S. Calculated as the difference between one perturbed trajectory (either using the EA method or the IS method) and the unperturbed (control) trajectory.The mean deviation is calculated as the mean of both the summer and winter period.trajectories, the mean deviation for the summer and winter period trajectories released at 850 hPa is calculated.Statistics are obtained for 5 different latitude bands (60 • N, 30 • N, 0 • N, 30 • S, 60 • S) expressed as the mean for the longitudes 120 • S, 0 • N and 120 • N. The deviation is calculated as the spherical distance at each time step between one perturbed trajectory (either using the EA method or the IS method) and the unperturbed (control) trajectory.The mean deviation as a function of trajectory time is plotted for 60 • N, 30 • N, 0 • N, 30 • S and 60 • S in Fig. 4. We see that the EA method start from zero deviation, as expected as there is no perturbation in the start point using the EA method.This is different to the IS method, which has an initial spread by definition.However, we see that for all latitudes, the EA method obtains a larger deviation with time due to a faster growth of the spread.This is most prominent for 60 • N and 60 • S, where the mean deviation is increased with 40% after five days.We also see that the IS cluster perturbed in both the horizontal and vertical direction obtains a larger deviation compared with the set of trajectories only perturbed in the horizontal plane.The shape of the evolution of the trajectory deviation appears to be exponential for the IS method.This behavior could be related to chaotic advection which is discussed in e.g.Bofetta et al. (2000).We propose a simple model of the deviation growth which is set up as where D is the deviation and αD represents the dispersion due to chaotic advection which is related to the Lagrangian Lyapunov exponent of the system.The parameter β represents the effect of the uncertainty in the analysis and is also a function of the distance.To verify this model the deviation speed is plotted as a function of deviation (Fig. 5).Since one deviation speed is obtained for each mid time step for each trajectory, a large amount of data points are obtained.Therefore the data is grouped in bins populated with 5000 data points.For short deviations, the deviation speed is higher for the EA trajectories, due to the effect of the perturbed analysis.One can view this as a displacement by β in Eq. ( 1).
The difference is approximately 0.5 ms −1 (cf.Table 1) and is of the same order of magnitude as the analysis perturbations (see Fig. 1).As the deviation grows, the relative effect of the EA method decreases and all methods have approximately the same deviation speed for deviations larger than 400 km.For the IS method, the deviation speed increases linearly, where the slope of the curve can be estimated as α, with the deviation distance for deviations up to 400 km.This indicates that the model is a reasonable approximation of the deviation growth.For larger deviations (>1000 km) the deviation speed saturates.This property of the deviation growth is discussed in Bofetta et al. (2000), and it can be viewed as two different regimes.One early regime with exponential deviation growth (for deviations shorter than a certain length scale) and normal diffusion for larger deviations.The length scale is dependent on the distance for which the wind field is correlated.A typical approximation of this length scale would be the Rossby radius of deformation.If the trajectories are located in different cyclone systems and separated by a distance larger than the Rossby radius, the error growth will be independent of the deviation.
In Table 1 the parameters α, obtained by linear regression, and β at D = 0 is presented for the winter and summer month separately, as well as for the trajectories released at 300 hPa.We can see a small seasonal dependence in β with slightly higher values in the winter hemispheres.This is probably related to the seasonal change in the position of the storm tracks and the higher uncertainty in the atmospheric analysis with more dynamic instability.This is also reflected in Fig. 1 where the initial perturbations are higher in the winter hemispheres.At 300 hPa the effect of β is higher compared to the 850 hPa level.However, in addition to a higher analysis uncertainty, the wind speeds at 300 hPa is also higher.The resulting relative difference in mean deviation between the EA and IS trajectories is therefore similar to the lower altitude cases.Note also that the highest β is found for the midlatitudes during both the summer and winter period.This is probably related to that the Ensemble Transform method generally underestimate the analysis error between 30 • N and 30 • S. The seasonal dependence can also be seen, to some extent, in α with the higher values found in the winter hemisphere of each period.This indicates a stronger influence from chaotic advection in the winter hemispheres while the advection is more laminar in the summer hemispheres and in the tropics.

Conclusions
The uncertainty in trajectory calculations due to errors in the atmospheric analysis is estimated by simulating uncertainties in the atmospheric flow.Errors in the analysis are unavoidable and these errors vary from day to day and is dependent on the number of available atmospheric observations and on the current flow situation.Two different methods to generate a set of trajectories have been compared.These are the Initial spread method (IS) and the Ensemble analysis method (EA).
To simulate analysis errors, the Ensemble Transform method from ensemble forecasting has been adopted.In ensemble forecasting the aim is to simulate the effect of analysis errors on forecast errors.In this study we simulate the effect on backward trajectories.By perturbing the analysis with perturbations of the same order of magnitude as the analysis uncertainties, we obtain (in our case) 20 equally plausible atmospheric states that are used for the trajectory calculations using the EA method.Each backward trajectory is using each one of these new states for each one of the locations.Every 6 hours a new set of perturbed analyses is obtained and used in the calculations.
The effect of performing trajectory calculations while sampling the uncertainty in the atmospheric analysis has been compared with the effect of running a set of trajectories, each with a perturbed receptor point (IS method).By studying a number of cases, we find that a large difference between the two methods appears under certain atmospheric situations.This occurs when the trajectories pass an area where the current flow situation results in analysis uncertainties.This is especially present in a saddle point of the flow, where small changes in the flow situation could lead to large differences in the trajectory paths.
We have also compared mean statistics calculated for 25 successive days during December 2005 and for 29 successive days during July 2005 for five latitudes, 60 • S, 30 • S, 0 • N, 30 • N and 60 • N, expressed as the average of 120 W, 0 • E and 120 • E. For a short trajectory time length we see that the deviation of the set of trajectories is larger due to the perturbations in the initial point (IS method).The set of trajectories using different analyses (EA method) start at the same point and therefore the initial spread is zero.But after a certain time, the mean deviation becomes larger for the EA method.To investigate the dynamics of the deviation growth in more detail, the growth rate has been studied as a function of the deviation distance.We found a linear dependence between the growth rate of the deviation and the deviation itself for small deviations between trajectories for the IS method.This is expected due to chaotic advection and implies an exponential growth of the deviation with time.For the EA method, the growth rate is much higher for small deviations compared to the IS method.This effect is due to the perturbation added to the analysis.When the deviation becomes large, the effect of the differences in the flow between the trajectory locations is more important than the analysis perturbations, and the growth rate becomes equal for both methods.
The fact that both methods have the same asymptotic properties in the error growth rate suggest that it should be possible to find an amplitude of the initial point disturbances yielding the same dispersion for long trajectories as the EA method.But that procedure has several disadvantages.Firstly, one overestimates the deviation for short trajectories.Secondly, but not less important, one will only sample uncertainties in the mean and may not capture cases where the analysis uncertainty plays a dominant role.
In this study we have used one analysis perturbation method and we cannot conclude that this method is superior to any of the other available methods.For this more research is needed.For example, the Ensemble Transform perturbations are designed to mainly catch the dynamically growing part of the analysis uncertainties, which is of interest for ensemble forecasting.But for the application of trajectory calculations, also the non-growing part is of interest.How to include this part needs further studies.The Ensemble Transform method does not include direct effects from observation uncertainties.Therefore one could expect a more reliable sampling of the uncertainties from analyses obtained from an ensemble data assimilation system using perturbed observations (Houtekamer and Mitchell, 1998;Buizza et al., 2008).Furthermore, we have only studied trajectories in the troposphere, where the analysis error is dominated by dynamically growing error structures.Trajectories calculated at higher altitudes could behave differently if the analysis error is not dominating the uncertainty in the trajectories.The stratosphere, for instance, is less chaotic than the troposphere which would decrease the error growth, but the number of observations in the stratosphere is at the same time fewer.For the same reasons as discussed above, a different sampling of the uncertainties in the stratosphere could help to better understand the effect of the analysis uncertainty on trajectories at these altitudes.Ensemble data assimilation, which add stochastic components to observations, could therefore also be a viable option in future studies of trajectories in the stratosphere.
The overall error in trajectory calculations due to uncertainties atmospheric motion will depend on many different factors and we have here studied one of these.Adding e.g. a chaotic component related to turbulence or convection would probably introduce an additional uncertainty, especially in the tropics.However, these errors act on a different scale than the analysis error and depend on parameterizations within the trajectory model.How this would impact the error growth speed needs to be studied further before any conclusions can be drawn.One could also imagine an ensemble with different variations of these parameterizations to estimate the uncertainty in a trajectory.How to construct such an ensemble would also need to be studied further.
The main conclusion from this study is that by perturbing the analysis consistent with the analysis uncertainties, both regarding perturbation amplitude and correlation length, the uncertainties in trajectory calculations can be more consistently estimated.Therefore we suggest that a set of analysis perturbations should be constructed for the reanalysis data sets used for dispersion calculations to obtain reliable estimates of transport trajectory uncertainty.This is especially necessary when performing case studies where the current flow situation could significantly increase the uncertainty in the trajectory calculations.

Fig. 1 .Fig. 1 .
Fig. 1.Initial ensemble standard deviation (black line), ensemble standard deviation hours (red, solid) and the RMS error of the ensemble-mean forecast at 48 hours (red, for the zonal wind speed during a) the northern hemisphere winter period and b) the n hemisphere summer period.

Fig. 2 .
Fig. 2. 12Z wind at 850 hPa (arrows) and the wind speed standard deviation in ms −1 (gray shading) between the ensemble members for four subsequent days included in the North Atlantic case study.a) 24th, b) 23th, c) 22th and d) 21th of December, 2005.Markers show the corresponding location of the EA method trajectories (red) and the IS method trajectories (green). 22

Fig. 2 .
Fig. 2. 12Z wind at 850 hPa (arrows) and the wind speed standard deviation in ms −1 (gray shading) between the ensemble members for four subsequent days included in the North Atlantic case study.(a) 24th, (b) 23th, (c) 22th and (d) 21th of December, 2005.Markers show the corresponding location of the EA method trajectories (red) and the IS method trajectories (green).

Fig. 3 .
Fig. 3. 12Z wind at 850 hPa (arrows) and the wind speed standard deviation in ms −1 (gray shading) between the ensemble members for four subsequent days included in the Tropical case study.a) 18, b) 17th, c) 16th and d) 15th of December, 2005.Markers show the corresponding location of the EA method trajectories (red) and the IS method trajectories (green). 23

Fig. 3 .
Fig. 3. 12Z wind at 850 hPa (arrows) and the wind speed standard deviation in ms −1 (gray shading) between the ensemble members for four subsequent days included in the Tropical case study.(a) 18, (b) 17th, (c) 16th and (d) 15th of December, 2005.Markers show the corresponding location of the EA method trajectories (red) and the IS method trajectories (green).

Fig. 5 .
Fig. 5. Trajectory mean deviation speed for the EA and the IS method as a function of the mean deviation for 60 • N, 30 • N, 0 • N, 30 • S and 60 • S. Calculated as the time derivative of the deviation shown in Figure 4.The mean deviation speed is calculated as the mean of both the summer and winter period.

Fig. 5 .
Fig. 5. Trajectory mean deviation speed for the EA and the IS method as a function of the mean deviation for 60 • N, 30 • N, 0 • N, 30 • S and 60 • S. Calculated as the time derivative of the deviation shown in Fig. 4. The mean deviation speed is calculated as the mean of both the summer and winter period.