Time-Reversibility in Atmospheric Dispersion

Due to the chaotic nature of atmospheric dispersion, small deviations, e.g., numerical errors in dispersion simulations, increase rapidly over time. Therefore, the accuracy of backward simulations is limited. In the paper, the degree of the fulfillment of time-reversibility over different time periods is investigated by a Lagrangian dispersion model at a global scale using pollutant clouds consisting of a large number of particles. The characteristics of the pollutant clouds in the backward simulation are compared to those in the forward simulation. In order to characterize the degree of time-reversibility, Lagrangian quantities, such as the fraction of particles that return to the initial volume, the center of mass and the standard deviation of the pollutant clouds, are determined. Furthermore, the overlap and the Pearson’s correlation coefficient between the forward and backward clouds are also investigated. Both a case study and global results are presented. Simulations reveal that the accuracy of time-reversibility decreases in general exponentially in time. We find that after the Lyapunov time of the dispersion (in our case three to four days), the results of the backward tracking become unreliable, and any sign of time-reversibility is lost.


Introduction
Backward calculation of the transport and dispersion in the atmosphere is a tool commonly used for identifying the sources of pollutants or other substances and/or to provide a dispersion model with input data optimized to all of the knowledge exploited from available measurements for future studies.Besides simple trajectory computation (see, e.g., [1]), more sophisticated techniques exist to estimate the location and distribution of substances at a certain time before the observation.These methods are described in detail in [2] and the references therein, including inverse modeling using a source-receptor relationship by means of a Lagrangian model [3] and adjoint versions of Eulerian dispersion models [4,5].
In this study, we focus on the degree at which the time-reversibility property is fulfilled in backward tracking of individual pollutant particles.For the investigation, we use a Lagrangian trajectory model applied in many research areas in meteorology, including calculation of the origin of particulate matter [6,7], pollen [8,9], moisture [10,11] or gases and passive tracers [12][13][14].In the simplest case, when only advection is taken into account, particle trajectories are determined utilizing 3D velocity fields of the atmosphere, and the motion of the particles is deterministic and described by three time-dependent differential equations in the horizontal and vertical directions.
In deterministic dynamical systems, the naive expectation suggests that if the system develops for a certain time forward, then time is reversed, and if the evolution of the system is followed backwards for the same time period, the initial state should be recovered.This is expected to hold in simple dynamical systems rather accurately, even if the development of the system is determined numerically.
However, if the deterministic dynamical system is nonlinear and is described by at least three first-order autonomous differential equations (or by two time-dependent differential equations), chaotic behavior can appear [15][16][17], the typical characteristics of which are sensitivity to initial conditions, irregular motion and the appearance of complex formations, like folded, elongated filaments and tendrils [18].Three-dimensional passive advection has three degrees of freedom, that is advection dynamics can result in chaotic motion even if the flow field by which the particles are advected is smooth and deterministic.The characteristics of chaos were also experienced in atmospheric dispersion problems of contaminants (see, e.g., [19][20][21]).
Chaotic behavior is the reason for the fact that the motion of a particle is unpredictable for long times as the trajectories of initially nearby particles diverge in an exponential manner.Small differences in the initial position and numerical approximations used in the calculation induce errors growing exponentially in time.In view of the sensitivity to slight deviations, instead of one trajectory, several adjacent particles should be tracked to get a more precise overview of the volume/area covered potentially by the contaminant.
In the dynamical systems literature, the breakdown of time-reversibility in chaotic systems due to the unavoidable numerical errors was first pointed out by Shepelyansky [22].An illustrative example with an ensemble of particles in an elementary problem is given in the textbook [17].
In this paper, we investigate the accuracy of backward trajectory simulations in atmospheric dispersion as a function of time and, thereby, the degree of the fulfillment of time-reversibility in the view of chaotic advection.As a first step, a forward trajectory simulation of a pollutant cloud consisting of several particles emitted instantaneously is carried out for a maximum of nine days using a Lagrangian trajectory model.Figure 1 provides a first insight into the basic phenomenon.The black dots in Figure 1b trace out the position of an initially localized particle cloud for time instant t = 7 days.Particles correspond to passive tracers, that is they represent solid or liquid contaminants of negligible size or a gaseous contaminant.The location of the pollutant clouds of the forward run started at day t = 0 serves as a reference for the succeeding backward simulations.Backward simulations start at t b = 1, 2, . . ., 9 days from the particle distribution of the forward simulation at t = t b as the initial condition and track back the particles till t = 0 (see the colored dots in Figure 1a).For clarity, an outline of the time schedule of the forward and backward simulations is given in Figure 2. In general, the accuracy of backward tracking as a function of the time interval ∆t b = t b − t of the backward simulation is assessed by comparing the forward cloud and the backward clouds at time instants t = 0, 1, . . ., 8 days.For example, we compare the forward cloud with nine backward clouds at t = 0, three of which are illustrated in Figure 1a, and with two backward clouds at t = 7 days, only one of which is represented in Figure 1b for better visibility, respectively.In this way, the "exact" initial conditions of the backward simulations at time t b are known within the rounding-off error of the computer.The forward pollutant cloud at time t is considered to be the "original" position of the particles where the particles should return.The deviations between the backward cloud and forward cloud are caused by the rounding-off errors and errors of the numerical scheme for integration of the equation of motion.These errors are amplified over time due to chaoticity. Figure 1 shows that the longer the time interval of the backward tracking, the worse the coincidence between the forward and backward clouds, e.g., the bad agreement of the forward and backward cloud can be illustrated by the red dots with ∆t b = 9 days of backward tracking in Figure 1a spreading from Asia to Europe, while the blue dots in Figure 1a and the red dots in Figure 1b both with ∆t b = 2 days of backward tracking are in reasonable agreement with the corresponding forward clouds (denoted by black dots), respectively.In order to analyze the matching of the appropriate backward and forward clouds, statistical parameters are determined.We investigate Lagrangian quantities, which depend solely on the position of the particles, and Eulerian ones, which require the particles to be projected into grid cells.In this way, the proportion of the particles returned back to the initial volume of the cloud is determined, as well as some parameters already applied in earlier studies [23,24].These quantities are the difference between the center of mass of the backward and forward pollutant clouds, the difference between the standard deviations of the particles in the clouds and the overlap and Pearson's correlation coefficient of the two clouds.The paper presents a case study and, furthermore, simulations distributed uniformly over the globe to get a general overview of the behavior of the statistical parameters.As a general rule, parameters reveal that the accuracy of backward tracking decreases exponentially in time, and usually after a few days or about a week, the results of the backward tracking become unreliable; any sign of time-reversibility is lost.
The paper is organized as follows.Section 2 provides a brief overview of the RePLaT model by which the particle trajectories are monitored.Section 3 describes the meteorological data used for the simulations and the initial distribution of the pollutants, while Section 4 introduces the statistical parameters applied for the characterization of the accuracy of backward simulations.The results of a case study and global studies are presented in Sections 5 and 6, respectively.An outlook is given in Section 7, and Section 8 summarizes the conclusions of the work.

