Atmospheric Chemistry and Physics Discussions Interactive comment on “ Ship plume dispersion rates in convective boundary layers for chemistry models

Reviewer specific comment 1: There are indeed some errors in the formulation of Equation (6). The obvious one, as spotted by the referee, is that the time-like arguments of p1 should be t − t′ and not t′ . We would like to mention that we wrote Equation (6) only to show how obtain a 3D concentration field of a “generic puff” at time t from the scattered simulated particles positions. This recombined concentration field C(x, y, z, t) is defined in a coordinate x, y, z that can be related to the simulated particles positions xp(t), yp(t), zp(t) following:

. Supported by measurements, it has been shown that ignoring or misrepresenting chemical conversion inside the plume can lead to important overestimation of NO x and O 3 emissions in global models von Glasow, 2003;Davis et al., 2001). Those systematic biases especially arise from dilution process during the early stage of plume dispersion that is before the plume has 10 been sufficiently diluted throughout the boundary layer (Esler, 2003;Chen et al., 2005;Schlager et al., 2006;Song et al., 2003;Poppe et al., 1998).
Although such dilution effects on plume chemistry simulations has been widely considered in recent studies (e.g. Esler, 2003, and references therein), detailed description of the dispersion regimes in realistic boundary layer is still lacking, and chemical mod- 15 elers have to rely on parsed observations (e.g. von Glasow et al., 2003) or simple theoretical approaches using Gaussian plume models and homogeneous turbulence (Poppe et al., 1998). The present study relies on detailed and realistic plume dispersion data generated by a Lagrangian Particle Dispersion Model (LPDM) coupled to an atmospheric Large Eddy Scale (LES) model. The objective is two-fold: first is to Introduction at ship exhaust stack. The same simulations are finally used to determine a simple constant dilution rate suitable for subgrid-scale parameterization of effective emissions in coarse chemistry transport model; conclusions are given in Sect. 5.

Models and methodology
2.1 LPDM and LES models 5 A solution to simulate realistic ship plume dispersion is to represent it by a large number of passive particle trajectories in a Lagrangian framework using time-dependent velocity fields from LES model outputs. Since the pioneering work of Lamb (1978), this method has been successfully employed in various boundary-layer and plume cases (e.g. Mason, 1992;Gopalakrishnan and Avissar, 2000;Weil et al., 2004;Cai et al., 10 2006). Spatial and temporal evolutions of ship plumes are simulated using DIFPAR (Wendum, 1998), a stochastic LPDM based on a Markovian "zeroth order" equation. The code has been developed by Electricité de France R&D and adapted to the scale of LES and ship emissions issue by CERFACS. The model uses periodic wind, turbulence 15 and thermodynamic input fields which are linearly interpolated in space and time, to compute trajectories of passive particles. Those fields are provided by LES of marine boundary layers using the Non-Hydrostatic atmospheric model Meso-NH (Lafore et al., 1998) The heat release from ship exhaust stacks represents an additional buoyancy flux that controls plume dispersion, especially close to the source. In some cases, the combined effect of momentum and buoyancy may lead to plume rise of 2 to 10 times the actual release height (Arya, 1999), thus reducing the maximum ground concentration 5 by a factor up to 100 (Briggs, 1984). In turn, the structure of the atmospheric boundary layer (ABL) can strongly modify the rise and shape of the plume as in case of slightly decoupled boundary layer (e.g. Liu et al., 2000). However, an exact description of the effect of buoyancy forces on particle motion is intrinsically not feasible in a Lagrangian approach because of the nature of entrainment processes that require to take into ac-10 count all fluid-flow interactions simultaneously (i.e. the interactions between fluid particles and between particles and background flow). Consequently, most of Lagrangian formulations use separate plume rise models that include the effect of entrainment, by means of analytical or semi-empirical formula of bulk plume properties, such as those derived by Briggs (1975) or similar approaches. The vertical velocity generated 15 by these hybrid Lagrangian-Eulerian models is then added to the Lagrangian particle motion (Luhar and Bitter, 1992;Anfossi et al., 1993;Hurley and Physick, 1993;Hurley, 1999). In our approach, the plume rise scheme added to DIFPAR model is based on the idea of Anfossi et al. (1993). The scheme is a simple parameterization that uses a generalized form of Briggs plume rise formula (Anfossi, 1985): where H is plume height above stack, u and S are, respectively, the local wind module and the stability parameter, provided by the LES model; while F is the initial buoyancy where w 0 is the initial exhaust air velocity; r is the stack radius; and T f and T a are the exit plume and ambient temperatures. Note that the estimation of F is usually complicated, due to its dependence on ship engine power and sailing configuration, 5 background air temperature, and eventual additional exhaust additional devices. For diesel engine-equipped, regular cruising ocean-going vessels, estimated values of F range from 80 to 250 m 4 s −3 with a typical mean value around 120 m 4 s −3 (Hobbs et al., 2000;Pingkuan et al., 2006;Moldanova, 2007). The basic idea is to assume that each particle is by itself a plume that rises according 10 to Eq. (1) independently on the others. In order to mimic mutual particle interaction effects on buoyancy, and to fit the empirical plume radius increase near the stack (Anfossi et al., 1993), a normally distributed F i is assigned to each i th particle such that: where δ i is a Gaussian random number with zero mean and unit variance. Finally, an 15 additional vertical velocity w b due to buoyancy forces is added to the particle motion at each time step: The plume rise scheme is applied to each particle until w b is less than the local vertical turbulent fluctuations.

