Escape rate : a Lagrangian measure of particle deposition from the atmosphere

Due to rising or descending air and due to gravity, aerosol particles carry out a complicated, chaotic motion and move downwards on average. We simulate the motion of aerosol particles with an atmospheric dispersion model called the Real Particle Lagrangian Trajectory (RePLaT) model, i.e., by solving Newton's equation and by taking into account the impacts of precipitation and turbulent diffusion where necessary, particularly in the planetary boundary layer. Particles reaching the surface are considered to have escaped from the atmosphere. The number of non-escaped particles decreases with time. The short-term and long-term decay are found to be exponential and are characterized by escape rates. The reciprocal values of the short-term and long-term escape rates provide estimates of the average residence time of typical particles, and of exceptional ones that become convected or remain in the free atmosphere for an extremely long time, respectively. The escape rates of particles of different sizes are determined and found to vary in a broad range. The increase is roughly exponential with the particle size. These investigations provide a Lagrangian foundation for the concept of deposition rates.


Introduction
There are several Eulerian and Lagrangian dispersion models that simulate and forecast the movement of air pollutants in the atmosphere by using meteorological data.Dry and wet deposition processes and their parameterization are some of the main issues in these models.First, we summarize these traditional approaches.
In Eulerian ones, based on the numerical solution of an advection-diffusion-sedimentation partial differential equation, the deposition processes are included generally in the form ∂c(r, t) ∂t = − (k d (r, t) + k w (r, t)) c(r, t), where c is the concentration of a pollutant at a given location and time, and k d and k w are the dry and wet deposition coefficients (the latter is also called the scavenging coefficient).
Lagrangian particle-tracking models can be divided into two classes: models in which an artificial mass is assigned to any particle and this mass depends on time (we refer to them in this paper as "ghost" or "computational" particles), and models that follow "real particles" with fixed, realistic size and density.
The latter ones (such as PUFF, Searcy et al., 1998 andVAFTAD, Heffter andStunder, 1993) are typically designed to predict the dispersion of volcanic ash as quickly as possible.Therefore, in these models some physical processes such as dry and wet deposition are neglected, as they are not very important higher in the atmosphere where the emitted ash mainly spreads.
In the case of ghost particle models (see, e.g., HYSPLIT, Draxler and Hess, 1998, 2004, FLEXPART, Stohl et al., 2005, NAME, Ryall and Maryon, 1998;Jones et al., 2007;Webster and Thomson, 2011, MLDP0, D'Amours and Malo, 2004, GEARN, Terada and Chino, 2008), the mass m carried by a ghost particle decreases due to dry and wet deposition if the particle travels through a region where these processes are present (e.g., close to the ground, and in clouds or in regions of precipitation).The impact of these processes can be described (in a similar spirit as in the Eulerian approach) along the path of a ghost particle as: Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
T. Haszpra and T. Tél: Escape rate in the atmosphere where C(r(t), t) can be a location and time-dependent dry or wet deposition coefficient.This equation, with a constant C, results in an exponential loss of the mass.In this approach a ghost particle is considered to be the center of mass of a great amount of adjacent pollutants.Since, however, the chaotic nature of the advection dynamics implies that an initially small, compact ball of particles becomes rapidly deformed into a complicated, filamentary shape of large extent in the atmosphere, the physical reality of ghost particle models remains questionable.
In this paper we, therefore, carry out simulations with real, spherical particles, but extend the above-mentioned real particle models to include boundary layer processes: to reckon with wet deposition as a stochastic process of individual particles and to take into account the effect of turbulence as a random walk.This extended approach is called the Real Particle Lagrangian Trajectory (RePLaT) model.Since the dynamics of aerosol particles is then a kind of dynamical system, concepts from chaos theory (like, e.g., topological entropy, Haszpra and Tél, 2013) can be applied.Here we propose the use of the escape rate as a measure of the speed of the deposition process from the atmosphere.
Our results reveal a considerable variance in the escape rate for particles with different radii, i.e., they unfold the importance of different gravitational settling velocities.Moreover, they make clear that the winds are essential, and the deposition dynamics is much more complicated than a simple settling that would only occur in motionless air.The impact of turbulent diffusion and wet deposition on the escape dynamics is also investigated.
The paper is organized as follows.In Sect. 2 we provide a brief definition of escape rates.Section 3 presents the equations of motion for aerosol particles advected in a given wind field.Turbulent diffusion is parametrized as a random walk process.The usual parameterization of wet deposition is extended here to real particles by simply allowing an aerosol particle to be converted in a much larger raindrop with certain probability depending on the local rain intensity.The data and numerical methods used in the simulations are given in Sect. 4. The results obtained for the deposition dynamics and for the dependence of the escape rate on particle radius, both in the free atmosphere and in the boundary layer, are given in Sect. 5. Section 6 presents a case study of a hypothetical eruption of the Merapi volcano.Discussion and outlook are given in Sect.7.

Escape rates
The state of the atmosphere is complicated both in space and time, and can therefore be called turbulent in Eulerian terms (Holton, 1992).The motion of an individual particle is described by an ordinary (or weakly stochastic) differential equation in which the wind velocity enters as a known input function.This is an example of a dynamical system, the solution of which is typically chaotic.The particle trajectory is a Lagrangian quantity.The advection of particles can be considered as a paradigmatic aspect of the chaoticity of the atmosphere.Concepts coming from chaos theory can thus be taken over with the hope of their successful applicability to the problem of particle dispersion.
Under certain circumstances chaotic behavior is of finite duration, i.e., the complexity and unpredictability of the motion can be observed over a finite time interval only.Nevertheless, there also exists in such cases a set in phase space responsible for chaos, which is, however, non-attracting.This type of chaos is called transient chaos and the non-attracting set is a chaotic saddle (for an introductory text see Tél and Gruiz, 2006).Since there are typically significant differences in the individual lifetimes, an average lifetime can be defined.To this end, it is worth following several motions instead of a single one: the study of particle ensembles is essential.To characterize the dynamics, one takes a preselected region, and starts n 0 1 trajectories in it.They escape the preselected region sooner or later, and the motion before escape appears to be chaotic.The number n(t) of trajectories that never left the preselected region up to time t is thus a monotonically decreasing function of t.After a sufficiently long time (for t larger than some t 0 ), the decay in the number n(t) of survivors is generally exponential (similar to the law of radioactive decay): n(t) ∼ exp(−κt), for t > t 0 . (3) Coefficient κ is called the escape rate (Ott, 1993;Tél and Gruiz, 2006;Lai and Tél, 2011).Its reciprocal value can be considered as an estimate of the average lifetime of chaos.A nonzero escape rate is thus a new, important chaos characteristic: the larger the value of κ, the faster the escape/sedimentation process.
There might be situations in which two chaotic saddles coexist.This leads to the appearance of escape rates κ 1 and κ 2 .In such cases, the number of survivors n(t) in the preselected region is the sum of two exponentials for t > t 0 (Lai and Tél, 2011): (4) For escape rates of non-infinitesimally differing values, this implies that the decay goes for large times with the smaller of the two escape rates and that the majority of particles escapes according to the larger one.
In the atmospheric context, the preselected region might be the entire atmosphere.The condition of escape is then the first arrival at the surface.We shall see that a separation of time scales is typical, and two considerably different escape rates characterize each deposition process.Our basic interest will be how the escape rates depend on the size of aerosol particles.Furthermore, we shall see that turbulence and wet deposition influence the escape dynamics.We claim that the escape rates provide a kind of Lagrangian characterization of the entire deposition process.1 3 Equations of motion

Free atmosphere
The motion of small, heavy spherical particles of radius r is determined by the sum of the gravity and the Stokes drag.Buoyancy is negligible since it is proportional to the ratio of the density ρ of air and ρ p of the particle (which is less than or equal to 1/1000 for typical atmospheric aerosol particles).The dimensionless form of Newton's equation (Maxey and Riley, 1983) is then where ṙp is the velocity of the particle and v(r, t) is the velocity of the ambient air at the location r of the particle at time t, while w term is the terminal velocity in motionless air, and n is a unit vector pointing upwards.
where L and U represent a characteristic distance and velocity, respectively, and ν is the kinematic viscosity of air.
Since we are interested in phenomena on length scales L ∼ 10-1000 km, and with wind speeds U ∼ 1-50 m s −1 , the Stokes number for ν ≈ 10 −5 m 2 s −1 , ρ p = 2000 kg m −3 , and r = 12 µm (the largest size we investigate) is St ≤ 3 × 10 −5 .The left-hand side of Eq. ( 5) can thus be neglected indicating that the motion takes place practically under the balance of the Stokes drag and gravity.In other words, the particle velocity becomes immediately equal to the terminal velocity superimposed on the wind velocity.
The dimensional form of the equation of motion can then be written in the form: The Stokesian terminal velocity for heavy particles of radius r (certainly valid for r ≤ 12 µm) is where g denotes the gravitational acceleration.
Since the meteorological fields used for the simulations (see Sect. 4) are given on pressure levels, we determine trajectories in pressure coordinates.The vertical component of Eq. ( 7) is żp (t) ≡ w p (t) = w(r p (t), t) + w term .
(9) Similar to Eq. ( 9), the vertical motion of a particle in pressure coordinates can be written as dp(r p (t), t) Using hydrostatic approximation, the terminal velocity ω term in pressure coordinates is found to be ω term = −ρgw term , and substituting w term from Eq. ( 8) we obtain as in Haszpra and Tél (2011) that Note that deviations from the hydrostatic approximation are known to be negligible on length scales larger than about 30 km (see, e.g., Kalnay, 2003), and since the typical spatial resolution of our database is 100 km (see Sect. 4), the use of hydrostatic relations is assured in our case.
The location dependence of the kinematic viscosity ν of air is due to temperature T and pressure p, and its simplest form can be represented by Sutherland's law (Sutherland, 1893) Here β 0 = 1.458 × 10 −6 kg m −1 s −1 K −1/2 is Sutherland's constant, T S = 110.4K is a reference temperature, and R d = 287 J kg −1 K −1 is the specific gas constant for dry air.
Since in the horizontal directions the use of spherical coordinates is appropriate, we solve Eq. ( 7) in the form: dp p dt = ω(λ p (t), ϕ p (t), p p (t), t) where λ p and ϕ p are the longitude and latitude coordinates, p p (t) ≡ p(r p (t), t) is the pressure coordinate of a particle along its path, and R E is the radius of the Earth.The limit of r = 0 can be considered as the passive advection dynamics for air parcels, since ω term = w term = 0, and hence from Eq. ( 7) the equation of motion of passive tracers, ṙp (t) = v(r p (t), t), follows.
T. Haszpra and T. Tél: Escape rate in the atmosphere

Planetary boundary layer
While far from the surface, on large scales, the effect of turbulent diffusion is typically negligible (for an estimate see Haszpra and Tél, 2013), this process plays an important role in the dispersion of particles in the planetary boundary layer (PBL).Similarly, precipitation in this region is also relevant.

Turbulent diffusion
In atmospheric dispersion simulations, constant horizontal diffusivity K h is often applied, while the vertical diffusivity K z has altitude dependence (Holtstag and Boville, 1993).We also apply this approximation, and based on Visser (1997) and Terada and Chino (2008), we take turbulent diffusion into account as a random walk process.The equations of motion are integrated by Euler's method and can be written as: We take R as a random process uniformly distributed between −0.5 and 0.5 (the standard deviation is σ R = 1/12), like in GEARN (Terada and Chino, 2008).Since K z depends on the altitude, Eq. (14c) compared to Eq. ( 14a) and (14b) has an additional drift term (the third term) advecting particles from regions with low diffusivity to regions with high diffusivity.According to the Monin-Obukhov similarity theory (Troen and Mahrt, 1986;Holtstag et al., 1990;Holtstag and Boville, 1993), the vertical diffusivity increases close to the surface and decreases when approaching the top of the boundary layer from below.Therefore particles are advected on average upwards from the bottom (where dK z /dz > 0) and downwards from the top (where dK z /dz < 0).
The equations of motion (13a)-(13c) in spherical and pressure coordinates solved by Euler's method in the presence of turbulent diffusion can be written as: where K λ and K ϕ are and K p is calculated from K z obtained from the Monin-Obukhov similarity theory.The altitude dependence of K z is calculated from the following form (Troen and Mahrt, 1986;Holtstag et al., 1990;Holtstag and Boville, 1993): K denotes the von Kármán constant, u * is the frictional velocity, and z PBL is the height of the planetary boundary layer.φ z L represents the similarity function that depends on the Monin-Obukhov length L: characterizing the stability of stratification.The dynamic temperature T * and frictional velocity u * are obtained from the following relationships (Högström, 1988;Wotawa et al., 1996): Here H represents the sensible heat flux, c p is the specific heat of air at constant pressure, τ EW and τ NS are the East-West and North-South surface stresses, respectively.For simplicity, we apply the similarity functions of Dyer (1974) 4 for unstable conditions (L < 0), and φ z L = 1 + 5 z L for stable conditions (L > 0) for any z ≤ z PBL .Because the vertical coordinate of the particles is the pressure level p p , not their altitude z p , the pressure level p PBL corresponding to z PBL has to be calculated.Assuming that ρ = const in the boundary layer and applying the hydrostatic approach again, where p s denotes the surface pressure, which is the pressure of the lowest level (1000 hPa) in the database used.(Taking the density changes also into account would lead to a 0.1 % difference in the value of p PBL only.)Since z = (p s −p)/(ρg) and z PBL = (p s −p PBL )/(ρg) are good approximations in the boundary layer, the vertical diffusivity in (17) can be rewritten as a function of p as: This is transformed into the vertical diffusivity in pressure coordinates by multiplying Eq. ( 21) by (ρg) 2 : If a particle reaches the lowest level (p s = 1000 hPa), it is considered to have escaped, to be deposited on the surface.There is thus no need to develop any dry deposition parameterization, as a natural consequence of our approach.

Wet deposition
In Eulerian models the impact of wet deposition is taken into account by means of Eq. ( 1) with k d = 0.This implies that after a short time t, a fraction 1 − exp(−k w t) ≈ k w t of the particles remains in the cell, where k w is the wet deposition coefficient (scavenging coefficient).
We use this relationship to incorporate wet deposition into our Lagrangian real particle model.We consider wet deposition as a random process that results in a particle being captured by a raindrop with probability p = 1 − exp(−k w t).Thereby the radius of the particle suddenly increases to the mean radius r rain of raindrops.Then in the equations of motion a new terminal velocity is calculated with r = r rain and ρ p = ρ rain = 1000 kg m −3 (the weighted mean of density and radius of an aerosol particle and a raindrop for the "new" particle is found to differ by less than 0.5 %).Raindrops containing aerosol particles can also collide, but such second collisions are unlikely events owing to the short time the raindrops spend in the atmosphere.Since for typical raindrops the Reynolds number evaluated with the terminal velocity is much larger than unity, w term and ω term are to be calculated from the quadratic drag force: ) where C D is the drag coefficient (≈ 0.4 for a sphere).
The scavenging coefficient is proportional to the number of aerosol particles collected by raindrops per unit time.Since r rain r and therefore |w term (r rain )| |w term (r)|, this coefficient is (Sportisse, 2007): where r 2 rain π is the collision area, E(r rain , r) is the collision efficiency and nrain is the number of raindrops per volume.The collision efficiency is defined as the fraction of the particles of radius r scavenged by raindrops of radius r rain in a volume.
The rain intensity can be expressed as (Sportisse, 2007) The scavenging coefficient k w thus appears in terms of P as: There are several parameterizations for the typical raindrop radius as a function of rain intensity P , generally in the form r rain = a P b (Sportisse, 2007).We use the Pruppacher-Klett parameterization (Pruppacher and Klett, 1998): where the unit of r rain is mm and the unit of P is mm h −1 .Assuming a typical constant collision efficiency E(r rain , r) = 0.1, and estimating the typical radius of the raindrops by the Pruppacher-Klett formula, the scavenging coefficient is where the unit of P is still mm h −1 .
To sum up, if there is precipitation of intensity P at the location of an aerosol particle, its radius suddenly becomes in the simulation r = r rain = 0.488 P 0.21 with the probability of p = 1−exp(−k w t), where k w is given by Eq. ( 28).With P = 1 mm h −1 , k w = 0.154 h −1 = 3.7 day −1 and with a time step of about 5 min, p = 0.0128.For simplicity, the effect of wet deposition is taken into account only below the 850 hPa level.
We have thus incorporated wet deposition in a simple manner in our RePLaT model, motivated by the Eulerian approach.Note, however, that the latter is independent of the winds, therefore the escape rate obtained in our Lagrangian picture might be different from k w .Moreover, the fact that only a small portion of aerosol particles becomes converted into raindrops makes the difference between k w and the escape rate even more pronounced.