The RePLaT Dispersion Model
The RePLaT (Real Particle Lagrangian Trajectory) dispersion model, first formulated in [25], is used to calculate the trajectories of particles that compose the pollutant clouds.RePLaT is a Lagrangian model tracking aerosol particles with a realistic radius and density.Ideal tracers and gaseous contaminants correspond to particles with radius r = 0. RePLaT is applied here in this limiting case.Although there exist a few Lagrangian trajectory models where the advection of a particle is determined by higher order integration schemes (e.g., TRAJ3D [26]), in most cases, this process is computed by integrating the equation of motion by using Euler's method (NAME [27], MLDP0 [28], SNAP [29], GEARN [30]) or second-order methods (FLEXPART [31], HYSPLIT [32]).RePLaT also calculates particle trajectories by means of the forward Euler's method: where r p (t) indicates the particle position at time t, v(r p (t), t) is the wind speed vector at the location of the particle and ∆t is the time step.In forward simulations, ∆t > 0, while in backward simulations, ∆t < 0. Since after a few days a particle cloud traces out the unstable foliation of a chaotic dynamics, the shape of this pattern is expected to be insensitive to the choice of the numerical code.The reason is that exponential divergence occurs along the unstable foliation only; in the perpendicular directions exponential contraction occurs, just like for chaotic attractors [17].Any numerical algorithm is thus able to easily converge to the unstable foliation.Where unpredictability plays a role is the position of a given trajectory endpoint along the unstable foliation; this, however, is not an issue in our simulations, since instead of individual trajectories, an ensemble is treated.RePLaT also reckons with the impact of turbulent diffusion on the particles as a random walk process, and the use of Euler's method is motivated by an easy treatment of this part of the code.However, since in large-scale simulations, its effect is much weaker than that of advection [19], as a first approximation, turbulent diffusion is neglected in the simulations.Turbulent diffusion will only be treated in Section 7.
In the present study, if an ideal tracer particle encounters the surface, it bounces back to the atmosphere with a perfectly elastic collision.
In order to calculate the trajectories of the particles using Equation (1), meteorological data given on a regular grid are interpolated to the position of the particles (using bicubic spline interpolation in the horizontal direction and linear interpolation in the vertical direction and in time).For more details about RePLaT, see [25].

Input Data
The particle trajectories are computed using the v = (u, v, ω) reanalysis fields of the ERA-Interim database of the European Centre for Medium-Range Weather Forecasts (ECMWF) [33], where u and v are the longitudinal and latitudinal velocity component of air, respectively, and ω is the vertical velocity in pressure coordinates.The meteorological variables utilized in the simulations are given on 32 pressure levels between 1000 and 100 hPa on a 1.5 • × 1.5 • horizontal grid with 6 h of time resolution.
The time step of the simulation is determined by the inequality ∆t < min ∆x u , ∆y v , ∆p ω , where ∆x, ∆y and ∆p are the resolution of the grid of the meteorological data [34].Based on this relation, the time step is chosen to be ∆t = 5 min.The time interval of particle tracking in both the case study and the global simulations is at most nine days beginning on 12 March 2011 at 00 UTC.The particles in the forward simulations are emitted instantaneously in the atmosphere at this time instant.
In the case study (Section 5), n 0 = 1.25 × 10 5 particles are initially distributed uniformly in a 1 • × 1 • × 100 hPa volume centered at λ 0 = 141 • E, ϕ 0 = 37.5 • N, p 0 = 500 hPa.It can be considered as a hypothetical sudden emission over Fukushima occurred on 12 March 2011 at 00 UTC, however, well above the real emission.We chose the time period because of the well-documented meteorological situation (see, e.g., [35]), but particles are placed higher in the atmosphere in order to be able to investigate the transport on large scales.
In the global simulation (Section 6), the initial position of the pollutant clouds is of the same size at the same altitude as in the case study.From λ 0 = 150 • W to 180 • E in 30 • increments and ϕ 0 = 80 • S to 80 • N in 20 • increments, 12 × 9 pollutant clouds are emitted instantaneously at the time mentioned above on a "grid" in the atmosphere.The results of Figure 1 arise from an instantaneous injection in one of these grid points.

