Parameterization of subgrid aircraft emission plumes for use in large-scale atmospheric simulations

Parameterization of subgrid aircraft emission plumes for use in large-scale atmospheric simulations A. D. Naiman, S. K. Lele, J. T. Wilkerson, and M. Z. Jacobson Aeronautics and Astronautics, Stanford University, Stanford, CA, USA Civil and Environmental Engineering, Stanford University, Stanford, CA, USA Received: 22 October 2009 – Accepted: 5 November 2009 – Published: 18 November 2009 Correspondence to: A. D. Naiman (anaiman@stanford.edu) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Commercial aircraft consume most of their fuel at a cruise altitude near the tropopause (Penner et al., 1999), where exhaust has a long residence time and conditions are conducive to the formation of line-shaped condensation trails (contrails).Under certain conditions, contrails persist, spread, and merge into fields of cirrus clouds.Even if contrails do not persist, aging aircraft exhaust can lead to formation of cirrus clouds (Minnis et al., 1998).The effect of contrails and aviation-induced cloudiness on climate is highly uncertain.The most recent IPCC estimate of the combined radiative forcing of linear contrails and induced cloudiness was 30 mW/m 2 with an uncertainty interval ranging from 10 to 80 mW/m 2 (Forster et al., 2007).This estimate also lists the level of scientific understanding as "low" for linear contrails and as "very low" for induced cloudiness.The range of estimates for these effects is large, and if the upper end of  et al., 2009).
Most studies of the climate effects of aircraft have been limited to estimating globalscale radiative forcing due to them (e.g., Minnis et al., 1999;Myhre and Stordal, 2001).Minnis et al. (2004) matched trends in cirrus coverage to atmospheric variables and aviation fuel inventory for historical data from 1971 to 1995.Ponater et al. (2002) and Marquart et al. (2003) parameterized potential visible contrail coverage as a function of atmospheric conditions and air traffic density.The resulting distribution was applied as a diagnostic model of contrail coverage in a global climate model.Additional tuning of the parameterization was required to match the distribution to contrail observations over a specific region.The Forster et al. (2007) estimate of radiative impact is based on Sausen et al. (2005), which updates previous studies such as these to account for increased traffic.
As noted by Burkhardt and K ärcher (2009), these estimates of radiative forcing suffer from certain conceptual limitations.Since these models combine potential contrail coverage with air traffic density, they are suited to predicting the formation of contrails.
They do not, however, provide prognostic information about how long these contrails persist or how atmospheric dynamics affect their properties.They also suffer from systematic error introduced by tuning global contrail coverage based on observations from a small region.Burkhardt and K ärcher (2009) introduced an alternative parameterization of contrail cirrus coverage, contrail length, and grid mean ice water mass mixing ratio based on physical processes.Contrail formation, transport, spreading, deposition, and precipitation processes were applied to these variables, which were tracked on the grid scale.The study emphasized the role of transport in distributing contrail coverage and found that contrail coverage scales more closely with supersaturation than with formation frequency.
One limitation of previous work is the treatment of aircraft exhaust plumes and contrails at the global grid scale.This requires parameterizations of the physical processes listed above, which for contrails take place on scales considerably smaller than the grid.For example, Burkhardt and K ärcher (2009) model (Roeckner et al., 1996), used a parameterization to determine deposition rates to contrail ice particles, as these clouds exist on a scale of 10 km compared to a grid spacing of 270 km.They also parameterized the effect of shear on contrail spreading, specifying both a contrail vertical thickness and a spreading constant.Tests in the study did not examine the sensitivity to the specified thickness, but did show significant sensitivity to the spreading parameter.
The evolution and dispersion of aircraft exhaust plumes are physically subgrid scale processes.We therefore take a Lagrangian approach in tracking individual aircraft exhaust plumes until they have dispersed, either by growing to grid scale or by diluting to background concentrations.Using a new aircraft emissions inventory that gives individual flight trajectories over a year (Wilkerson et al., 2009) and a global climate model that models aerosol processes on the individual contrail level (Jacobson et al., 2009), the parameterizations noted above are exchanged for physical models.In the case of the plume transport and spreading processes, this treatment requires a model that can assess the evolution of the plume location, volume, and shape based on grid scale variables.
This paper presents a model of aircraft plume dynamics that is intended to fulfill this role in a large scale atmospheric simulation.It provides prognostic equations for the advancement of the volume and width of a plume based on variables provided by the atmospheric simulation on the grid scale.Although the equations are simple, comparison with a high fidelity model of plume dispersion shows that they are adequate to describe plume dynamics as compared to the level of fidelity of a large scale atmospheric simulation.
Section 2 presents a new model of emission plume dynamics that uses global grid scale variables to advance the plume.Section 3 compares the new model to a large eddy simulation of an aircraft exhaust plume and to the analytical solution to the dynamics of a sheared Gaussian plume.Finally, Sect. 4 argues that the new model also provides a reasonable model of line-shaped contrail dynamics and suggests how it might be used as such in a global climate model.