Simulations set-up
The Cloud-Topped Marine Boundary Layer observed during First ISCCP Regional Experiment NH, using initial conditions provided by the EUROCS model comparison exercise (Duynkerke et al., 2004). FIRE is a well documented experiment of a typical convective marine boundary layer, and a benchmark validation test for the atmospheric modellers community. It has been successfully simulated by Meso-NH in LES mode, for validation against measurements and within intercomparison studies (Sandu et al., 2005;Sandu, 5 2007). The simulation has a spatial resolution of 50×50×10 m 3 and horizontally cyclic domain size of 20×5×3 km 3 , large enough to prevent side effects and constraints on the development of convective cells. Atmospheric profiles exhibit a neutral boundary layer of about 600 m depth capped by a strong inversion ( Fig. 1 left panel), and moderate to strong ABL wind with no shear. This simulation will be hereafter named FIRE1 10 case. Similarly, two other simulations named FIRE2 and FIRE3 cases, are performed with the same set-up than before, but with modified initial profiles of liquid water potential temperature θ L and total water mixing ratio r tot , in order to modify the boundary layer height. To that end, θ L and r tot are kept constant in the boundary layer, but the temperature and humidity jump characterizing the inversion layer is set at the desired 15 altitude, above which the profiles follow the initial free troposphere slope ( Fig. 1 left panel). The resulting boundary layer heights are 400 m for FIRE2 case and 800 m for FIRE3 case. All FIRE simulation cases run for one hour after 2 h of spin-up, in order to reach a fully developed, steady, well-mixed boundary layer. In addition, a fourth kind of boundary layer is simulated based on the convective 20 steady state conditions observed during the Barbados Oceanographic and Meteorological Experiment (BOMEX; Holland and Rasmusson, 1973). This case served as a basis for a well known LES intercomparison study, and is fully described in Siebesma et al. (2003). In the present study we used the same domain size and resolution of FIRE cases, with a 2 h simulation run. The atmospheric profiles (see Fig. 1 right panel) 25 show a 600 m neutral boundary layer, capped by conditionally unstable mixing layer (weak inversion). The moderate to strong wind is constant in the ABL, with weak shear above.
Each of the four simulations is identify by the two key parameters that directly control Interactive Discussion plume dispersion: the boundary layer depth z i and the turnover time scale t * =z i w * , with w * representing the convective velocity scale defined here as in Eq. (5): where θ v is the horizontally averaged virtual potential temperature and wθ v is the mean vertical buoyancy flux. The values of these two parameters are presented in Table 1, 5 and are used as scaling factors to characterize plume dispersion.