Statistical Parameters
In order to quantify the accuracy of the backward simulation, the pollutant clouds of the backward and forward simulations at the appropriate time instants are compared using different statistical quantities, which are detailed in the subsequent subsections.The proportion of the particles n/n 0 that returned to the initial volume, the deviation of the center of mass and the difference of the standard deviation of the pollutant clouds are solely based on the position of the particles, while the figure of merit in space (overlap) and the Pearson's correlation coefficient (suggested by, e.g., [24]) demand the use of an artificial grid over which concentration can be calculated.
The differences between the backward and forward pollutant clouds (BWC and FWC) are considered as a function of the time interval of the backward simulation ∆t b = t b − t, where t b is the initialization time of the backward simulation (t b = 1, 2, . . ., 9 days) and t is the time instant of the comparison (t = 0, 1, . . ., 8 days).Since we have backward simulations beginning at t b = 1, 2, . . ., 9 days, e.g., for ∆t b = 1 day, nine values are available for each measure (except for n/n 0 ), because the BWC with t b = 1 day can be compared to the FWC at t = 0, the one with t b = 2 days at t = 1 day, and so on.Therefore, for ∆t b = 9 days, we have only one value, since only the backward cloud with t b = 9 days can be compared to the forward cloud at t = 0.

Proportion of Particles Returned to the Initial Volume
For each t b , the number n of the particles that returned to the initial 1 • × 1 • × 100 hPa volume is determined.The proportion n/n 0 describes how accurately a backward simulation of ∆t b = t b can track back the pollutant cloud to the source.

Center of Mass and Standard Deviation of the Pollutant Clouds
The horizontal coordinates of the center of mass CM h,f and CM h,b of the FWC and BWCs and the standard deviation σ h,f and σ h,b of the particles in the clouds are calculated at t = 0, . . ., 8 days, and the deviation between the backward and forward simulations is determined.In order to eliminate the effect of the convergence of meridians, CM h,f and CM h,b are not computed by averaging the longitudinal and latitudinal coordinates of the particles, but by converting their position to 3D Cartesian coordinates before averaging ((0, 0, 0) being the center of the Earth).The averaged coordinates are projected back to the surface of the Earth in latitude-longitude coordinates.The standard deviation of the particles in the clouds and the differences between the center of masses of the BWCs and FWC are calculated using spherical distances: • 180 π 111.1 km (3) where λ i , ϕ i are the latitudinal and longitudinal coordinates of the i-th particle in the cloud studied and CM h,f,λ , CM h,b,λ and CM h,f,ϕ , CM h,f,ϕ are the same for the center of mass of the FWC and BWC, respectively.The spherical distance between two locations given by latitudinal and longitudinal coordinates is traditionally given in radians [36].In order to represent it in kilometers, radians are changed to degrees (multiplication by 180 • /π), and then, they are multiplied by 111.1 km/1 • , the average length of one degree on the surface of the globe.
The difference ∆CM h between the center of masses characterizes the average shift between the FWC and BWC, while the difference ∆σ h = σ h,b − σ h,f between the standard deviations describes the typical difference in the extension of the BWC and FWC.

Figure of Merit in Space
The percentage of overlap (called the figure of merit in space) between the BWC and FWC is also calculated.It is defined as the ratio between the intersection and the union of the two pollutant clouds [24].Besides the conventional 2D overlap FMS 2D (i.e., overlap of air columns), the 3D overlap FMS 3D is also determined in the study: where A b and V b are the area and volume of the backward pollutant clouds, respectively, and A f and V f correspond to those of the forward clouds.If A b = A f , there is an exact overlap between the two clouds, that is FMS 2D = 100%.When FMS 2D < 100%, the BWC and FWC do not cover each other completely due to their different shape and/or shift following from their different location.
As RePLaT tracks individual particles, in order to estimate A b , A f and V b , V f particles need to be projected onto a grid.In the horizontal direction, instead of the obvious regular longitudinal-latitudinal grid, we use a modified grid detailed in [19] in order to eliminate a strange feature of the regular grid, namely that a cell of ∆λ × ∆ϕ has a much smaller area close to the poles than near the Equator.In the modified grid, the latitudinal side ∆ϕ of a cell is equal at any latitude, and the longitudinal size ∆λ changes depending on the latitude, so that the area of the cells is almost equal.In this paper, A b , A f and V b , V f denote the number of such cells with at least one particle.Note that in the value of FMS 2D and FMS 3D , rarely occupied cells are over-weighted, since a cell consisting of only one particle is taken into account with exactly the same weight as a cell with a large number of particles.