Subgrid Plume Model
In this section, we describe our approach to developing the Subgrid Plume Model (SPM).We describe the basic approximations of the model, derive generalized equations, and then apply them to a specific parameterized plume shape.The resulting model equations can be integrated to analytical solutions with additional minor assumptions, as presented in Sect.4.3.

Approach
Commercial jet aircraft in cruise emit exhaust that spreads at a much lower rate than the flight speed, resulting in long, slender plumes along its flight path.Our approach in the SPM is to idealize these plumes as high aspect ratio, linear structures with a crosssection that may vary temporally and spatially.Thus, if the aircraft flight trajectory is split into segments, we can represent a corresponding segment of exhaust plume with a length and a cross-section.
The SPM has three goals: to capture the important aspects of aircraft emissions plume evolution under a variety of conditions, to maintain physical properties by obeying conservation laws, and to provide an accurate model of plumes at low computational cost.This third goal is important, as it is the property that allows the SPM to be useful as a subgrid process to track tens of millions of flights during a simulated year in a global model.
To this end, several approximations are made before deriving the SPM equations of motion.First, the SPM treats plumes as tracers of the atmospheric fluid without internal dynamics.Second, the atmospheric disturbances that affect the plume are aggregated into three processes that capture major modes of development: advection due to mean wind, distortion due to wind shear, and dilution due to turbulent mixing.Third, these processes are treated as though they are decoupled.The implications of these approximations with respect to contrails will be examined further in Sect.4.2.The three processes treated by the model are applied to the plume segment repre-Introduction

Conclusions References
Tables Figures

Back Close
Full sentation.Advection moves the endpoints of the segment, changing its length.Shear and diffusion act in the plane perpendicular to the segment (due to the slender plume approximation), changing its cross-sectional shape.Quantities needed by the climate model can then be calculated -for example, plume volume is calculated by simply multiplying cross-sectional area by length.This is the approach taken in deriving the SPM equations of motion.

General equations of motion
Figure 1 shows a plume cross-section at a segment endpoint defined by the position vector, x, in the global reference frame.The change in the position of the segment endpoint over time is: where u i is the mean wind component in the i direction.This equation applies advection to the plume segment.
Figure 1 also shows a relative position vector, ξ, in a reference frame with its origin at x and the same orientation as the global frame.The change in the position of the piece of plume cross-section at ξ is, to first order, (2) This equation represents the kinematic deformation of a material element of the plume and applies the effect of wind shear to the plume cross-section.
The effect of mixing on the plume cross-section is applied using a one-dimensional diffusion equation.The change in the position of the piece of plume cross-section at ξ is: where D i is the diffusion coefficient in the i direction.This equation is meant to indicate that the change in each component of ξ is related to the diffusion coefficient in that direction, so there is no sum over i on the right hand side of Eq. 3.