Simulating a characteristic plume
One major issue in the simulation of ship plumes -e.g. for validation of parameterizations or comparison with observations -is that the shape of the plume strongly depends on the local flow conditions around each particle trajectory. For example, the rise and 10 fall of a plume are influenced by extremely scattered and random events such as the updrafts and downdrafts occurring in a turbulent convective boundary layer. Besides, the exact values of the ship speed and direction (relative to the wind) are generally unknown. All this uncertainties and the wide range of spatial and temporal variability of the parameters make any single simulation of the plume a particular case among 15 huge number of possible combination. On the other hand, simulation of all conceivable cases is obviously unaffordable. A reasonable compromise made in this study is to reconstruct the properties of a generic ("idealized") plume that is somehow representative of any plume evolving in a given ABL. This done in two steps as detailed next.

20
The first step is to consider that a plume is the superposition of puffs released at different time and different locations. The above mentioned generic plume model can then be treated using one single puff simulation. The second step is to model such a puff -characteristic of the whole boundary layer -so as to get rid of the dependency 6800 Introduction Interactive Discussion on local time-changing dynamical events. The idea is to release a large number of particles equally distributed in space and eventually in time at the same ship stack level over all the modeled boundary layer. The time evolution of the normalized concentration field of this "generic puff", characteristic of the given boundary layer, can then be recombined from the particle positions 5 at time t following (adapted from Lamb, 1978): where Q is the source strength function which can be time dependent;x p , y p and z p identify the particle position at time t; x 0 and y 0 identify the initial particle position at its release time t ′ ; U, V are horizontal mean wind speed components of the ABL. p 1 10 is the position PDF for particles found at time t in a coordinate system that depends on the particle initial position and the mean advection. The method is thus particularly suitable for practical implementation in atmospheric chemistry models which rely only on general boundary layer parameters. In our cases, the LES fields are stored every minute, provided off-line to the Lagrangian model and interpolated in time and space (taking advantage of the cyclic boundary conditions of the LES domain); 50 000 particles are released instantaneously, equally distributed horizontally over the domain, at an initial altitude of 60 m, matching the upper limit of ship height above sea level in international PANAMAX standard for major marine vessels. The mean initial buoyancy flux at ship stack F (Eq. 2) is a tunable parameter for the puff dispersion simulation. The plume rise scheme is applied to each particle following Eq.
(3). For each boundary layer case, seven plume simulations are performed with 5 a mean initial buoyancy flux F =0 m 4 s −3 (no buoyancy), 50, 100, 120, 150, 200 and 250 m 4 s −3 . Although the initial buoyancy flux can have a major impact on plume dispersion as discussed above, the heat and humidity release at ship stack can not significantly raise the in-plume temperature and water mixing ratio above background level, even as close 10 as 200 m from the stack. In rare and specific cases, only a small temperature increase of less than 0.4 • K has been documented (Hobbs et al., 2000). In the simulations, the temperature difference between plume and boundary layer background can be estimated using the plume rise scheme integrated in the Lagrangian model. For each particle and at each time step, a vertical buoyancy acceleration is obtained by taking time derivative of Eq. (4). This acceleration can be related to temperature difference between the plume air particle and the surrounding environment via a buoyancy force: where θ vref is the virtual potential temperature of environment, provided by the LES 20 fields. However, the heat flux at ship stack has little impact on the estimated plume temperature, which drops quickly to background value. This is shown in Fig. 2 that reports the maximum temperature excess with respect to background as a function of time and altitude for a initial buoyancy flux at ship stack of 250 m 4 s −3 ; almost twice our estimated mean value for a cruising large ocean-going vessel (Moldanova, 2007 This result is in agreement with in-situ measurements of temperature close to the ship stack. Although the buoyancy force is active only during early stage (a few minutes) of puff evolution, it indirectly affects the dispersion of the plume on time scales of the order of several convective turnover times of the boundary layer. This is shown in Fig. 3 that   5 reports the evolution of the normalized mean vertical concentration, defined as in Cai et al. (2006) by for various initial buoyancy fluxes in the BOMEX case. In this case, the high sensitivity of the simulations to this parameter can be explained by the fact that the boundary 10 layer is capped by a weakly stable inversion layer. Thus, if particles get sufficient buoyancy or are driven by strong convective burst, they can overshoot the boundary, and eventually stay trapped in the upper layer (Liu et al., 2000), as clearly shown in Fig. 3 for F =250 m 4 s −3 . This possibility is prevented in the FIRE simulations due to the strong capping inversion, so that all particles remain inside the boundary layer. 15 However, the initial buoyancy flux has again a strong impact on the plume dilution pattern. In fact, up to three convective turnover times are needed to reach a final state independent on F (not shown). Note that differences are not significant for values larger than 100 m 4 s −3 .

Influence of boundary layer height and turnover time scale on plume dispersion 20
The height of the boundary layer has two effects on the evolution of the plume and its shape. The first (direct) effect is that it sets a global length scale for dilution -and defines different regimes of dilution in the case of convective boundary layers. As shown in Fig. 4, a typical exhaust buoyant puff first disperses quickly into the boundary layer until its upper portion reaches the inversion height. As a result, the dilution process Introduction Interactive Discussion then slows down. On the other hand, the lower portion of the puff falls down and touches the ground soon after, which again reduces the dilution rate of the mean puff concentration. Note that there are situations, depending on the initial buoyancy flux, the height of exhaust release and the inversion height, where the first two dilution regimes cannot be clearly separated. Once the puff is roughly well mixed in the vertical 5 direction, its dilution is mainly driven by horizontal dispersion. The second (indirect) effect of boundary layer depth is that the inversion height controls the size and turnover time of convective cells and then impacts the dispersion process of the plume. Figure 5 shows the normalized mean vertical concentration for the FIRE1, FIRE2, FIRE3 and BOMEX cases. In the FIRE2 case, with t * ≈12 min, the 10 puff dispersion in the vertical direction is driven by the boundary layer convective pulsation. This feature is also present in the BOMEX case, although the inversion is less permeable (see the discussion in previous section). This pattern is not clearly seen in the FIRE1 case where t * ≈22 min, i.e. only one third of the total simulation run-time, and much larger than the onset of the third dilution regime. That is to say that dilution 15 is too fast compared to larger-scale turbulent convection, whose action is outweighed by small-scale dispersion. (Note that, for the same reason, in the FIRE3 case there are no initial "waving" patterns). For all simulated cases, the second indirect effect becomes negligible after a few convective turnover time scales, when the plume can be considered as "well mixed".

Estimate of the plume dilution rate
In chemical box models (that use an explicit plume model, see e.g. von Glasow et al., 2003) or in chemistry transport models (where the plume is parameterized at subgrid scale level), mixing in generally represented using two air masses where chemical species are assumed to be homogeneously distributed: the well mixed background 25 (e.g. the boundary layer) and the plume itself. The dilution process and the way it influences the plume chemistry is then controlled by the exchange rate between these ACPD 8,2008 Ship plume dispersion rates in convective BL where C p is the mean in-plume concentration of a given chemical species (or aerosol species); C a is the mean homogeneous background concentration, while F (t) is the (unknown) dilution rate to be determined.

Dilution rate for chemical box models
In our model, the plume is considered as a superposition of puffs, and mixing occurs only at edge of the puff via entrainment of background air in the crosswind plane (see e.g. von Glasow et al., 2003). Thus, we neglect the puff dilution in the downwind direction (the so-called "slender plume approximation"), which is fully justified in our 10 simulations since wind is not nil and we do not consider dispersion very close to the source. In this case the change in concentration within the plume is related to the change in puff surface A p on the cross-wind plane such that: Three different definitions of the puff surface have been examined depending on the 15 representation of the puff/background interface: the first one consists of the surface that contains 98% of the emitted particles; the second one takes the core of the plume that contains 60% of the particles; and the last one defines the surface as: where σ y and σ z are the dispersion parameters in crosswind and vertical directions, 20 respectively ( Dosio et al., 2003). However, we did not find significant differences in the results (not shown) using these three definitions, thus we retained the last one, Eq. (11) as it provides smooth variation of plume area, and is more easily related to the widely used Gaussian plume models in the literature. Figure 6 presents the natural logarithm of the plume crosswind surface for all simulations as a function of non-dimensional time t/t * .

5
This scaling eases comparisons among boundary layers with their own convective time scale t * , which is a common parameter characterizing the convective ABL in generic dispersion studies (see e.g. Fedorovich, 2004).
The dilution rate function F (t) is obtained from crosswind surface data using Eq. (10) and is presented in Fig. 7 for all simulations. The figure suggests that dispersion in the 10 boundary layer is characterized by the competition between the buoyant plume dilution regimes and convection processes, as discussed above. This results in abrupt changes in the derived instantaneous dilution rates that is clearly visible in the logarithmic scale, especially for FIRE1 case (note that data in Fig. 6 were not smoothed out for the derivation). Meanwhile note that the computed dilution rates were found to be 15 almost insensitive to the initial buoyancy flux, despite the rather scattered nature of the plume vertical dispersion patterns shown in the previous section.
One of the goals of this paper is to provide the characteristic dilution rate for convective BL that can be used in a simple and straightforward way into chemistry transport models. To that end, we propose a power-law function that best fits the simulation 20 results: Table 2 reports the values of the fitting parameters a and b, while Fig. 7 presents the evolution of best-fitted dilution function rate in Eq. (1), together with the curves extracted from the simulations for various initial buoyancy fluxes. Note that derived dilution rates F (t) are in minutes. Both Table 2 and Fig. 7 confirm the minor effect of initial buoyancy flux on plume dilution rate estimates. In chemical plume models, such those included in meso-scale, regional or global circulation models, the plume is usually parameterized as a subgrid-scale process. Indeed, a detailed representation of the plume dispersion would require grid resolutions that are not affordable by present supercomputers. One must find out a more simple de-5 scription of the dilution rate function F (t) that does not change within a characteristic time smaller than the time step of such models (or simply independent on time). Moreover, in those models, the plume is not considered as a superimposition of puffs, but rather a single puff which diffuse slowly within the grid cell. We then look for a plume dilution rate in a given boundary layer in the form:

ACPD
where the characteristic dilution time scale τ is found by fitting the mean in-puff particle concentration derived from a single puff simulation, for each initial buoyancy flux. This dilution rate should not depend of the initial setup of the puff diffusion and the early stage of the dispersion when buoyancy is still dominant and the puff has not yet reached 15 a steady evolution (i.e. the first two dilution regimes described in Sect. 3).
To do so, we linearly fit the logarithm of the normalized mean in-puff particle concentration evolution starting at a time equal to the convective turnover time t * . This is justified from the linear behaviour of this mean concentration, as seen in the right panel of Fig. 4. The absolute value of the line slope thus gives the characteristic dilution rate 20 1 τ. Figure 8 presents the characteristic dilution time τ scaled by t * as a function of the initial mean buoyancy flux, for all simulation cases.
The figure indicates the same evolution for all boundary layers situation: initially τ decreases as F increases, and quickly converges to a constant value for F >100 m 4 s −3 . This is due to the fact that initial buoyancy flux act as an accelerator of plume dilution until it reaches a "well-mixed" state throughout the boundary layer, that is when dilution relies mainly upon convective related processes. In other words, for small F , this state 6807 Introduction can be not fully reached within a period t * , as discussed in Sect. 3. This clearly shows the importance of using more sophisticated methods to represent dilution in boundary layers, compared to standard Gaussian plume models, where plume rise is either totally neglected or taken into account by simply lifting the initial emission source. Figure 8 also shows that, for initial buoyancy fluxes of typical ocean-cruising vessels 5 values (F =80 to 250 m 4 s −3 ), all scaled dilution times exhibit similar values of about τ t * ≈4.12±0.47. This interesting result suggests a general (far-field) dispersion law since the data shown in Fig. 8 come from very different boundary layers (especially BOMEX versus FIRE cases), and very different plume dilution patterns. Furthermore, the range of ship stack buoyancy flux considered in this study also suggests that vessel 10 type or size has a minor impact on plume dispersion. Nevertheless, this has to be confirmed by additional simulations and observations (for example using the simulation results and tank experiments reported in Fedorovitch, 2004), and/or by theoretical work to derive a more systematic parameterization.
The result is however encouraging, as it provides a simple parameterization of dilu-15 tion in chemical transport models. Indeed, the boundary layer height z i , the surface convective velocity w * , and then t * , are directly available from standard outputs of most atmospheric models. Note that this is only justified in convective boundary layer, which is a common situation over ocean surfaces.
The proposed parameterization should be extended to cover the initial dilution 20 regime (the dilution time scale is effectively constant only after the plume gets mixed with background). In such a parameterization, one should take into account dispersion processes from the source to this steady dilution regime, using, for example: (i) the detailed results of dispersion presented in previous section (but the computational gain of a parameterization would be lost); (ii) instantaneous release of ship exhaust 25 injected directly at well-mixed stage throughout the BL (but bias are expected in net chemical products); (iii) using equivalent or effective emissions (Esler, 2003;Paoli et al., 2008 1 ;Cariolle, 2007)  Interactive Discussion stage (additional CPU cost is then only limited by prior off-line computations).

Summary and conclusion
The present paper was motivated by the increasing concern about the environmental impact of ship traffic emissions. Although it is now recognized that early stage of plume dispersion has a strong impact on plume chemistry, chemistry models generally rely 5 on rough parameterizations of plume dispersion or simple Gaussian plume models. The aim of this study was to provide a realistic description of the entrainment-mixing processes between plume and background air that can be easily incorporated into chemical box models or large scale chemistry transport models. To that end, detailed ship plume simulations were performed within various convective boundary layers, us-10 ing a Lagrangian Particle Dispersion Model (LPDM) driven by high resolution atmospheric Large-Eddy Simulations (LES). The analysis focused on early stage of plume dispersion, i.e. the first one-two hours after release. The LPDM contains a plume rise scheme to improve the plume representation and explore the impact of the initial buoyancy flux at ship exhaust funnel on plume dispersion. 15 The different plume dilution regimes up to the well-mixed plume steady state were clearly identified in the simulations. The impact of plume rise scheme, boundary layer height and characteristic convective turnover time scale on plume dispersion was also discussed. Analytical dilution rate functions were derived as they are classically used in chemical box models. These rates were found to depend on the convective turnover 20 time scale of the boundary layer, but are almost insensitive to the initial buoyancy flux, despite its striking effect on plume dispersion patterns.
Finally, a subgrid plume dispersion model for chemistry transport models was proposed. The plume dispersion was parameterized by a constant dilution rate characteristic of the well-mixed dilution regime. This parameter was found to be weakly 25 tions and emissions from concentrated sources into global models, C. R. Mechanique, submitted, 2008.  Lamb, R. G.: A numerical simulation of dispersion from an elevated point source in the convective planetary boundary layer, Atmos. Environ., 12, 1297-1304, 1978. Liu, Q., Kogan, Y. L., Lilly, D. K., Johnson, D. W., Innis, G. E., Durkee, P. A., and Nielsen, K. E.: Modeling of Ship effluent transport and its sensitivity to boundary layer structure, J. Atmos. Sci., 57, 2779Sci., 57, -2791Sci., 57, , 2000