Pearson's Correlation Coefficient
Pearson's correlation coefficient [24] is also determined for air columns (PCC 2D ) and for air volumes (PCC 3D ) to study whether the backward and forward clouds are correlated in space.The PCC 2D/3D is defined as: where K is the number of 2D or 3D cells where the backward and/or forward clouds are present (here, cells to which at least one particle belongs), c b/f (i) is the number concentration (the number of particles in the cell divided by n 0 ) of the backward/forward cloud in cell i and c b/f is the mean concentration of the backward/forward cloud.

Case Study
In order to study the accuracy of time-reversibility in the dispersion of pollutant clouds with deterministic Equation ( 1), as a first exercise, a forward simulation of nine days and the corresponding backward simulations are carried out.The particles of the FWC are initiated in a volume of 1 • × 1 • × 100 hPa at the time and location mentioned in Section 3. Figure 3 illustrates the FWC (black) and the BWCs (colored) at t = 0 and 2, 5 and 7 days after the emission.The meteorological situation during the studied days is detailed, e.g., in [35].The black FWC leaves Japan eastwards with the strong flow of a jet, while due to the shearing in the wind field, it stretches more and more.On 15 March (not shown in Figure 3), the elongated filament gets into a flow region dominated by a cyclone at the western coast of North America due to which it begins to form a spiral.Some days later, as was also experienced in the case of the real emission from Fukushima, some fraction of the FWC reaches Europe.The complex, folded and rather stretched shape of the initially compact pollutant cloud is the consequence of the chaotic nature of the advection mentioned in the Introduction.
A fairly good agreement and almost total overlap appear on 12 March between the FWC and the dark blue BWC (t b = 1 day, ∆t b = 1 day), on 14 March for the cyan and light blue BWC (t b = 3, 4 days, ∆t b = 1, 2 days), on 17 March for the green and yellow BWC (t b = 6, 7 days, ∆t b = 1, 2 days) and on 19 March for the orange BWC (t b = 8 days, ∆t b = 1 day).Additionally, on 12 March, the BWC with t b = 2 days (∆t b = 2 days), on 14 March, the BWC with t b = 5 days, on 17 March the BWC with t b = 8 days and on 19 March, the BWC with t b = 9 days (∆t b = 3 days in these cases) do not differ too much from the corresponding FWC either.The observation suggests that for ∆t b = 1-3 days, the backward simulations provide reasonable agreement, with the FWCs considered here as the "measured" position of the pollutants.For larger values of ∆t b , however, the agreement is poor, and a clear breakdown of time-reversibility is found.
The agreement of the FWC and BWCs can be also studied not only by the naked eye, but in a quantitative way.In Figure 4, the example of the horizontal standard deviation σ h,f and σ h,b of the clouds are shown.Figure 4a reveals that after a few days of exponential growth (the exponent is Λ ≈ 0.8 day −1 ), σ h,f exhibits a power-law behavior with an exponent of γ = 0.75.The crossover between an initially exponential growth to a power-law is a general observation of studies of chaotic advection [37].The crossover is set by the largest characteristic size of the flow, beyond which particles carry out a random walk in-between these large-scale structures.In our case, this size is that of the cyclones, and indeed, σ h,f is on the order of 1000 km when the crossover takes place.The value of exponent γ is also in harmony with previous studies; see, e.g., [38].Huber et al. experienced superdiffusivity (roughly ballistic dispersion) in mid-and high latitudes on time-scales of 10 days and found that in some mid-latitudinal locations, an exponential growth for short time intervals was present.Figure 4b shows the dependence of σ h,b on the time interval of the backward simulation ∆t b = t b − t.That is, the data in Figure 4a are projected to the vertical axis of Figure 4b, because they are the initial conditions of the backward simulations with t b = t.Dashed lines connect the values corresponding to the time instants t when the properties of the FWC and BWCs are going to be compared.One can see that in the first 3-4 days of the backward simulations, σ h,b decreases in a similar manner as σ h,f increases in the forward mode.The almost horizontal nature of the dashed lines show that the values of σ h,b are close to σ h,f .Later, however, the decrease slows down, and within 5-6 days, it switches to an increase again.This increase is a clear sign of irreversibility.To quantify the differences between the BWCs and FWCs, the statistical parameters introduced in Section 4 are calculated.Figure 5 illustrates the Lagrangian ones when the calculation is based only on the position of the particles.
The horizontal distance ∆CM h between the center of mass of the FWC and the BWCs increases as a function of ∆t b .The chaotic nature of the advection in the atmosphere implies that initially close particle pairs move away from each other exponentially.Therefore, the naive expectation emerges that the difference between the center of mass of the FWC and the BWCs can be described by an exponential function, as well.Figure 5a confirms this behavior: the mean data indeed seem to fit well to an exponential function ∼ exp(Λ∆t b ) − 1 with Λ = 0.481 day −1 .The reason of the presence of term −1 is that in this way, the curve has to go through (0, 0), as is expected for ∆t b = 0.The BWCs with t b of warm colors (larger t b ) usually perform worse than the others.This is the consequence of the fact that, e.g., for t b = 9 days and ∆t b = 4 days (i.e., the BWC at t = 5 days), the extensions of the BWC and the FWC are larger than that for t b = 6 or 7 days with the same ∆t b (t = 2, 3 days), and the meteorological impacts on the pollutant clouds distributed on a much wider geographical area influence the position of the center of masses more strongly.From ∆t b = 1 to about five days, a widening distribution of the ∆CM h values of the BWCs can be seen due to the longer backward simulation.The observed distribution happens to narrow for ∆t b > 5 days, since it consists of less and less data.