Data and methods
Wind components u, v and ω, boundary layer height z PBL , sensible heat flux H , precipitation P , temperature T and East-West and North-South surface stresses τ EW and τ NS are taken from reanalysis fields of the ERA-Interim database of the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011).The meteorological variables are available at 22 pressure levels between 1000 and 100 hPa on a 1.5 • × 1.5 • horizontal grid with a 6 h time resolution.The time period 1 January to 31 December 2010 is considered.
In order to compute trajectories, the data on the regular grid are interpolated to the location of the particles, firstly using linear interpolation in time, then linear interpolation in vertical, and bicubic spline interpolation in horizontal.The only exception is the height of the boundary layer for which in the full afternoon period 12:00-18:00 UTC the value of 12:00 UTC is used, since the collapse of PBL is rather fast and the use of linear interpolation would lead to a strong underestimation of z PBL in this period (Stohl et al., 2005).Particle trajectories are determined from Eq. ( 15a)-(15c) with a time step of t = 5.625 min.In each time step, if at the location of a particle with p p > 850 hPa precipitation is P > 0, the method described in Sect.3.2.2 is applied.

Deposition dynamics, separation of time scales
In order to determine global escape rates, we distribute n 0 = 2.5×10 5 particles uniformly over the globe on different pressure levels on 1 January 2010.They are tracked up to their escape, but at longest for 1 yr.To study the dependence of the escape rate on the particle size and on the initial altitude, simulations are run with radii of r = 0, 1, 2, . . ., 12 µm and initial altitudes of p 0 = 500, 700, 850 and 900 hPa.Note that the radius of a particle can suddenly change according to Sect.3.2.2 if the particle is captured by a raindrop.The limiting case of a "particle" with r = 0 µm may be considered as a gaseous contaminant in the atmosphere that can be dissolved in raindrops with some probability, which, in a first approximation, does not depend on the properties of the "particle".
To compare different effects, simulations are carried out in three setups that take into account: 1. advection, turbulent diffusion and precipitation, 2. advection and turbulent diffusion, and 3. only advection.
As a first example, Fig. 1 exhibits the number of survivors vs. time for different initial altitudes p 0 and for various particle radii r in the three setups.Panels a, b and c correspond to free atmospheric initial conditions above the 850 hPa level.As the aerosol particles are initially far from the surface, the curves start with a plateau: no outfall from the atmosphere takes place within the first few days.This phenomenon is present also for particles initiated lower in the atmosphere if turbulent diffusion and precipitation are switched off, since these effects cannot influence deposition (Fig. 1d, green curve).After a short transition following the plateau (for t > t 0 ≈ 1-15 days), an approximately exponential decay can be seen in all of the three setups for a few days (for an example see the dashed line belonging to days 2-5 in Fig. 1b).After some time, however, a crossover takes place and a slower exponential decay sets in (see, e.g., Fig. 1b for t > 10 days).Thus, we can speak of a short-term and a long-term exponential decay taking place with different exponents.The corresponding escape rates will be denoted by κ s and κ , respectively.The κ s and κ values extracted from the data of Fig. 1 are summarized in Table 1 in the unit of day −1 .Escape rates can be used as measures of the deposition process.It is indeed striking to see that any escape rate is at least 10 times larger for large aerosol particles (9 or 10 µm) than for small ones.The deposition process is thus very fast for large sizes.At any given size, the long-term escape rate is at least half or smaller than the short-term one.Since this difference appears in the exponent, we can safely speak about a separation of time scales in the deposition process.
For particles initiated below 850 hPa, precipitation and/or turbulent diffusion influences their motion from the very beginning.Thus, no plateau can be found (Fig. 1d, blue and red lines) as there are always regions on the globe where precipitation takes place and/or the particles are in the boundary layer where they are also subject to turbulent diffusion.
Figure 1a, c and Table 1 demonstrate that for small particles initiated at any height in the atmosphere, the effect of rain and turbulent diffusion plays an important role in the deposition process and intensifies the outfall.The effect is also present for small particles in the boundary layer (not illustrated here).For r 5 µm emitted above 850 hPa, the differences between the curves disappear (Fig. 1b).For particles with r 10 µm initially below 850 hPa, the quickest depletion is found in setup 3, without turbulent diffusion (Fig. 1d).These phenomena will be explained in the following paragraphs.
All our findings illustrate that the naive expectation coming from dynamical systems theory according to which the global emptying is a random process described by a single exponential decay does not hold.In the atmosphere, instead, a short-term and long-term dynamics can be identified, characterized by two different approximately exponential decays.