SPM equations
To apply these equations to a model plume, we specify a particular cross-section.The SPM uses an ellipse with three degrees of freedom (two radii and a rotational angle).This choice reflects the results of studies of plumes and contrails at late times under turbulent and shear conditions (D ürbeck and Gerz, 1996;Schumann et al., 1998;Chlond, 1998;Lewellen and Lewellen, 2001;Huebsch and Lewellen, 2006).In particular, Schumann et al. (1995) fit a two-dimensional Gaussian plume to their observations of aircraft exhaust plume cross-sections.The SPM is compared to the analytical solution used in that study in Sect.3. Since the plume is represented by a high-aspect-ratio segment, shear and diffusion are limited to act on the cross-section only in the plane perpendicular to the segment (the slender plume approximation).Figure 2 shows a schematic of the plume crosssection.The z-axis is the same as the global coordinate, while a new horizontal s coordinate is defined as orthogonal to the z-axis and in the plane of the cross-section.This figure also defines the three degrees of freedom of the cross-section: a, the initially vertical radius of the ellipse; b, the initially horizontal radius of the ellipse; and θ, the rotational angle of the ellipse.The angle, θ, is defined as the clockwise angle between the z-axis and a and is initially zero.
At cruise altitudes, vertical wind shear dominates the other terms in the velocity gradient tensor, ∂u i /∂x j .Equation (2) can therefore be written for a: θ geometrically.Similar manipulations using Eq. 2 result in ODEs for the magnitudes of the radii, a and b.The set of ODEs describing the effect of shear is: Note that these ODEs conserve the cross-sectional area of the plume ellipse.Diffusion coefficients in the SPM are estimated in the a and b directions as where D v and D h are the vertical and horizontal diffusion coefficients, respectively.Applying Eq. ( 3), the set of ODEs describing the effect of mixing are: In summary, Eqs. ( 1), ( 5)-( 7), and ( 8)-( 9) are advanced in the SPM to determine the location, volume, and shape of a plume segment over time.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion compare it to these previous studies of exhaust plumes and show that it produces very similar results.Instead of a Gaussian distribution, the SPM plume represents a quasi-uniform distribution of exhaust.Observations and simulations of contrails (e.g., Huebsch and Lewellen, 2006) show that real plumes are somewhere in between these two distributions.In Sect.4.2 we will further discuss the differences between exhaust plumes and contrails and the significance of the differences between the SPM and the Konopka analytical solution.D ürbeck and Gerz (1996) conducted a three-dimensional large eddy simulation (LES) study of late (wake-free) aircraft exhaust plume development.The goal of the study was to determine effective diffusion coefficients by fitting two-dimensional Gaussians to plume cross-sections computed using LES at different levels of background turbulence and shear.The comparison to the current SPM is thus a simple matter of advancing the SPM equations using the prescribed shear and calculated diffusion parameters from the study.
Figure 3 compares the SPM to the results from the D ürbeck and Gerz study and to the analytical solution to the diffusion of a Gaussian plume in a uniform shear flow (Konopka, 1995).The quanities plotted are σ 2.2σ v = acosθcosr z − bsinθsinr z , (13) Figure 3 shows results for the baseline case (case 2 in D ürbeck and Gerz (1996)) and d u x /d z = 0.003 s −1 .Both the SPM and Konopka analytical solutions match the LES results closely in the growth rates of horizontal and skewed variance.Neither provides a very good match for the growth rate of vertical variance, which, as noted by D ürbeck and Gerz, evolved linearly in the LES over an intial time and a late time, but at two different rates.This quantity also showed the most variation between LES plumes (three plumes were simulated for each case to check the effect of the variability of mixing processes within the flow).Figure 4 compares contours of exhaust concentration from the D ürbeck and Gerz LES to the SPM and Konopka solutions under the same conditions as noted above.The solutions for all three models are plotted at three times during the seventy minute simulation.Both the SPM and the Konopka model capture the spread of the plume as it is sheared and diffused.These plots illustrate that the difference in the vertical variance development between the models has little effect on the plume extent and horizontal spreading, which is dominated by the vertical shear.
Figure 5  cases, whereas the Konopka solution says that the plume area increases as a square root function under low shear and a quadratic function under high shear.The LES results showed some of this sensitivity to shear, but were generally more linear than the Konopka solution would suggest.Over the time scale presented by D ürbeck and Gerz, the SPM more closely matches the LES results quantitatively for the high shear case, with the plume area increasing by a factor of eight over the seventy minute simulation.