Figure 5b illustrates the difference ∆σ
between the horizontal standard deviations of the BWCs and FWC.The panel reveals that in most cases, the BWCs cover a larger region than the corresponding FWC (∆σ h > 0).This is also a sign of the chaotic advection of the pollutant particles: small deviations grow rapidly over time (in this case, as a function of the time ∆t b of backward tracking), and therefore, particles spread to a larger region than where they were at the same time instant in the FWC (see, e.g., t b = 9 days (red) in Figure 3).For ∆t b = 1-4 days, an exponential growth with Λ = 0.376 day −1 can be identified.Based on the characteristics of the growth in Figure 4a, over larger ∆t b values, a power-law function is fitted to the mean data.However, note that we cannot expect either a clear exponential or a clear power-law behavior, since exponential and power-law regimes of σ h,f and σ h,b do not extend through the entire nine-day interval of investigation in ∆t b .The power-law fit turns out to have an exponent of about γ = 4.3.This γ is much larger than the γ for σ h,f (t), which implies that an exponential fit to ∆σ h over the entire interval of interest might also be appropriate.The exponential function ∼exp(Λ∆t b ) − 1 fitted to the mean data has a similar exponent as the one in the case of ∆CM h (see Table 1).
Table 1.The exponent Λ (day −1 ) and the error of the exponential fits to the curves ∆CM h , ∆σ h , n/n 0 , FMS 2D , FMS 3D , PCC 2D and PCC 3D versus ∆t b of the case study in Figures 5 and 6.
-Λ ∆CM h 0.481 ± 0.032 ∆σ h 0.467 ± 0.050 n/n 0 0.339 ± 0.032 FMS 2D 0.466 ± 0.048 FMS 3D 0.456 ± 0.040 PCC 2D 0.447 ± 0.028 PCC 3D 0.437 ± 0.027 The proportion of the particles that are able to return to the inside of the initial cuboid (Figure 5c) is n/n 0 = 0.8-0.9 for ∆t b = 1-2 days, in agreement with Figure 3a; then, it falls off quickly to 0.1 for ∆t b = 8-9 days.The decrease can again be approximated by an exponential decay: n/n 0 ∼ exp(−Λ ∆t b ), where Λ = 0.339 days −1 .Although, the data are somewhat scattered around the exponential behavior, the type of the fit will be confirmed in the next section.
Figure 6 represents the 3D figure of merit in space FMS 3D and Pearson's correlation coefficient PCC 3D .The resolution of the grid used for the calculation is about 20 km × 20 km × 20 hPa (∆λ = ∆ϕ = 0.2 • at the Equator).The alteration of the grid resolution, e.g., to ∆λ = ∆ϕ = 0.5 • , ∆p = 50 hPa, does not change the whole picture significantly, i.e., the trends and the fitted exponents are the same.Obviously, too small cells (compared to the number of released particles n 0 ) can result in "split-to-particle" clouds that have a lower chance to overlap, even if their visual expression suggests that.The limiting case is when the size of the cells tends to zero, and no overlap between the appropriate clouds can be found, since infinitely small "cells" (i.e., the individual particles) have a negligible chance to coincide.The 2D and 3D figures of merit and Pearson's coefficients behave similarly with a shift of the measured data towards smaller values for 3D.The reason for the shift lies in that not only the horizontal, but also the vertical distribution of the clouds has to be similar.As expected, the longer the time period ∆t b of the backward simulation, the worse the overlap and correlation of the BWCs and FWC.After ∆t b ≥ 5 days of backward tracking, the BWC becomes so much extended compared to the original FWC at that time that there is less than a 10% overlap between the clouds.The shift of the clouds (which is quantified in Figure 5a as ∆CM h ) also contributes to the lack of the coincidence of the two pollutant clouds investigated.
Although, exponential decay (∼exp(−Λ∆t b )) of the FMS 3D and PCC 3D cannot be explained directly by particle pair separation, since not only the shift and the change in extension of the clouds, but also the shape of them contribute to these quantities, Table 1 shows a reasonable agreement between the Λ of FMS 2D , FMS 3D , PCC 2D , PCC 3D and ∆CM h , ∆σ h .
Note that the exponent fitted to n(∆t b )/n 0 is 0.339, which is different from those of the other quantities.We can discover here the signature of the so-called Lyapunov exponent [16] by measuring the horizontal distance between the initial position of a particle and the position of the returned particle of the backward simulation at t = 0 for all of the particles.The mean distance exhibits again an exponential growth as a function of ∆t b (∼exp(Λ2∆t b ) where 2∆t b stands for the total time period over which the particles are advected).The Lyapunov exponent obtained this way is found to be Λ = 0.364 ± 0.056 days −1 , which coincides reasonably well with the one obtained for n/n 0 .This Λ is in reasonable agreement with the atmospheric Lyapunov exponents of advection found in earlier studies (e.g., [21,39]) for a similar time span.The Λ exponents of the other statistical parameters can also be considered to be Lyapunov exponents; however, in those cases, meteorological situations of different days all contribute, the time spans are different and the contributions are thus not weighted equally.This explains why their Λ values deviate from that of n/n 0 .The order of magnitude of all 1/Λ values is a few days.1.
Although particular values mentioned here might depend on the resolution of meteorological data used, on the interpolations applied and on the numerical scheme of the trajectory integration, one feature is independent of these details, namely that the deviation between the appropriate BWC and FWC becomes significant if the integration time exceeds the Lyapunov time of advection (the reciprocal of the Lyapunov exponent of advection) characteristic to the region and to the atmospheric conditions.This is a general property of any kind of chaotic processes and is in harmony with the observation that the Lyapunov time is the characteristic time interval beyond which predictions become unreliable [17,22].