Dependence of the escape rates on the radius and initial altitude
Figure 2 shows the long-term escape rate κ (obtained as linear fittings to the curves ln(n/n 0 ) vs. t in the asymptotic linear regime) as a function of particle radii r for the four initial pressure levels in the three setups investigated.Although the atmospheric decay dynamics obviously vary with the initial altitude, the slopes of the ln(n/n 0 ) curves, i.e., the escape rates, do not seem to be dependent on the initial level in either setup.The reason for this phenomenon can be the fact that particles surviving a long time in the atmosphere become well mixed.The independence of p 0 indicates that there exists a global atmospheric chaotic saddle, and the long-lived particles reflect properties of this set underlying the deposition dynamics.κ (r) is thus a global atmospheric characteristic of particles of size r.The atmospheric saddle is likely to be time-dependent, and the κ (r) values are characteristic of the time period investigated.
It is remarkable that κ ranges over about two orders of magnitude although the radii vary over one decade only.The dependence is thus strongly nonlinear.The best approximate fit appears to be exponential Exponent k is found to be k ≈ 0.33-0.38µm −1 for setups including rain and/or turbulent diffusion (setups 1, 2), and k ≈ 0.44-0.47µm −1 otherwise.A comparison of Fig. 2c with Fig. 2a and b reveals that precipitation and turbulent diffusion strongly enhance the proportion of the outfalling particles with small radii; the escape rate grows by a factor of 2-3.It is worth comparing the scaling of Eq. ( 29) with a naive estimate.The time needed to pass a fixed vertical distance Z with the terminal velocity Eq. ( 8) is Z/ | w term |.Since the terminal velocity is proportional to r 2 and the reciprocal of this time corresponds to the escape rate, this estimate results in a scaling proportional to r 2 .The fit of this functional form to the data is much less satisfactory than that provided by Eq. ( 29).The difference between the power law behavior and the observed exponential one can only be interpreted by realizing that atmospheric winds play an essential role in the deposition process.
The short-term escape rates κ s are also determined for different particle sizes (Fig. 3).The dependence on the radius seems to remain basically exponential, but contrary to the long-term escape rate, κ s also depends on the initial pressure level.The particles responsible for the short-term behavior fall out rapidly; they have no time to visit the global chaotic saddle of the atmosphere, and experience a chaotic saddle characteristic of the individual initial altitude.
It is not surprising that, in setup 1 including rain, κ s is larger for lower initial levels (900 hPa and 850 hPa) than that for free atmospheric initial levels (700 hPa and 500 hPa), since close to the surface, precipitation enhances the deposition from the very beginning.In this setup (Fig. 3a) exponent k is much smaller (k ≈ 0.196 and 0.262 µm −1 ) for the two lower levels than that for the higher levels (k ≈ 0.336 and 0.322 µm −1 ).Precipitation accounts for this feature because it has the same effect for particles of any size in our approach, therefore these escape rates are significantly influenced by the escape rate of raindrops 2 , hence no strong 2 We note that the naive estimate for the time of raindrop outfall from the atmosphere is Z/ | w term | with Eq. ( 23) as the terminal velocity.It yields 1-10 min for Z = 1.5 km for r rain = 0.1-5 mm raindrops with | w term |∼ 1-10 m s −1 .Our simulations with a probability 1 (p = 1) conversion of particles into raindrops at locations with precipitation fully support this estimate and the corresponding ∼ r rain 1/2 scaling of the escape rate.Moreover, for smaller raindrop radii (less than 0.1 mm) where Eq. ( 8) provides the terminal  dependence on the particle radius can show up.Free atmospheric initial conditions in setup 1 are necessarily also affected by the scavenging process, but the overall influence is weaker since not all particles "feel" the effect, only those that reach the 850 hPa level, and experience precipitation at their location.In agreement with this, for setups 2, 3 without precipitation exponent k is not smaller or not much smaller for the lower two levels than for the two free atmospheric levels (Fig. 3b and c).
The influence of turbulent diffusion is similar, apart from the fact that turbulent diffusion does not necessarily intensify deposition.This is a consequence of the dK p /dp term in Eq. (15c).As mentioned below Eq. (14c), particles are subvelocity, the ∼ r rain 2 scaling is also recovered in the numerics.All this indicates that the naive estimate would also be correct for aerosol particle sizes of 0.1 mm or larger.The mentioned strong effect of winds is thus only present for aerosol particles of radii less than a few times 10 µm.
ject to an uplift from the bottom of the boundary layer (where dK z /dz > 0, dK p /dp < 0), and particles are advected downward from the top of the boundary layer (where dK z /dz < 0, dK p /dp > 0).This combined effect can be responsible for the fact that without rain, but with turbulent diffusion, k has some variability between the lower and the upper levels (Fig. 3b), but this difference is less considerable than in Fig. 3a.The reason for the smaller deviation can also be due to the fact that the vertical extension of the boundary layer has a diurnal cycle, so PBL is very shallow over a large area of the globe, and therefore few particles can get into it.Figure 3c shows that if neither precipitation nor turbulent diffusion is taken into account, k has no significant height dependence: the ratio of the ks at different initial levels varies between 0.935 and 1.145.
Precipitation and turbulent diffusion enhance again the proportion of the outfalling particles with small radii: the escape rate grows by a factor of 7-13 in setups 1 and 2 relative to setup 3 (see Fig. 3c   effect is thus more significant than for the long-term escape process.It is, however, interesting that for large particles, the escape rate ratio mentioned above becomes 0.1-0.25 for the lower levels, implying that for r 5 µm, turbulent diffusion reduces the number of outfalling particles.The reduction can again be due to the term dK p /dp.Since larger particles have larger terminal velocities (ω term ≈ 7 × 10 −2 -1.5×10 −1 Pa s −1 for r = 5-10 µm), these particles have more chance to approach the ground, and because of the sign and magnitude of dK p /dp (which is ∼ (−0.5)-(−1)Pa s −1 ), close to the surface they are advected upward again.Hence their deposition process slows down compared to the case when they are only subject to advection, and can leave the atmosphere quickly owing to their relatively large terminal velocity.Since smaller particles have smaller terminal velocities, which are less than or of the order of the vertical velocity component of air; they approach the ground slower than larger particles in the lack of turbulent diffusion.For small particles, the term dK p /dp plays an important role in the upper part of the boundary layer, and advects them much faster towards the surface.The random term in Eq. ( 15c) also gives small particles a chance to fall out, so their deposition process is accelerated compared to the case without turbulent diffusion.
It is an interesting question, how particle ensembles with an initially vertical distribution deposit.Our data also provide insight into this aspect, since the superposition of the 4 ensembles (with initial distributions on approximately spherical surfaces, on levels of fixed pressures p 0 ) models an initial ensemble also distributed vertically in the atmosphere.Adding up the four functions for the number of survivors in the most complete setup, 1, a new global ensemble is obtained in which a new feature shows up: the plateaus disappear (not shown).The reason is the lowest level where escape Table 2. Short-term escape rate κs for a particle ensemble with vertical extension (described in the te the average κs and standard deviation σ of the short-term escape rates of different initial pressure leve has been found to be present from the very beginning.In this unified ensemble, we again find a separation of time scales and a short-term and a long-term escape rate, κs and κ , respectively.The latter one coincides with κ investigated so far, since the escape process from the well-mixed state of the atmosphere does not depend on the initial level.The shortterm unified escape rate, κs , differs, however, strongly from any of the κ s determined earlier.The results for different particle sizes are summarized in Table 2.The fact that the global results are consistently much smaller than the individual ones can be explained by the disappearance of the plateau.The decay process starts earlier and we observe it to be immediately exponential (t 0 is practically zero) with a necessarily smaller slope.The table also indicates that the global short-term escape rate κs is therefore by far not the average κs of the 4 individual ones.It is often a full standard deviation away from the average.We can thus say that the global picture does not contain important details visible in ensembles initiated on individual pressure levels.When, after some t 0 , the deposition process starts in simulations with these ensembles, the decay is much stronger than in the global ensemble.The latter one, however, properly follows the overall emptying process of the atmosphere, and might be relevant for estimating the average lifetime of chaos.The size dependence of κs (r) is found again to be approximately exponential: with k ≈ 0.297 µm −1 .Besides the escape rates, we also determine the average residence time for different particles in the atmosphere.In  (a, b, c, d, e) illustrate the dispersion 1, 3, 5, 7 and 9 days after the eruption.Colorbar indicates the pressure level of the particles in hPa.(f) shows the total precipitation of the last 12 h at 00:00 UTC on 8 November 2010.Colorbar indicates here the total precipitation in mm/12 h. the language of dynamical systems theory, this is the average lifetime of chaos for particles of radius r. Figure 4 shows the statistics for the four initial pressure levels.The global values τ (r) (not shown) follow from the averages of the individual ones.We find that the average residence times τ (r) vary between 16 days (for r = 1 µm) and 22 h (for r = 12 µm).Moreover, they are given by the reciprocal of the global short-term escape rate, with good accuracy: with k as given above.The dependence of the residence time on the initial altitude (p 0 ) is more complex, as shown in Fig. 4. For different initial levels, the average residence time for small particles varies between a few days and about 40 days, and for large particles between 0.2 and 2 days.The difference between the maximum and minimum decreases with the radius and with the initial height.The reason for the latter is the fact that some fraction of the particles initialized low in the atmosphere can fall out quickly due to either turbulent diffusion or precipitation.For guiding the eyes, dotted and dashed lines www.nonlin-processes-geophys.net/20/867/2013/ Nonlin.Processes Geophys., 20, 867-881, 2013 represent exponential fits to ln τ ave and ln τ med vs. r, and scaling ∼ exp (−kr) is found for the different initial pressure levels with k ≈ 0.207-0.283µm −1 and 0.116-0.228µm −1 , respectively.We find that in cases without any plateau, 1/κ s provides a good estimate of the residence time.For cases when a relatively large portion of particles survives the shortterm deposition (i.e., for small particles initiated in the free atmosphere), the residence time is approximately 1/κ .

A case study
Mount Merapi in Indonesia had long-lasting eruption series in 2010 from late October to November, which caused disruption to air traffic in the surroundings, besides which several people had to be evacuated (see, e.g., Surono et al., 2012).To study the outfall dynamics of aerosol particles, instead of the continuous eruptions, we simulate only a single volcanic ash puff of column shape of size 1 • × 1 • × 400 hPa, centered at λ 0 = 110.44• E, ϕ 0 = 7.54 • S and p 0 = 500 hPa.
Figure 5 demonstrates the horizontal dispersion of the ash cloud including n 0 = 2.16×10 5 particles of r = 10 µm emitted at 00:00 UTC on 1 November and tracked over 10 days.In the first few days, the particles spread northward, then with an anticyclonic flow they start to move in a westerly and southerly direction.A small fraction of the particles leave the atmosphere close to the source in this period after the hypothetical "eruption" when they happen to reach the 850 hPa level (below which precipitation is taken into account) due to the frequent rainfall events above Indonesia and the surroundings.In the period of days 6-8 (6-8 November), particles reach a region of a cyclone with strong precipitation, therefore a large amount of particles are scavenged out by rain in this period.A comparison of Fig. 5d and f shows that the particle distribution on the surface (marked by brown) is strongly correlated with rain intensity.Within 8 days, the majority of the particles falls out from the atmosphere.Indeed, only small changes can be seen in the last two days in the deposition pattern, as Fig. 5e indicates.
Figure 6 shows the proportion n/n 0 of the survivors over a somewhat longer time interval than in Fig. 5. Two setups are studied: setup 1 with the effect of turbulence and rain (blue) and setup 2 without the effect of precipitation (red).As expected, precipitation has an important role in the deposition dynamics: the simulation with rain results in much quicker depletion of the particles than the one without rain.
In the period 6-8.25 days the relation between n/n 0 and t seems to be exponential, and we find κ s = 2.295 day −1 with and κ s = 1.393 day −1 without rain.These differ from the global values 1.026-1.891day −1 for 500-700 hPa.After day 11 another exponential decrease takes place with κ = 2.066 day −1 and κ = 1.739 day −1 with and without rain, respectively.These are at most formal analogs of the longterm escape rates discussed in the previous section, since there is no time to reach a globally well-mixed state before 25 Fig. 6.Proportion n/n 0 of the number of survivors of the eruption in Fig. 5 vs. time in simulations, including the effect of rain (setup 1) (blue) and without rain (setup 2) (red), with exponential fittings to the data (black dashed lines).The slope of the fitted curves are κ s = 2.295 day −1 and κ = 2.066 day −1 for the blue curve, and κ s = 1.393 day −1 and κ = 1.739 day −1 for the red curve.the outfall in this example.Their existence indicates, however, a clear time scale separation again.The considerable deviations from the global behavior studied in the previous section are due to the fact that the decrease of n/n 0 is strongly affected by the local events, since particles are initiated in a relatively small volume corresponding to a volcanic eruption, and expand to an area of only approximately 4000 km × 100 km after 10 days.
It is insightful to look at the vertical distribution of the particles over the time span followed.This can be seen in the form of a histogram for setup 2 in Fig. 7.The initially columnar shape is deformed into a Gaussian one that spreads as its center moves downwards.This behavior was also observed in a simple cloud model with aerosol particles (Drótos and Tél, 2011).It is remarkable, however, that after the center of the Gaussian distribution reaches the surface, and the majority of the particles is deposited, the small fraction of particles remaining aloft is distributed widely in the different layers.It is the fraction of these extreme survivors that is responsible for the second, long-term exponential decay observed.We believe that this wide altitudinal distribution of the extreme survivors is also the physical background of the time-scale separation described in the previous section (although the average deposition process is much slower there than in this particular case).
Figure 8 illustrates the dispersion from the volcanic eruption with the same initial condition, but with smaller particles (r = 5 µm).As expected, such particles spread and reach very different regions in the atmosphere.Entering at different vertical levels, they become subjected to different horizontal winds.The strongly localized ash cloud on the 3rd day (Fig. 8a   (Fig. 8b).It is worth mentioning that despite the simplifying one-puff assumption, this figure shows good agreement with the satellite image of sulfur dioxide tracers in the period 4-8 November (http://earthobservatory.nasa.gov/NaturalHazards/view.php?id=46881).20 days after the hy-pothetical emission, the particles initialized in a small volume cover a huge area and are well mixed in the midlatitudes of the Southern Hemisphere.Therefore the long-term escape rate for this case (κ = 0.103 µm −1 ) is almost the same as the global escape rate κ for r = 5 µm particles.A remarkable T. Haszpra and T. Tél: Escape rate in the atmosphere feature of Fig. 8c is that the distribution of the deposited (brown) particles is fractal-like.There are large regions without any outfall, and the overall pattern is filamentary.The set of particles on the surface appears to trace out the intersection of the unstable manifold of the atmospheric chaotic saddle with the surface.This saddle might in principle be time-dependent, and what we see here is the set of these intersections in the period 7 to 20 days.

Discussion
We have illustrated that escape rate, a concept well known from the theory of transient chaos, can be usefully applied to the understanding of the deposition process of atmospheric aerosol particles.We have found a time-scale separation with two different escape rates.The short-term value characterizes the majority of particles falling out in a bulk, while the longterm escape rate characterizes the subensemble of extreme survivors.Escape rates can be used to estimate average residence times in the atmosphere.
Let us finally turn to the relation between escape rates and deposition rates.In our RePLaT model we can directly measure the number of deposited particles per unit time.This is much more natural than associating a mass-loss rate with a volume of air or with a ghost particle.In view of the long-term decay n ∼ exp(−κ t) found in our simulations, dn/dt ∼ −κ n follows.One might naively think that this is exactly Eq. ( 1) with k d + k w = κ .This is, however, not the case, since Eq.(1) applies not to the number of all particles in the atmosphere, but rather to those below a certain level, say 900 hPa only.To check the consistency of Eq. (1) with our approach, we evaluate n/(n 900 t) as a function of time for different particle sizes.Here n is the number of deposited particles over time step t ( n/ t ≈ −dn/dt) and n 900 is the number of particles below the 900 hPa level, the level where particles are initiated with uniform distribution over the globe.From Eq. (1) we expect that this quantity is k d + k w , and this deposition rate is at most weakly dependent on time as a possible effect of the temporal change of precipitation over the globe.The results are shown in Fig. 9 for radii r = 1 and 12 µm.
One can see that n/(n 900 t) is not constant and, even after smoothing, there is a pronounced time dependence.This is, however, of a different character for the two types of particles.This indicates that the reason cannot be the overall precipitation alone.The average of n/(n 900 t) (with t = 338 s) in the time period investigated is about 3 × 10 −5 s −1 = 2.59 day −1 for the large particles.This provides an analog of the deposition coefficient of the standard parameterization.The corresponding value for the small particles is about one third of this.These values happen to be of the order of the used deposition coefficients (Sportisse, 2007), but as we see, there are strong deviations from the average.parameters are the same as in Fig. 5).Panel a, b, c illustrate the dispersion 3, 7 and 20 days afte Colorbar indicates the pressure level of the particles in hPa.All in all, the parameterization of even wet deposition in the form of Eq. ( 1) is not without problems.With the increase of computer power we think that our RePLaT approach can overcome the ghost particle approach with the strange feature of an artificially defined mass.Some aspects can be improved, of course, thus, e.g., a more detailed in-cloud and below-cloud parameterization of the wet deposition process remains a subject for future work within the RePLaT model.
In summary, we have found that the emptying process of aerosol particles cannot be characterized by a single exponential decay.The global emptying process, from any height of the atmosphere, is governed by two temporal periods in which different exponential functions appear defining two different escape rates.The reciprocal value of the short-term escape rate provides an estimate of the average residence time of typical particles.The analogous quantity belonging to the long-term escape rate characterizes exceptional particles that become convected or remain in the free atmosphere for an extremely long time, respectively.It is interesting to note that the escape rates of particles of different sizes are found to vary in a broad range rather rapidly, roughly exponentially with the particle size.These investigations provide a Lagrangian foundation for the concept of deposition rates.

Fig. 1 .Fig. 1 .
Fig. 1.Proportion n/n 0 of the number of survivors in the three setups described in the legends.Initial conditions and particle radius: n 0 = 2.5•10 5 particles uniformly distributed over the globe (a) on p 0 = 500 hPa with r = 2 µm, (b) on p 0 = 500 hPa with r = 9 µm, (c) on p 0 = 700 hPa with r = 4 µm, (d) on p 0 = 900 hPa with r = 10 µm.Dashed lines in panel b illustrate the short-term and long-term escape rates.

Fig. 4 .
Fig. 4. Statistics of the residence time for different initial levels for n0 = 2.5 • 10 5 particles including t of rain and turbulent diffusion (setup 1).Minimums and maximums (lower and upper lines), upper an quartiles (boxes), medians (lines in the boxes) and averages (stars) are indicated.Dashed and dot indicate exponential fittings to median and average vs. r for given initial pressure levels.

Fig. 5 .Fig. 5 .
Fig. 5. Dispersion of a volcanic ash column of initial size 1 • × 1 • × 400 hPa from Mount Merapi consisting of n 0 = 2.16 • 10 5 particles with radius r = 10 µm initialized at 0 UTC on November 1, 2010.The initial position of the column center is λ 0 = 110.44• E, ϕ 0 = 7.54 • S, p 0 = 500 hPa, and particles are distributed uniformly in the column.Panel a, b, c, d, e illustrate the dispersion 1, 3, 5, 7 and 9 days after the eruption.Colorbar indicates the pressure level of the particles in hPa.Panel f shows the total precipitation of the last 12 hours at 0 UTC on November 8, 2010.Colorbar indicates here the total precipitation in mm/12 h.

Fig. 6 .Fig. 7 .
Fig.6.Proportion n/n0 of the number of survivors of the eruption in Fig.5vs.time in simulations includin the effect of rain (setup 1) (blue) and without rain (setup 2) (red) with exponential fittings to the data (bla dashed lines).The slope of the fitted curves are κs = 2.295 day −1 and κ = 2.066 day −1 for the blue curv and κs = 1.393 day −1 and κ = 1.739 day −1 for the red curve.

Fig. 7 .Fig. 7 .
Fig. 7. Vertical distribution of the proportion of the particles dispersed in Fig. 5 in vertical layers of size 50 hPa for 0, 1, 3, 5, 7 days after the hypothetic eruption.The dashed horizontal line represents the surface.

Fig. 8 .
Fig. 8. Dispersion of the volcanic ash containing r = 5 µm particles from the Mount Merapi eruption (all other parameters are the same as in Fig. 5).Panel a, b, c illustrate the dispersion 3, 7 and 20 days after the eruption.Colorbar indicates the pressure level of the particles in hPa.

Fig. 8 .
Fig. 8. Dispersion of the volcanic ash containing r = 5 µm particles from the Mount Merapi eruption (all other parameters are the same as in Fig. 5).(a, b, c) illustrate the dispersion 3, 7 and 20 days after the eruption.Colorbar indicates the pressure level of the particles in hPa.

Fig. 9 .Fig. 9 .
Fig.9.Quantity ∆n/(n 900 ∆t) vs. time in a simulation including the effect of rain and turb (setup 1) for n 0 = 2.5 • 10 5 particles initiated with uniform distribution on the 900 hPa level.∆ number of deposited particles over a time step ∆t, n 900 represents the number of particles below above the surface (p s = 1000 hPa).The lower (upper) curve belongs to particle size r = 1 µm (

Table 1 .
Short-term (κ s ) and long-term (κ ) escape rates for the cases represented in Fig.1.

Table 2 .
Short-term escape rate κs for a particle ensemble with vertical extension (described in the text), and the average κs and standard deviation σ of the short-term escape rates of different initial pressure levels.r[µm] κs [day −1 ] κs [day −1 ] σ [day −1 ] Statistics of the residence time for different initial levels for n 0 = 2.5 × 10 5 particles, including the effect of rain and turbulent diffusion (setup 1).Minimums and maximums (lower and upper lines), upper and lower quartiles (boxes), medians (lines in the boxes) and averages (stars) are indicated.Dashed and dotted lines indicate exponential fittings to median and average vs. r for given initial pressure levels.