Using the Subgrid Plume Model
The SPM is intended to be used as a subgrid model of aircraft exhaust dynamics, including line-shaped contrails.In this section, we argue that the plume model is suitable for modeling aircraft exhaust and contrail development in the context of a global climate model.We also give an example of how it might be used in such a global model, including analytical solutions to the model equations.

Climate model interface
In order to discuss the implications of the approximations made in deriving the SPM, it is necessary to refer to the usage of the SPM in a particular large-scale atmospheric model, though these arguments could be made for other models of similar scale.The atmospheric model that the SPM was designed for is GATOR-GCMOM, a nested global-to-regional climate model that treats time-dependent gas, aerosol, radiative, dynamical, cloud, land, and ocean processes (Jacobson, 2001a(Jacobson, ,b, 2002(Jacobson, , 2003(Jacobson, , 2004)).The model is being used to investigate the global impact of aviation on climate (Jacobson et al., 2009) using an emission inventory that specifies individual flight trajectories (Wilkerson et al., 2009).The climate simulation tracks the volume and shape of emission plumes from individual aircraft over time using SPM segments.For each segment, the simulation uses plume volume to calculate the dilution of plume components for a microphysical model Introduction

Conclusions References
Tables Figures

Back Close
Full and optical property calculation.Segments are tracked individually until they are diluted, at which point their aerosol and water vapor components are added to the grid scale distributions.The optical properties and shape of each plume segment are used in calculations of radiative transfer through plume-occupied portions of the climate simulation grid cells.Details of these calculations are described by Jacobson et al. (2009).
When a SPM segment is added to the climate simulation, its initial cross-section is specified based on an estimate of plume sizes at the end of vortex descent from the literature (e.g., Lewellen and Lewellen, 2001).In the first implementation of the SPM, the values a = 120 m, b = 65 m, and θ = 0 are taken as typical for all aircraft in the simulation.This implementation will be refined.An ongoing LES study is developing a database of contrail sensitivities to atmospheric and aircraft parameters.The results from the study will be used to adjust the initial conditions in future implementations.

Applicability to contrails
In Sect.2.1, we noted several approximations that we used to reduce the subgrid plume model to a more tractable problem.We first treat the plume as a tracer of atmospheric dynamics.Even for a passive plume, this notably neglects the vortex dynamics found in the wake of an aircraft.The vortex wake descends until the Crow instability causes it to break up, generally within two minutes after the aircraft passes (Crow and Bate, 1976).This initial descent is essential in spreading the exhaust plume vertically, since the thermal buoyancy of the jet exhaust causes some to detrain from the vortices, leaving a vertical curtain of exhaust (e.g., Lewellen and Lewellen, 2001).The time scale of this vortex descent, however, is much smaller than a single time step of the atmospheric simulation (60 min) and can be accounted for simply by taking this initial descent into account when initializing the plume size (see Sect. 4.1).
In the case of a plume that is not passive, i.e., a contrail that contains significant ice particle density, the tracer approximation also neglects other effects.Line contrail ice particles have been observed to grow to effective particle diameters of 20-30µm within several hours of contrail formation, with induced cirrus particles up to 200µm observed 24766 Introduction

Conclusions References
Tables Figures

Back Close
Full (e.g., Minnis et al., 1998;Heymsfield et al., 1998).This full range of particles has sedimentation velocities from 1-10 cm/s (Seinfeld and Pandis, 1998).Sedimentation therefore removes large ice particles from the contrail or contrail-induced cirrus core and increases the vertical extent of the plume directly.In the presence of vertical wind shear, horizontal spread is dominated by the vertical extent of the plume, so sedimentation (to the extent that it occurs) also increases the spread of the plume in the horizontal direction.
The importance of this effect with respect to the climate model, however, is small enough to neglect.The observations noted above and simulations of line contrails and induced cirrus such as Unterstrasser and Gierens (2009) indicate that the number of particles on the smaller end of the size distribution is much higher than the number of large particles.Further, over the lifetime of a long-lived contrail (e.g., five hours), an average particle with a diameter of 30µm falling 3 cm/s would descend approximately 500 m, which is the vertical resolution of the climate model.If larger particles, which fall faster, form in the line contrail plume of the climate model, they are too few in number to have a large effect on the radiative balance (see for example, Unterstrasser and Gierens (2009) Fig. 1) and the spread of these particles away from the plume is therefore insignificant.In the context of GATOR-GCMOM, these particles can still grow to large sizes while contained in the SPM plume segment if the line contrail is sufficiently long-lived.In that model, when line contrails have dispersed, their vapor and aerosol particles are added to the grid scale, where they affect cirrus.Once the plume aerosol components are added to the grid scale, they will be tracked with the grid scale aerosol and cloud model, where sedimentation is accounted for explicitly as a function of particle size and composition.
Another approximation made in deriving the SPM was to aggregate atmospheric disturbances into three processes (advection, vertical shear, and turbulent mixing) and apply them uniformly to a plume segment.In GATOR-GCMOM, as in other global models, subgrid scale variations in parameters that drive these processes in the SPM are captured by diffusion coefficients, which are determined from grid scale shear and Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion stability characteristics.One area of interest is the transition from contrail clouds to cirrus.This transition implies loss of the characteristic linear shape of contrails, which is not accounted for by the SPM.This is related to subgrid perturbations of the velocities that cause vertical shear and turbulence.This transition is accounted for by increasing the diffusion coefficients passed to the the SPM to the levels typically experienced on the aircraft exhaust plume scale (Schumann et al., 1995) as described in Jacobson et al. (2009).Once the contents of the SPM reach a certain level of dilution relative to ambient atmosphere, they are added to the grid scale, where they can trigger cirrus cloud formation.
Figures 6 and 7 show a final comparison of the SPM and Konopka solutions under the case 1 conditions noted in Sect.3, advanced over ten hours of simulation time.
The two models again match each other closely in terms of plume shape and spread in the horizontal and vertical directions.After ten hours, the SPM predicts the plume area has grown by a factor of 28 compared to the initial area, whereas the Konopka solution predicts a factor of 24 growth.This difference is relatively insignificant in the context of the global atmospheric model.

Solutions to SPM equations
This section contains numerical solutions to the ODEs presented in 2.3 that have been applied inside a global climate simulation (Jacobson et al., 2009).The climate simulation treats the dilution of, microphysics in, and radiation through subgrid line contrail plumes from flights worldwide, but does not currently treat the advection of line contrail position.Although advection of line contrails is not directly treated, contrail components are advected once added to the grid scale following contrail dispersion, where they can induce cirrus cloud formation.In principle, Eq. ( 1) can be solved using a variety of techniques to directly advect contrail line segments, but the solution to this advection equation should be matched to the model time stepping scheme in which it is implemented.The remaining SPM equations are solved using the method of operator splitting.First, Eq. ( 5) is solved analytically by assuming a constant s = ∂u s /∂z over the time step: where θ 0 is the value of θ at the beginning of the time step.Next, Eqs. ( 6) and ( 7) are also solved analytically by again holding s = ∂u s /∂z constant and using Eq. ( 15): where a 0 and b 0 are the values of a and b at the beginning of the time step.These equations are used in the following discrete form: where superscripts refer to the time level and ∆t is the time step.These intermediate solutions, ã and b, are then used in the analytic solution of Eqs.
8 and 9: where D a and D b are calculated as noted in Sect.2.3, using an average value of θ = (θ n+1 + θ n )/2 over the time step.
The quantities used by GATOR-GCMOM can be calculated using these values.
A n+1 p = πa n+1 b n+1 , (23) where A p is the cross-sectional area of the plume, W p is the projected width of the cross-sectional ellipse, and Multiplying A p by the length of the plume segment gives the plume volume and multiplying W p by the length of the plume segment gives the top-view area of the plume, or the projected area of the plume onto the ground.

Conclusions
We have presented a new model of aircraft exhaust plume dynamics that is intended to be used as a subgrid scale model within a large scale atmospheric simulation.It provides prognostic equations for the evolution of individual exhaust plumes based on parameters provided by the large scale simulation.Although the equations and their analytical solution are simple to implement, the model shows good agreement with the results of high fidelity, three-dimensional simulations of exhaust plume development.
The model presented has been used as a plume dilution model within a large scale atmospheric simulation, which simulates the evolution of emissions from individual aircraft.Specifically, the dilution has been applied to the aerosol microphysical model to compute the formation and persistence of line-shaped aircraft contrails with in the 24770 where a = a s ŝ+a z ẑ and u s is the projection of the global velocity u onto ŝ.The ordinary differential equation (ODE) for θ is derived by relating the components of a to the angle Introduction horizontal, skewed, and vertical variances of the plume respectively.For the D ürbeck and Gerz plots, these quantities were calculated from the LES results and presented in the reference.Konopka gives analytical expressions for these quantities given a constant vertical shear and diffusion tensor.For comparison, the SPM quantities a, b, and θ have been converted to effective plume variances by projecting the plume ellipse onto the horizontal and vertical axes and solving for the skewed variance (as in Schumann et al. h = asinθcosr x + bcosθsinr x , shows the evolution of the plume cross-sectional area normalized by the initial plume area for two of the D ürbeck and Gerz (1996) cases.Case 1 has d u x /d z = 0.001 s −1 and case 4 has d u x /d z = 0.007 s −1 , with other parameters as noted for case 2 above.Both the SPM and Konopka analytical solution are plotted for comparison to the D ürbeck and Gerz results (their Fig. 8).Both solutions match the general trend of the LES results, with higher vertical shear causing faster area increase.The SPM solution displays nearly linear behavior (with a small quadratic component) in both Introduction 21)b n+1 = b2 + 2D b ∆t,

Fig. 1 .
Fig. 1.Schematic of a generalized plume cross-section at a segment end point in the global reference frame.

Fig. 1 .Fig. 2 .
Fig. 1.Schematic of a generalized plume cross-section at a segment end point in the global reference frame.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Horizontal, skewed, and vertical plume variances compared from a computational study and two analytical solutions: D ürbeck and Gerz (1996) LES results (left, their figure 7), the SPM equations (right, solid lines), and the Konopka (1995) solution to the diffusion of a Gaussian plume in a uniform shear flow (right, dashed lines).SPM quantities have been converted to plume variances for comparison.

Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 5. Normalized plume area compared from a computational study and two analytical solutions: D ürbeck and Gerz (1996) LES results (left, their figure 8), the SPM equations (right, lines with symbols), and the Konopka (1995) solution (right, lines with no symbols).The case legend on the right plots matches the original case legend from D ürbeck and Gerz (1996).

Fig. 7 .
Fig. 7. Comparison of the areas predicted by the SPM and Konopka skewed Gaussian solutions over ten hours of simulation time.
, in the framework of the ECHAM4 climate Figures In the future, this model could be used to model other processes, such as dilution for plume chemistry calculations, and for other stationary and moving point sources of emissions.The simplicity of this model makes it a good candidate for such future work, with a low computational cost that would allow it to be used in the Lagrangian tracking of exhaust plumes from many such sources within a large