Global Results
The case study provides an overview of the possible behavior of the investigated statistical parameters in a specific case.In order to illustrate the general validity of the dependence of the parameters on ∆t b and to gain a global picture and a hint of the typical accuracy of backward simulations, several pollutant clouds are initiated at different geographical locations as described in Section 3. In the next figures, besides the global distribution of the statistical parameters, the statistics of the tropical region (represented here with ϕ 0 = 0 • , 20 • N/S) and that of the mid-and high latitudes (defined as ϕ 0 = 40 • N/S, 60 • N/S, 80 • N/S) are also presented separately.
Figure 7a shows that the ∆CM h of all three investigated regions exhibits a clear exponential increase over ∆t b .The positive difference between mean (denoted by diamonds) and median values (horizontal slashes) refers to extreme cases when backward simulations give much more inaccurate results than for most of the simulations with the particular ∆t b .2.
A rapid growth can also be seen in the ∆σ h (∆t b ) diagram.For at most three-day backward simulations, ∆σ h remains below 10-50 km, which implies generally a good agreement in the extension of the BWC and FWC considering that σ h,f ranges between 33 and 5480 km with an average of about 1100 km.As in Figure 5, the increase of ∆σ h for small ∆t b values can be approximated by an exponential function.Similar to the case study, the long-term behavior can be fitted by a large-exponent power-law function, but an exponential function also fits to the entire interval.Here, we only show the latter fit.
The mean values of n/n 0 in Figure 7c exhibit an exponential decay as in the case study.The decay of the n(∆t b )/n 0 function is the slowest for the tropical region.The slower decrease is in agreement with the fact that the usual wind field of the tropics differs from that of the mid-and high latitudes, where permanently generated and dissipated cyclones dominate.The weaker shearing in the wind field in the tropics and, thus, the weaker chaoticity of advection results in more particles being able return to the initial volume.For small ∆t b for the tropics, negative deviation appears between the means and medians (mean minus median); while for large ∆t b for the mid-and high latitudes, positive deviations are the rule.This implies that, although, in the tropics, most of the pollutant clouds are advected relatively accurately back to their initial position, occasionally, n/n 0 can be small.Analogously, for a long time backward simulation, the chance of a particle of a pollutant cloud initiated in the mid-latitudes to get back to its initial cuboid is usually very small (note the narrow extension of the blue boxes in Figure 7c for ∆t b ≥ 5 days); there exist, however, specific cases when n/n 0 can even reach 0.85 (five days) and 0.55 (nine days) (not shown in Figure 7c).These features and the large uncertainty in n/n 0 in the tropics might be explained by the fact that, on the one hand, pollutant clouds in the tropics can become mixed into the extratropics (see, e.g., Figure 1) and vice versa, and on the other hand, clouds dispersing even in their original region can be subject to more intense or weaker circulation than the typical one in that region.
Figure 8 gives an overview of the statistical parameters FMS 3D and PCC 3D .As in the case study, FMS 3D and PCC 3D seem to decrease approximately exponentially over ∆t b .The property of simulations with excessively bad coincidence of FWC and BWC for the tropics for small ∆t b (where the difference between the mean and median is positive) and ones with better agreement than the average for mid-and high latitudes for large ∆t b (where the mean minus median difference is negative) appear also in Figure 8.The overlap of the BWC and FWC reduces below 50% after six days of backtracking for the tropics and after three days for mid-and high latitudes.For large ∆t b , the FMS 3D , as well as the PCC 3D is small; however, there is a factor of about eight and six between the values of tropical and non-tropical values, respectively.Table 2 presents the Λ values of the fitted exponential functions to the mean data.It reveals that the time-dependence of the different statistical parameters can be described by similar Λ values.2.
Table 2.The exponent Λ (day −1 ) and the error of the exponential fits to the curves ∆CM h , ∆σ h , n/n 0 , FMS 2D , FMS 3D , PCC 2D and PCC 3D versus ∆t b for the globe (GL), for the tropical region (TR) and for mid-and high latitudes (M/H) in Figures 7 and 8. -

Outlook
Although, in large-scale dispersion simulations, the effect of small-scale 3D turbulent diffusion on particles (taking effect below the grid-size of the meteorological data) can generally be neglected compared to that of advection, in order to obtain insight into how turbulent diffusion alters our findings, we repeated the simulation of the case study with a typical constant horizontal diffusivity of K h = 50 m 2 •s −1 [31].Adding turbulent diffusion to the trajectory calculations can be considered as adding noise to the system.The effect of turbulent diffusion in the forward simulation over time t is enhanced further over an additional time interval of ∆t b in the backward simulation.It broadens the extension of the pollutant cloud in backward mode, too, as also found in other dispersion models (see, e.g., [31]).
Table 3.The exponent Λ (day −1 ) and the error of the exponential fits to the curves ∆CM h , ∆σ h , n/n 0 , FMS 2D , FMS 3D , PCC 2D and PCC 3D versus ∆t b of the case study described in Section 5, but with horizontal turbulent diffusivity K h = 50 m 2 •s −1 .The deviation ∆Λ of the Λ values of the case study with turbulent diffusivity from the case study without turbulent diffusivity is also shown.

-Λ ∆Λ
∆CM h 0.498 ± 0.026 0.017 ± 0.058 ∆σ h 0.481 ± 0.050 0.014 ± 0.100 n/n 0 0.379 ± 0.024 0.040 ± 0.056 FMS 2D 0.514 ± 0.037 0.048 ± 0.085 FMS 3D 0.538 ± 0.023 0.082 ± 0.063 PCC 2D 0.464 ± 0.023 0.017 ± 0.051 PCC 3D 0.447 ± 0.022 0.010 ± 0.049 Table 3 gives the exponents Λ of the curve fittings gained from the case study when adding turbulent diffusivity to the simulations.∆Λ indicates the difference from the Λ values of Table 1.As expected, e.g., for n/n 0 , a somewhat faster decrease shows up, and the n/n 0 values are found 0.02-0.1 smaller at each ∆t b than in the case without turbulent diffusion.Nevertheless, the change in Λ even for this quantity falls in the error bar of the fit.There are nearly no changes, neither in the data nor in Λ for ∆CM h , as turbulent diffusion has no direct impact on the center of mass.Only a slight increase is found for ∆σ h with about 30 km (∆t b = 4 days) to about 100 km (∆t b = 9 days) in the data.This is due to the fact that turbulent diffusion enhances the extension of both the forward and the backward cloud; however, in the latter case, diffusion acts on a time interval longer by ∆t b .Nevertheless, neither does this significantly alter the value of Λ.It confirms the assumption that K h = 50 m 2 •s −1 has negligible effect on large-scale dispersion.The mean values of PCC 2D/3D change about ±0.03, while for FMS 2D/3D , a definite decrease of 2%-10% appears for backward simulations of 1-5 days duration.
The results obtained in the presence of weak horizontal turbulent diffusion are practically the same as without any diffusion.This is consistent with the well-known property of low dimensional chaos according to which a weak noise washes out fine details below a characteristic length scale, but does not influence at all properties beyond that scale [16,40].
If the particles represent not passive tracers, rather aerosol particles with realistic radius r and density ρ p (e.g., r = 1 µm, ρ p = 2000 kg•m −3 ), Equation (1) contains an additional term describing the gravitational settling of the particles [25].In this case, particles can deposit on the ground in both the forward and the backward simulation, and during the backward one, particles that deposited in the forward run are treated as particles emitted from the surface in the appropriate time instant and location.This additional process shortens the time that the particles spend in the atmosphere before deposition (i.e., the time of the forward and backward tracking); nevertheless, our preliminary results show that the characteristics of the accuracy of backward tracking do not change.

Summary
In order to study the accuracy of backward trajectory simulations for locating the initial position of pollutant clouds precisely, we carried out, besides a case study, several simulations of pollutants distributed over the globe.Some characteristics of the backward and the forward pollutant clouds were determined at the same time instant.Forward pollutant clouds were considered as "real" clouds with which the results of the backward simulation are expected to coincide under full time-reversibility.
The results reveal that for ∆t b = 1-3 days of backward simulations, a reasonable agreement of the backward and forward clouds can be found, that is the backward tracking for such time periods performs reasonably well.This agreement can be seen in the figures of the simultaneous plots of the backward and forward clouds, and statistical parameters also support this observation.The difference between the center of mass/standard deviation of the two pollutant clouds grows exponentially with the time interval of the backward simulation at any initial location over the globe, which can be considered to be a consequence of the chaotic motion of the particles.The accuracy of backward simulations proves to decrease exponentially, too, by studying the proportion of the particles that are able to get back to their initial volume, the figure of merit in space or the Pearson's correlation coefficient between the backward and forward cloud.A faster decrease characterizes the mid-and high latitudes compared to the tropical region.The above-mentioned observations are also valid for simulations that take into account the effect of turbulent diffusion, as well.
Note that the time interval over which the accuracy of time-reversibility is found to be reasonable is on the order of magnitude of the predictability time of atmospheric dispersion (the reciprocal of the Lyapunov exponent), and furthermore, this is on the same order as the predictability time of the wind field [41].
We conclude that aside from rare occasions with calm weather conditions and weak wind shear (especially in the tropics), with the resolution of the meteorological data used here, backward tracking of pollutant clouds is not reliable beyond a few days.On time scales longer than this, simulations may yield an error of about several 100 km in the position of the center of mass, and the overlap and Pearson's correlation coefficient also fall below e −1 , usually considered as a limit of predictability.
Although results in the coincidence of the forward and backward clouds may be improved by applying more sophisticated numerical schemes and meteorological data with finer resolution at the expense of computational cost, it is not expected to basically alter the dynamics of the time-dependence of the quantities investigated, due to the unavoidable exponential degradation of the accuracy of dispersion simulations in our chaotic atmosphere.

Figure 1 .
Figure 1.Forward (black) and backward (colored) pollutant clouds emitted instantaneously and consisting of n 0 = 1.25 × 10 5 particles at (a) t = 0 and (b) t = 7 days.The forward pollutant cloud is initiated on 12 March 2011 at 00 UTC at 120 • E, 20 • N, 500 hPa in a 1 • × 1 • × 100 hPa cuboid beyond Indonesia.Colored dots denote the location of the particles whose backward simulation started at time t b = 2 days (blue), t b = 7 days (yellow) and t b = 9 days (red), respectively.In the inset of panel (a) instead of black dots the initial position of the particles is illustrated by a black rectangular frame.This simulation is an element of the global study in Section 6.

Figure 2 .
Figure 2.An outline of the time schedule of the forward and backward simulations.Forward simulation starts at Time 0 and terminates after nine days.Backward simulations start at time t b (where the value of t b can be 1, 2, . . ., 9 days) and terminate at Time 0. The comparison of the forward and backward cloud is carried out at time t (t = 0, 1, . . ., t b − 1 days) as a function of the time period ∆t b of the backward tracking.

Figure 3 .
Figure 3. Forward (black) and backward (colored) pollutant clouds emitted instantaneously and consisting of n 0 = 1.25 × 10 5 particles of the case study at (a) t = 0, (b) t = 2, (c) t = 5 and (d) t = 7 days.The forward pollutant cloud is initiated on 12 March 2011 at 00 UTC at 141 • E, 37.5 • N, 500 hPa in a 1 • × 1 • × 100 hPa cuboid beyond Fukushima.Dots colored according to the color bar on the right denote the location of the particles whose backward simulation started at time t b (day).In the inset of panel (a), instead of black dots, the initial position of the particles is illustrated by a black rectangular frame.For better visibility, only every fifth particle of the clouds is plotted.

Figure 4 .
Figure 4. (a) The horizontal standard deviation σ h,f of the forward cloud (FWC) in Figure 3 as a function of time t with an exponential fit over t = 0-2 days (dashed line) and a power-law fit over t = 2-9 days (dash-dotted line); (b) the horizontal standard deviation σ h,b of the backward clouds (BWCs) in Figure 3 as a function of the time interval of the backward simulation ∆t b = t b − t.Solid lines represent the time evolution of the BWCs initiated at t b (see color legend) backward in time; dashed lines connect the values corresponding to time t.

Figure 5 .
Figure 5. (a) The horizontal difference ∆CM h between the center of mass of the BWCs and the FWC depending on the time period ∆t b of the backward simulation; (b) the difference ∆σ h between the horizontal standard deviations of the BWCs and FWC; (c) the proportion n/n 0 of the particles returned back to the initial cuboid.Colored triangles indicate the starting time t b of the backward simulations (see color legend); black diamonds (in some cases in overlap with triangles) are the mean values belonging to a given ∆t b (denoted by "m").Dashed lines indicate exponential fits to the mean values; the exponents are given in Table 1.The dash-dotted line in panel (b) indicates a power-law fit to mean values over the interval ∆t b = 3-9 days with an exponent of 4.3.

Figure 6 .
Figure 6.(a) The figure of merit in space FMS 3D between the BWCs and FWC depending on the time period ∆t b of the backward simulation; (b) Pearson's correlation coefficient PCC 3D between the BWCs and FWC.Colored triangles indicate the starting time t b of the backward simulations (see the color legend); black diamonds (in some cases in overlap with triangles) are the mean values belonging to a given ∆t b (denoted by "m").Dashed lines indicate exponential fits to mean values; the exponents are given in Table1.

Figure 7 .
Figure 7. (a) The horizontal difference ∆CM h between the center of mass of the BWCs and the FWC depending on the time period ∆t b of the backward simulation; (b) the difference ∆σ h the horizontal standard deviations of the BWCs and FWC; (c) the proportion n/n 0 of the particles returned back to the initial cuboid.The results are shown for the tropical region (TR, red), for midand high latitudes (M/H, blue) and for the globe (GL, black).Mean values (diamonds), medians and lower and upper quartiles (boxes with horizontal lines) are indicated.Dashed lines indicate exponential fits to the mean values; the exponents are given in Table2.

Figure 8 .
Figure 8.(a) The figure of merit in space FMS 3D between the BWCs and FWC depending on the time period ∆t b of the backward simulation; (b) Pearson's correlation coefficient PCC 3D between the BWCs and FWC.Mean values (diamonds), medians and lower and upper quartiles (boxes with horizontal lines) are indicated for the tropical region (TR, red), for mid-and high latitudes (M/H, blue) and for the globe (GL, black).Dashed lines indicate exponential fits to the mean values; the exponents are given in Table2.