Ice crystal number is a critical ingredient in the potential climate impact of persistent contrails and contrail-induced cirrus. We perform an extensive set of large-eddy simulations (LES) of ice nucleation and growth within aircraft exhaust jets with an emphasis on assessing the importance of detailed plume mixing on the effective ice-number emission index (EIiceno) produced for different conditions. Parameter variations considered include ambient temperature, pressure, and humidity; initial aerosol origin (exhaust or ambient), number, and properties; and aircraft engine size. The LES are performed in a temporal representation with binned microphysics including the basics of activation of underlying aerosol, droplet growth, and freezing. We find that a box-model approach reproduces EIiceno from LES well for sufficiently low aerosol numbers or when crystal production is predominantly on ambient aerosol. For larger exhaust aerosol number the box model generally overestimates EIiceno and can underestimate the fraction from ultrafine aerosol. The effects of different parameters on EIiceno can largely be understood with simpler analytic models that are formulated in low and high aerosol-number limits. The simulations highlight the potential importance of “cold” contrails, ambient ultrafine aerosols, crystal loss due to competition between different-sized crystals, and limitations on reducing EIiceno. We find EIiceno insensitive to engine size for lower aerosol numbers, but decreasing with increasing engine size for higher aerosol numbers. Temporal versus spatial representations for jet LES are compared in an appendix.
Given projections for significant increases in air traffic and concerns about its potential effect on Earth’s climate, aircraft contrails and contrail-induced cirrus have been the subject of many observational, modeling, and simulation studies in recent years [see, e.g., Kärcher (2018), Schumann and Heymsfield (2017), and Paoli and Shariff (2016) for recent reviews]. The water vapor and aerosol in aircraft exhausts can nucleate ice crystals that for supersaturated ambient conditions can grow and spread over time to form natural-looking cirrus in skies that otherwise would have remained cloud free. The number of ice crystals per length of contrail Nice plays a particularly critical role. In prior work based on simulations and scaling analysis, we argued that in the absence of external effects limiting contrail impact (e.g., large-scale subsidence or higher level clouds) the lifetime-integrated radiative impact of a contrail is governed, to lowest order, by the product of Nice and the depth of the supersaturated layer below flight level (Lewellen 2014). While other factors such as increased ambient relative humidity with respect to ice (RHi) or temperature may strongly affect optical depth (by producing larger crystals) that influence is canceled to leading order in the time-integrated effects by a decrease in lifetime (from faster sedimentation).
The initial formation of ice crystals depends on complex processing in the exhaust plume typically completed within the first second after emission (Kärcher et al. 1996). The basic pathway involves activation and droplet growth on aerosol when conditions in the plume exceed water saturation, followed by homogeneous freezing. Exhaust soot is believed to be the dominant aerosol in this process for current aircraft. Exhaust products are most conveniently normalized per mass of fuel burned to form an “emission index” EI. For soot number this is typically in the range EIsoot ~ 2 × 1014 to 2 × 1015 kg−1 fuel for modern jetliners (e.g., Anderson et al. 1998; Schröder et al. 1998; Moore et al. 2017). Ultrafine exhaust aerosol (UFA; e.g., few-nanometer-scale sulfuric acid droplets) and ambient aerosol entrained into the plume contribute as well. The former is typically much more plentiful than soot, EIUFA ~ 1016 to more than 1017 (e.g., Anderson et al. 1998; Schröder et al. 2000) or even greater if chemi-ions play a significant role (Yu and Turco 1997), but harder to activate because of its small size. The ambient aerosol is secondary to soot because it typically occurs in much smaller numbers in the early plume. With future changes in engine technology and fuel blends the relative importance of these contributions may change significantly. One of the motivations for the present work is to better assess how changes in EIsoot, EIUFA and ambient aerosol may impact the effective emission index for ice number EIiceno.
It is known that the mean plume mixing rate is important in determining EIiceno, but not whether the details of that mixing process are critical in any regimes. While LES has frequently been used for later stages of contrail evolution (see, e.g., Paoli and Shariff 2016), the ice initiation stage has been treated most extensively by coupling a complex microphysics model with a simple box-model plume dilution (Kärcher and Yu 2009; Wong and Miake-Lye 2010). In the present work we address this issue by using LES of the exhaust plume to simulate the formation and subsequent evolution of EIiceno. This extends our prior LES treatment, where values of EIiceno after the formation stage were simply imposed (Lewellen et al. 2014, and references therein), backward in time to include nucleation of ice crystals on exhaust or ambient aerosol. Simplified microphysics and fluid mechanics are employed so that simulations over a large parameter space could be performed. Paoli et al. (2004) performed LES of contrails in the jet regime but with more limited ice-only microphysics, implemented along sample Lagrangian parcels tracked in the flow, and for only a few cases. This study was extended in Paoli et al. (2013) to consider variations in EIsoot and engine configuration and by including a simple activation criterion from soot to ice when water saturation is first exceeded in a parcel. As in those works we employ a “temporal” rather than a “spatial” representation for the simulations: following a length of exhaust plume in time in the frame of reference of the ambient atmosphere rather than including the entire downstream development of the jet in the engine frame of reference. As such we do not attempt to reproduce the visual appearance of the near-field contrail (i.e., within a wingspan or so of the engine), but do try to realistically represent mixing rates affecting the microphysics development. Khou et al. (2017) perform contrail simulations in a true spatial representation, though for Reynolds averaged Navier–Stokes (RANS) rather than LES, restricted to the very near field, with different compromises on the microphysics, and again only for a few cases.
Our LES model and simulation procedures are discussed in section 2, with some details deferred to the appendix. Sections 3 and 4 present the principal simulation results and a more detailed theoretical analysis in some limits to help explain those results. Section 5 includes further discussion of connections to prior work and some conclusions.
2. LES model
Our starting point is the contrail LES model described in Lewellen et al. (2014). It has been used extensively for exploring aircraft wake and contrail dynamics from the period beginning several wingspans behind the aircraft after ice crystal formation through to contrail demise possibly many hours later (Lewellen and Lewellen 1996, 2001a,b; Lewellen et al. 2014; Lewellen 2014). Bin-resolved microphysics is used for the diffusional growth and sublimation of ice crystals, using subroutines adapted from the NASA Ames Community Aerosol and Radiation Model for Atmospheres (CARMA) (e.g., Jensen et al. 1994; Ackerman et al. 1995). The LES model uses several refinements for efficiency and accuracy including: a TKE based subgrid model including rotational damping effects; multiple run segments with different grids and domains; stretched grid spacing and active grid tracking; ice bins defined locally in time and space; and dynamic substepping of microphysics as needed. To extend the treatment to ice crystal formation in the exhaust jet, several additions and compromises were made. We consider those involving the fluid mechanics, initialization, and microphysics in turn.
a. Fluid mechanics
In the free exhaust jet, pressures rapidly revert to the environmental value pe within a few jet diameters, but density and temperature variations remain significant for much longer. To model these changes within an incompressible flow solver we implemented a variable-density low-Mach-number approximation (see, e.g., Bell and Marcus 1992). Essentially the local density is evolved prognostically as a conserved scalar, temperature is diagnosed from it via the ideal gas law, and the divergence free condition ∇ ⋅ u = 0 is retained on the velocity. The hot jet core is essentially sonic (i.e., Mach 1) at the engine exhaust plane; however, as shown by Papamoschou and Roshko (1988), compressible effects on mixing within a shear layer depend on the local “convective Mach number” Mc—that is, the Mach number in a coordinate system convecting with the velocity of the dominant structures of the shear layer. For the B737-sized engine flying at Mach 0.8 that we use as our base case below, Mc = 0.145 in the core-bypass1 shear layer, and Mc = 0.218 in the bypass-ambient shear layer initially, dropping with time as the exhaust jet mixes. At these levels of Mc the experimental data shown in Papamoschou and Roshko (1988, their Figs. 16 and 17) imply that the shear-layer mixing rates would be the same in the compressible case as they would be in the incompressible with the density variation accounted for. Since deviations to it scale as , retaining ∇ ⋅ u = 0 is not unreasonable for our purposes. A density variable based on the conserved total specific enthalpy is employed in the LES to approximately include effects of kinetic energy dissipation and latent heating, as described in the appendix.
A more serious compromise is our use of a “temporal” rather than a “spatial” representation of the jet for the simulations (following, e.g., Paoli et al. 2003). In the temporal representation a length of exhaust plume is followed in time in the frame of reference of the ambient atmosphere, with periodic boundary conditions in the axial direction. In contrast the spatial representation would be in the engine frame of reference, with downstream domain extended to encompass the plume ages of interest. The temporal approximation presumes a jet with a small expansion angle2 (as is true here); it requires ~2–3 orders of magnitude less computing resources in the present case than the spatial representation at comparable resolution. There are, however, significant differences that must be considered. Conservation laws differ. For a conserved tracer the integrated concentration per jet length is conserved in time in the temporal case but increases over time in the spatial case as ~Vflt/Va(t), the ratio of flight velocity to axial jet velocity: there is effectively an axial compression of the jet. For the spatial representation it is instead the mean flux of the tracer across a plane normal to the jet axis that is conserved as one moves downstream. Further, the relation between exhaust parcel age since emission and distance downstream from the exit is different in the two representations: the age can vary across the plume at a fixed downstream station in the spatial case but will be constant within the domain in the temporal. Mixing between parcels of different ages could potentially impact microphysics results. To quantitatively assess these differences we performed a spatial-representation LES without microphysics for comparison. Details and further discussion are given in the appendix.
Jet simulations were initialized using results from a 2D compressible axisymmetric RANS calculation of an exhaust plume from a single uninstalled engine of the approximate size and thrust needed for a B737-sized aircraft, provided to us by Boeing and using a generic engine cycle model, generic exhaust nozzle and commercial code (CFD++). Conditions are for ambient pressure p* ≡ 238.4 hPa, temperature T* ≡ 218.8 K, cruise velocity of Mach 0.8 (237 m s−1), fuel flow ff = 0.32 kg s−1, specific heat of combustion Q = 42.9 MJ kg−1 and emission index for water . Analyzing the fields for mass, momentum and energy conservation we deduced a thrust of 18.84 kN and propulsion efficiency of η = 0.325. The initial jet core and bypass diameters are ≈0.6 and 1.2 m, respectively. The RANS pressure field rapidly equilibrates to p* downstream as expected. Nonetheless, there are pressure (and accompanying temperature and velocity) fluctuations in space in the first 5–10 m downstream of the engine exit. Variable combinations that are invariant under adiabatic pressure changes [potential temperature , where cp is the isobaric specific heat and Ra the gas constant for air; “potential” density ρs ≡ pe/Raθ and “potential” velocity us ≡ (ρ/ρs)u] are free from these fluctuations. We then use the RANS θ, ρs, and us fields 1 m downstream from the engine exit plane3 to initialize the temporal LES T, ρ, and u fields to be consistent with the variable density approximation. Small velocity perturbations are added using a random streamfunction to preserve continuity. The RANS included a conserved exhaust tracer, which we scale appropriately to initialize LES aerosol and moisture fields. Interpolations from the axisymmetric RANS to the 3D LES are done so as to exactly preserve core and bypass integrals of mass, axial momentum, and exhaust tracer. An alternative initialization based on the CFM56–3 jet-exit-plane conditions given in Garnier et al. (1997) (adjusted to ambient , , and Mach number from the RANS) was also implemented. Multiple LES test cases with each showed quite consistent results between the two. Other LES tests showed no great sensitivity in final EIiceno to details of the initial aerosol spatial distribution (e.g., restricting them to the inner or outer jet core), fraction of heat transferred from the core to the bypass flow within the engine, or precise thrust level. Procedures to initialize with environmental pe and Te differing from the base values and or with an engine size scaled by a factor sE are described in the appendix.
For longer simulation times the wake vortices influence the mixing4 and hence potentially the microphysics. We include these as discussed in Lewellen et al. (2014), using velocity fields from 2D RANS wake runs of a two-engine jet scaled to a B737-300 wingspan (28.88 m), provided by Boeing. In these fields the axial exhaust jet velocity and turbulence are zeroed (using the exhaust tracer for identification) to be replaced by the jet-only LES results. The single jet temporal LES are performed out to 0.53 s (to best correspond to development five wingspans downstream; see appendix). These results are then axially compressed as appropriate for a far-field temporal representation and periodically continued onto a longer domain (see appendix), properly positioned within the wake field, and reflected for the second engine. Ambient stratification is added at this stage as well. Throughout this procedure the local thermodynamic state, and aerosol numbers and properties, are conserved. The simulations can then be continued as in Lewellen et al. (2014) (with most performed in the Q3D mode discussed there, defined by more limited downstream domains), with the choice of microphysics depending on whether the formation of ice crystals has concluded yet or not. Initial sensitivity tests were performed to determine modest but still adequate grid resolutions and domain sizes for all run segments; defaults are given in Table 1 and discussed further in the appendix. Note that the length scales here are much smaller than those usually associated with a cloud-resolving LES. Sample runs with potentially the most demanding resolution requirements were later repeated with refined (×2) resolution and showed quite consistent results (e.g., those included in Fig. 6 below).
The complex microphysics within the exhaust jet must necessarily be simplified for use within the LES. We assume (and vary for sensitivity studies) effective emission indices for aerosol content as present in the exhaust plume shortly before saturation conditions are encountered. Initial formation and chemical processing of the aerosol within the engine or exhaust plume at higher temperatures is not modeled. Particle scavenging/coagulation is omitted as well, assumed significant for our constituents only for initial aerosol formation in the hot jet (not modeled) or on longer time scales not impacting the ice-number formation itself (see, e.g., Fig. 3 in Kärcher 1998). Ice crystal formation is modeled as a two-step process: droplet growth on the postulated aerosol, followed by homogeneous freezing. A size-resolved (“binned”) representation is used for both the droplets and ice crystals (each assumed spherical), so size spectra are determined by the underlying equations without any spectral-shape assumptions. A mass ratio of 2 between neighboring bins is used for the simulations included here. Given the temperature range involved, local T dependence was added as needed for some parameters normally treated as constants within CARMA.
Aerosol activation and droplet growth is modeled via the κ–Köhler parameterization of Petters and Kreidenweis (2007) in which the equilibrium saturation ratio with respect to water over a solution droplet is
for droplet radius r, water activity aw, Kelvin parameter for solution droplets , solution–air interface surface tension σs/a, water vapor gas constant Rυ, and liquid water density ρw. Because of competing solute and curvature effects Swd will have a maximum at some critical radius r*. The critical saturation ratio for aerosol activation, , is the value of Swd at r = r*. The aerosol properties are defined in this approximation solely by the dry aerosol core radius rd and hygroscopicity parameter κ.
For the present study we mostly consider three idealized aerosol types: a coated aircraft “soot” with κ = 0.005 [following the choice of Kärcher et al. (2015)] and rd = 20 nm; an ambient CCN with rd = 20 nm and κ = 0.9 [consistent with the value for H2SO4 given in Petters and Kreidenweis (2007)]; and an ultrafine aerosol (UFA) with rd = 2 nm and κ = 0.9. For reference the at 218.8 K for these three choices is respectively 1.065, 1.0099, and 1.347. For our base conditions quite large Sw are potentially reached, as shown in Fig. 1a; thus even quite small aerosol with well above 1 can potentially activate, grow, and freeze. The for the nominal “CCN” is sufficiently close to one that further increase in rd from the chosen 20 nm has little practical effect (i.e., it can stand in for larger ambient aerosol as well, justifying the CCN label). Aerosol numbers are specified via emission indices for the exhaust aerosol and environmental number densities for the ambient. The numerical cost of binned microphysics within the LES is high because of the number of 3D fields it requires. This precludes, for example, using a continuous size spectrum of dry soot in a simulation (requiring then an entire family of droplet distributions). Nonetheless, we can still meaningfully test the sensitivity of EIiceno to the shape of the core-size spectrum with simpler examples (see the discussion of Fig. 2 in section 3a below) and find that it is not large. Likewise, while not matching the complexity of a broad set of aerosol that may be encountered, the use of only two values of κ (one on the low side for relevant aerosol, one on the high side) still permits an assessment of the sensitivity of EIiceno to κ. It is important to note that limiting to a unimodal dry core size does not imply a unimodal droplet or ice size spectrum even when restricted to single grid cells of the simulation: the binned representations for both droplets and ice permit continuous spectra, which the mixing between different parcels with different histories then naturally produces.
For homogeneous freezing we assume the probability for a droplet containing a volume of water Vw to freeze in time interval dt is 1 − exp(−JhomVwdt). For Jhom we use the water activity dependent parameterization of Koop et al. (2000) with aw given by (2). Riechers et al. (2013) argue that values for Jhom in prior studies are significantly too high. The temperature dependence they find for pure water (the only case they consider) agrees with that in Koop et al. (2000) but at a magnitude two orders lower. Accordingly, we added a scaling factor fJH to Jhom in the LES and use fJH = 0.01 as a default value; extensive sets of runs with fJH = 1 were also performed to assess sensitivities to homogeneous freezing rate in different regimes. Generally we have found that sensitivity to be low: Jhom increases extremely rapidly with decreasing T so even a large change in fJH corresponds to only a small shift in time. Typically, essentially all of the activated aerosol freezes so that the activation step mostly governs the number of ice crystals formed. Within the LES, parcels can mix with both ambient and core air allowing both cooling and rewarming. Accordingly, we allow for ice to melt to solution droplets that can later activate and freeze again. Test simulations showed that this “aerosol recycling” can be important at times to the final EIiceno.
d. Box-model approximation
It is useful to have a box-model representation that may be directly compared to the LES and used for more extensive parametric studies when it performs well. This requires a precise definition of the box-model volume and averaging. Exhaust products (aerosol and water vapor) play a critical role for contrail formation so we choose an exhaust-weighted average to define box-model quantities. A conserved exhaust tracer concentration field CE is included in the LES. We then define the exhaust-weighted average of a local field ϕ as
We interpret box quantities as approximations to the exhaust-weighted averages, ϕbox ≈ ⟨ϕ⟩E. From this association the proper definition for box plume volume follows by demanding , whence
In the box model the fluid mechanics contributes only through this dilution history, which can be obtained from a separate LES. Then for a field ϕ that is conserved by the fluid mechanics and has environmental value ϕe, the contribution to the evolution of the box-model variable ϕbox from dilution alone follows from its conservation to be
Its evolution due to microphysical processes is simulated exactly as in the full LES. The use of the exhaust-weighted average represents a far more robust way of determining the plume volume than trying to define plume boundaries (e.g., with some threshold) and integrating. Applied to the limit of a uniformly mixed plume (i.e., a top-hat distribution for CE), ⟨ϕ⟩E and as defined give the traditional box-model values. When comparing box model and LES results for a specific turbulence realization we use a dilution history obtained from the latter for use in the former. For more general purposes we average over multiple LES dilution histories (see Fig. 1b for an example). Note that it is not necessary to run new LES to obtain dilution histories for changes in engine size (sE) or Te; an existing history can be scaled consistently to account for these changes.5
3. Principal simulation results
An extensive set of simulations of contrail formation were performed varying ambient conditions (Te, RHi, pe, stratification dθ/dz), aerosol (rd, κ, EIsoot, EIUFA, or combination EIaero ≡ EIsoot + EIUFA, environmental number density Ne), engine size scaling (sE), and numerical treatment (LES vs box model, resolution, turbulent realization, freezing rate). Ranges for ambient conditions and aerosol numbers were chosen to include typical contrail conditions (see, e.g., Kärcher 2018) as well as potentially informative extremes. A sampling of results is provided here to illustrate some of the main factors involved in determining EIiceno. Where unspecified, base conditions were used of sE = 1, pe = 238.4 hPa, Te = 218.8 K, RHi = 110%, dθ/dz = 2.5 K km−1, and fJH = 0.01.
a. Exhaust aerosol
Figure 2 shows sample EIiceno evolution beginning from soot aerosol, varying the dry-size specifications, turbulence realization, and model treatment. A relatively large EIsoot was chosen for added sensitivity (for reasons that will become clear). Three features are apparent that we have found to be generally true in our simulations. First, with the modest domain sizes used and importance of turbulent mixing, there can be significant variation of EIiceno with turbulence realization. This provides a basic uncertainty level that impacts any potential conclusions to be drawn. Second, the structure of the aerosol dry-core-size spectrum (whether monodisperse or bimodal in the case illustrated) is of at most secondary importance, less than the spread from turbulence realizations. Third, when the ice crystals originate primarily from exhaust aerosol the box model EIiceno can deviate significantly from the corresponding LES result. Different parcels in the plume produce ice crystals at different times. As these parcels interact via mixing, the moisture consumption from early crystal growth in some parcels can prevent the activation and freezing of soot aerosol in other parcels that later mix in, leading to only partial conversion of EIsoot to EIiceno. This dynamic is absent in the box-model representation where all of the exhaust aerosol experiences the same supersaturation conditions at the same time. In that case different aerosol sizes or types can still lead to competition for available moisture favoring the conversion of larger or more hygroscopic aerosol to ice crystals, but that competition is generally less effective in preventing conversion than is the presence of crystals that had a head start on their growth. Thus, the much narrower time window for ice number production seen in the box model leads to greater EIiceno than in the LES. And similarly, for the LES results in Fig. 2 the turbulence realizations corresponding to longer mixing scales provide a greater suppression of EIiceno production than the competition between different sized cores in the bimodal dry-core-size distribution provides. Figure 3 shows a related example. The eventual EIiceno from the LES is approximately the same whether the exhaust aerosol is all soot, all UFA or a 50/50 mixture. Further, in the mixed case although only ~35% of the aerosol is converted to ice, a significant fraction of the ice still originates on the UFA even though the larger soot aerosol is more easily activated. This is contrary to what is found in the box model (not shown) where if only a modest fraction of aerosol is converted to ice it will all be preferentially on soot rather than UFA. Again, the primary competition limiting the conversion of aerosol to ice in the LES involves the mixing of parcels with different histories (early crystal formation versus late), more so than between different aerosol types within the same parcel.
Figure 4 shows summary results from a set of LES of contrail formation on exhaust aerosol. Given the findings noted above these have all, for efficiency and simplicity, been initialized with only a single dry aerosol size and type for each run. Values of EIaero much larger than are typically expected have been included to clarify behavior in that limit. Several results are apparent beyond those already noted. First, for Te = 225 K EIiceno is significantly reduced and is highly sensitive to effectively all variations considered, including turbulence realization. This is near the contrail formation threshold (i.e., the ambient temperature at which the exhaust plume will just become supersaturated with respect to water at some point during its dilution, which for the cases here is 225.2 or 226.3 K for RHi = 110% or 130%, respectively). Second, away from threshold conditions, all EIaero is converted to EIiceno if EIaero is small enough, irrespective of the details of aerosol type or ambient conditions. As EIaero is increased, and with it the competition for available moisture in the plume, the fraction of aerosol converted to ice reduces. The degree of reduction is less for lower Te (at fixed RHi), consistent with the greater supersaturation levels in the plume that are present (cf. Fig. 1a). The effect of ambient RHi is slight, except near threshold conditions, again consistent with its impact on supersaturation in the plume (Fig. 1a). As with the specific cases illustrated in Figs. 2 and 3 we see, for the broader simulation set, the impact of turbulence realization, lesser sensitivity to aerosol type, and the overprediction of EIiceno from the box model relative to the corresponding LES.
In Fig. 5a as EIaero increases, EIiceno appears to saturate toward a Te dependent limit. Moreover (Fig. 5b), ice crystal losses due to the sublimation of the smallest crystals through the wake vortex regime of contrail evolution further compress the range of EIiceno that is encountered [a result shown already in Lewellen et al. (2014), where these are referred to as “in situ” crystal losses to distinguish them from losses due to sedimentation, which do not occur on the short time scales considered here]. This is because these losses increase with increasing ice crystal number density (Lewellen 2012). While these losses are often considered to be a consequence of adiabatic heating from the descent of the aircraft wake forcing contrail conditions below ice supersaturation (Sussmann and Gierens 1999; Lewellen and Lewellen 2001a), this condition is not required. For large enough ice number density, in situ losses will occur even given a slow growth of ice mass within a parcel, driven by Kelvin-effect-dependent asymmetry in competition for moisture between larger and smaller ice crystals (Lewellen 2012). Thus, for large enough EIiceno these losses ensue well before the wake descent becomes significant,6 as shown in Fig. 6. This in situ loss dynamic implies an upper limit on EIiceno at the end of the jet-dominated regime (~10 s) regardless of how large EIiceno is after initialization, or what mechanisms produce it; we estimate this “in situ loss limit” in section 4 below. The crystal losses are augmented by the adiabatic heating from wake descent as it becomes significant in lowering growth in ice mass. This is responsible, for example, for the lapse rate dependence seen in Fig. 6 after ~50 s: an increased lapse rate reduces the maximum depth of wake descent and therefore reduces the potential for crystal loss. The fraction of EIiceno lost through the wake vortex regime depends on multiple physical factors. In our simulation experience this includes strong sensitivity in at least some regimes to Te, RHi, aircraft, stratification, and ice crystal numbers. Empirical parameterizations based on limited sampling of parameter space [see Unterstrasser (2016) for the most extensive example] should be interpreted accordingly.
b. Ambient aerosol
Figure 7 shows sample EIiceno evolution for cases with only ambient aerosol. As in the cases with exhaust aerosol (and for the same reasons), colder Te leads to greater EIiceno. Two contrasts relative to the exhaust aerosol cases are worth noting. First is the dramatically extended period over which the ice crystal initialization sometimes takes place.7 This is a consequence of the continual addition of aerosol to the exhaust plume and of the extended period where potentially Sw > 1 in the plume at colder temperatures (cf. Fig. 1). Second is the greater degree to which the box-model results mirror the LES, a result we find quite generally when EIiceno arises dominantly from ambient aerosol. Because of the gradual incorporation of aerosol into the exhaust plume in this case the box model here does include a fair representation of the competition for moisture between ice initiated early in the plume evolution and aerosol mixed in only later.
The dependence of EIiceno on Ne is more clearly illustrated by plotting the ratio EIiceno/Ne as in Fig. 8. Values of Ne above and below typical expectations in the upper troposphere are included to clarify the behavior in these limits. EIiceno(t) is bounded above by the total ambient aerosol number contained within the plume segment at time t (for unit fuel consumption). This gives a bound on EIiceno/Ne that is just a measure of evolving plume volume, set purely by plume dilution. With increasing Ne, EIiceno/Ne breaks off from this limit at progressively earlier times as moisture within the plume is more rapidly consumed by ice growth. For small enough Ne, EIiceno/Ne breaks from this limit purely due to dilution reducing Sw within the plume (cf. Fig. 1) rather than ice growth; below this point the ratio becomes independent of Ne.
Figure 9 summarizes results from a set of LES with ambient aerosol only. It shows again the relative success of the box model in these cases. There is negligible crystal loss seen through the wake-vortex regime (ending by 300 s) in these cases, a consequence of the lower ice crystal number densities present (relative to most cases in Fig. 5). This is also partly due to the B737-sized aircraft being simulated: a larger aircraft leading to greater wake-vortex descent can give rise to crystal loss starting from smaller initial EIiceno. Given the success of the box model it is possible to more thoroughly explore EIiceno production from ambient aerosol over a large parameter space. Figure 10 provides examples varying Ne, Te, altitude and aerosol size. Some of the dependencies seen, particularly in the low and high Ne limits, can be understood from simpler models considered in section 4 below.
c. Exhaust plus ambient aerosol
If EIaero is sufficiently large that only a fraction of the activatable aerosol converts to ice, this conversion will generally complete before plume expansion incorporates enough ambient aerosol to significantly impact EIiceno. On the other hand, if EIaero is well below this threshold then even fine ambient aerosol could contribute significantly (Fig. 11). In this case we find that the box model reasonably represents the LES results: correctly predicting the EIiceno contribution from the exhaust aerosol for simple reasons (for low-enough EIaero, all is converted) and that from the ambient aerosol for the reasons discussed above. The box model then allows a more complete exploration of this parameter space, as shown, e.g., in Fig. 12. With both Ne and modest EIaero present the resulting EIiceno will lie below the sum of the EIiceno produced by each if present in isolation, but generally only moderately so. There can be exceptions, however, such as for Te = 205 K and Ne ~ 108 m−3 in Fig. 12. There the addition of EIaero of 1013 kg−1 leads to less EIiceno than the ambient aerosol would have produced alone. At this low Te and modest Ne the EIiceno production is spread over ~10 s (cf. Fig. 7); over this long a time the early ice production from the addition of EIaero can effectively prevent much of the activation of ambient aerosol.
d. Scaling with engine size
Figure 13 summarizes results from a series of LES scaling engine size by a factor sE as discussed in the appendix. For low-enough EIaero all the aerosol converts to ice and EIiceno is independent of sE. For larger EIaero, the EIiceno after formation falls with increasing engine size (sE). The larger diameter jet from a larger engine dilutes on a longer time scale.8 This provides more time for ice crystals produced early in the process to consume moisture, thus preventing later aerosol from activating, depressing EIiceno. Figure 14a shows related results from box-model simulations of EIiceno production from ambient aerosol (where the box-model results have proved trustworthy). Again for larger Ne a larger engine produces a smaller EIiceno than a smaller (but otherwise comparable) engine. Further there is a clear similarity in the effects on the falloff of EIiceno/Ne in Fig. 14a due to a decrease in engine size versus due to a decrease in Ne (curves for different sE have similar shapes, differing by a shift in Ne). This fall off depends on the relative magnitudes of the time interval during which a parcel potentially activates aerosol (call it τmix since it depends on plume mixing rates), and the time interval required for a population of ice crystals to grow enough to significantly deplete the available moisture [for which we can use the phase relaxation time τfc introduced by Khvorostyanov and Sassen (1998)]. Restricting to variations with sE and Ne while keeping other factors (T, RHi, κ, etc.) constant, τmix scales as sE and τfc approximately as .9 If the dependence of EIiceno/Ne on τmix and τfc is only through their ratio (scaling as ) then plotting the results against any function of this ratio [e.g., ] should collapse the sE dependence in this regime. Figure 14b confirms this to reasonable approximation.
4. Regimes and idealized limits
The simulation results may be grouped into three general regimes: near threshold (NT), low aerosol number density (LND) where essentially all the activatable aerosol (from either exhaust or ambient) contributes to EIiceno, and high aerosol number density (HND) where competition between the aerosol limits ice conversion to only a fraction. We develop a quantitative formulation in the LND regime first, followed by an upper limit in the HND regime; the NT regime we have found to be dependent on too many parameters to be amenable to simplified analysis.
We define the LND limit to be when ice numbers remain small enough and mixing rates fast enough that ice growth never becomes a significant moisture sink during the ice nucleation window. We then approximate the initial EIiceno production as follows. We assume that aerosol nucleates ice if and only if it activates (i.e., encounters ). Consider an infinitesimal plume parcel following the products from burning a mass of fuel fp, diluting in an environment with total water mixing ratio qte (assuming no condensate there). The evolving mass of water in the parcel (Mwp) is related to the growing parcel air mass (Map) as
Combustion provides an increase in enthalpy in the parcel that is conserved thereafter, δHp ≡ Q(1 − η)fp for specific heat of combustion Q and propulsion efficiency η, so that the parcel’s temperature evolves over time as
The evolving saturation ratio with respect to water within the parcel is then
with ϵ ≡ Ra/Rυ and es(T) the saturation vapor pressure over liquid water at temperature T. Setting in (8) for our aerosol to ice transition assumption, using (7) and recombining gives an equation for the parcel temperature at the aerosol activation point ,
As in the closely related Schmidt–Appelman criterion for contrail formation (Schumann 1996) this depends critically on the slope G of the mixing line for the parcel in a plot of partial pressure of water vapor versus temperature.
In general for contrail conditions there are two solutions of (9) for , corresponding to when the parcel first reaches and when it last does so, which we denote Tp1 and Tp2, respectively. In the LND limit under our assumptions the number of ice crystals produced from ambient aerosol in the parcel will be given by NeMap/ρe with Map given by (7) with Tp = Tp2. Different parcels reach this point at different times with different mixing histories. Nonetheless, given the LND limit assumption that the growth of condensate during this time period is not a significant moisture sink, the total ice number production is a simple superposition over all parcels, independent of the details of the plume mixing process. Summing parcels for unit fuel burn and adding in exhaust aerosol (but small enough to remain in the LND limit) this gives the prediction
The same equation with Tp2 replaced by Tp1 is in principle relevant in the opposite limit, i.e., Ne large enough and mixing time scales long enough that moisture depletion sufficient to cut off activation occurs before Tp dilutes much below Tp1. This is closer in spirit to the approximation for ice crystal number from ambient aerosol considered in Kärcher et al. (2015). In the present simulations we approach but do not actually reach this limit, even for the largest Ne considered; nonetheless, as discussed below, it explains some behavior encountered for large Ne.
Both solution branches are plotted in Fig. 15 for sample aerosol and ambient conditions. To more explicitly illustrate the parameter dependencies we provide an approximate closed form solution for (11) and (9) arrived at by expanding the parameterization for es(T) given by Buck (1981) about T = Te and solving; after some algebra,
where Swe is the environmental water saturation ratio. In the term in curly brackets Te should be in kelvins and the bracketed term is dimensionless. The approximation (12) agrees with (11) at colder temperatures but loses accuracy nearer to threshold conditions (e.g., for the conditions in Fig. 15 it is good to 1% or better below 216 K but off by as much as 10% for some cases at 223 K). The dependence on rd and κ is solely through . The first of the two terms in square brackets dominates the second, particularly at lower temperatures. If the second term is dropped as a further approximation the result becomes independent of Q, η, and pe.
For cases with exhaust aerosol alone the Eq. (11) LND prediction is simply EIiceno = EIaero, already seen in the LES results (e.g., the green lines in Figs. 5 and 13). More nontrivially, much of the results for EIiceno produced from ambient aerosol (e.g., Figs. 10 and 14) can be understood by comparing with Fig. 15 and Eq. (12). The plateaus seen in EIiceno/Ne for low Ne and their dependence on Te and rd are in good quantitative agreement with the corresponding LND limit model predictions, particularly for the smaller aerosol (see Fig. 10a). To a lesser extent the “Tp1” branch of the solution in Fig. 15 illuminates the behavior of the simulation results for large Ne. In particular the switchover in Fig. 10a from EIiceno/Ne increasing with decreasing Te at low Ne to increasing with increasing Te at high Ne is predicted by the corresponding switch in behaviors seen in the “Tp2” and “Tp1” branches in Fig. 15. The LND limit model is not in exact agreement with the simulations even for low Ne. At colder temperatures for larger aerosol the LND limit model tends to slightly underpredict EIiceno/Ne from the simulations (up to 15% below has been seen), and the simulations show a small dependence of EIiceno/Ne on plume mixing time scale even for low Ne that is not predicted by the LND model (cf. Fig. 14). A more detailed examination shows that these arise from the LND model assumption that is necessary and sufficient for ice crystal production. That it is not necessary leads to the first discrepancy (in fact, as for natural cirrus production, for large enough RHi some hygroscopic aerosol can grow enough to freeze even without water supersaturation). That it is not sufficient leads to the second (for shorter mixing time scales some activated aerosol may not have sufficient time to grow and freeze before mixing lowers Sw sufficiently to drive enough evaporation to make the droplet too small to freeze).
For large enough aerosol number, moisture depletion from early ice growth limits EIiceno to only a fraction of that implied by the activatable aerosol present. We will not try to formulate an analytic estimate for EIiceno in the HND regime here [see, e.g., Kärcher et al. (2015) for recent attempts to do so]. Instead we consider the upper limit at the end of the jet regime from quite different physics, the in situ crystal loss limit mentioned in the discussion of Fig. 6 earlier. We employ the exact solution to the diffusional growth equations derived in Lewellen (2012) for the evolution of a population of ice crystals in the “slow growth” regime.10 We assume a uniform plume parcel in which, by time t0, all ice crystals have been produced and slow-growth evolution has commenced. Then from Eqs. (28) and (8) of Lewellen (2012) (and normalizing to form an emission index),
with the analog for ice of in (2), qsi the ice saturation mixing ratio, Dυ the water vapor diffusion constant, ME(t) the equilibrium ice mass for unit mass of fuel burn, and Λ a parameter depending on the growth rate of ME(t). For t ≫ t0 the plume T will have relaxed to good approximation to Te and ME(t) will be given by
with A(t) the plume cross-sectional area and fm the fuel mass per length of exhaust plume. To get an upper limit we assume arbitrarily large initial EIiceno(t0), so that is negligibly small leaving
For a plume mixing into a supersaturated environment in the slow-growth regime with no heating of the parcel, Λ will be in the range ~1/6 to 1/2. Several complications discussed in Lewellen (2012) could also be swept into this parameter (e.g., kinetic theory corrections for small crystals, uncertainties in the ice deposition coefficient, bin-resolution effects when comparing with simulations), many effectively lowering the value of Λ. For an order of magnitude estimate we simply choose a nominal value of Λ ≈ 0.2 here. Using LES results for A(t)/fm in (15) this gives the limit lines shown in Fig. 6. These lie above the LES results as expected but at ~10 s only by about a factor of 4 and capturing the same Te dependence. The results from (15) have not been carried to later times in the figure because neither the effects of parcel descent or of growing plume inhomogeneity (both increasing in importance in the LES over time) are included. The former leads to a reduction in ME(t), which is largely responsible for the increased falloff rate seen in the LES results as the wake vortices descend most rapidly.
5. Discussion and conclusions
We have, to our knowledge for the first time, performed a large set of LES of contrail initialization, focusing on the impact of mixing on EIiceno for different aerosol choices and ambient conditions. The results support many of the findings of the box-model study of Kärcher and Yu (2009) (e.g., the near-complete nucleation of soot to ice for EIsoot ~ 1015 or less away from near-threshold conditions; the only partial conversion of soot to ice for much larger EIsoot with the fraction decreasing with increasing temperature; the rising (and significant) contribution of ambient aerosol to EIiceno with decreasing temperature and EIsoot), but not all (e.g., their conclusion that UFA do not contribute to EIiceno in the soot-rich regime of EIsoot ~ 1015 or above). We have shown that well below the contrail threshold temperature and for sufficiently low numbers of activatable exhaust or ambient aerosol, EIiceno is essentially independent of any details of plume mixing and can be computed analytically to a good approximation. That analytic expression [(11) and (12)] makes explicit the dependence of EIiceno on the ambient atmospheric state, aircraft, aerosol, and fuel properties in the low-aerosol-number limit. Given larger aerosol numbers we find that faster plume mixing rates (e.g., arising from smaller engine size) lead to larger EIiceno. Even for larger aerosol numbers we have found that LES results for EIiceno are approximated well by a box-model representation if EIiceno arises dominantly from ambient aerosol. On the other hand, for large enough numbers of activatable exhaust aerosol we find that a box-model representation misses the critical competition between ice crystals formed early and exhaust aerosol mixed with them later. This causes both an overestimation of EIiceno relative to the LES results and an underrepresentation of the fraction of EIiceno that may arise from ultrafine aerosol (UFA) even in a soot rich regime.11 One consequence of the latter (and a potentially testable result) is that given copious UFA a reduction in EIsoot alone may not significantly change EIiceno: additional early ice-number production from UFA simply replaces the lost early production from soot in nearly one-to-one fashion (see, e.g., Fig. 3).
Our results for EIiceno are at least broadly consistent with estimates from in situ measurements. For example, the compendium of measurements collected in Schumann et al. (2017) taken over a large range of conditions show an upper limit on EIiceno determinations of a few times 1015 kg−1 (their Fig. 8), consistent with expectations here for aged contrails in Fig. 5b. More recent measurements of Kleine et al. (2018) (for Te = 215.5 K and a comparable aircraft to what we consider here) suggest an EIsoot ~ 5 × 1015 and EIiceno of ~3 × 1015 initially after formation, dropping to ~1 × 1015 with age as the wake descends, consistent with expected ranges here in Fig. 5 at Te = 215.5 K and EIaero = 5 × 1015.12 Prospects for more telling quantitative comparisons are severely limited by the sheer number of important physical parameters involved that would need to be accurately measured (particularly for conditions approaching contrail threshold) and by the significant inhomogeneities in the plume being sampled.13
Our simulations suggest that it will be difficult to determine the type of aerosol responsible for the ice nucleation solely from an examination of a contrail of a few minutes age or older. Even scenarios differing dramatically in initialization such as in Fig. 16 [showing a range of ~1000 in EIiceno(t) at 0.1 s] may show far less difference later (less than a factor of 4 spread at 300 s). In fact in our simulations, we do not find dramatically larger values for EIiceno after the wake-vortex regime regardless of the initial choice for aerosol number or their type, due largely to the importance of in situ crystal losses (cf. Fig. 5). This suggests, for example, that measurements in the wake-vortex regime showing comparable EIsoot and EIiceno are not definitive evidence for the box-model prediction in the soot-rich regime that ice numbers are actually controlled by soot numbers. In the simulations for EIsoot in the range measured for current jetliners (~2 × 1014 to 2 × 1015 kg−1 fuel) EIiceno also ends up in a comparable range (see Fig. 5b), even with initial ice nucleation sometimes on a different mix of aerosol or in higher numbers. Further, spatial and spectral analysis of our simulations has so far failed to identify any robust (even if complex) signature of the ice-number-formation stage that remains at the end of the wake-vortex regime. This suggests that remote optical measures in the near field behind the aircraft (where, for example, differences in Fig. 16 are still pronounced) may provide clearer evidence for the source of ice nuclei.
The simulation results suggest that a dramatic reduction in both EIsoot and EIUFA has the potential to reduce EIiceno from typical current levels (~1014–1015 kg−1) by a factor ~10 for persistent contrail conditions assuming ambient aerosol number densities of hundreds to thousands cm−3. The potential reduction is not greater in part because of the extended time over which ice originating from ambient aerosol can be nucleated, and could in fact be less because a much broader class of aerosol may be involved. Ultrafine few-nanometer-scale ambient aerosol are generally ignored in cirrus cloud studies because they cannot be activated by natural processes, but given the large saturation ratios that can occur in aircraft exhaust plumes they may become significant for contrail initialization. More extensive measurements of ultrafine aerosol concentrations in the upper troposphere (UT) and how they might correlate with ice supersaturated regions seems to us a critical need in assessing the practical lower bound on EIiceno that may be obtainable. Processes such as ion-induced nucleation from cosmic rays may enhance UFA concentrations in the UT (see, e.g., Kulmala et al. 2014). One of the few measurements of UFA in the UT that we are aware of, Mirme et al. (2010), shows median total ultrafine aerosol concentrations at 10-km altitude of ~2000 cm−3 (their Fig. 8) and maxima in the UT of ~3 × 104 cm−3 (cf. their Fig. 5). This motivated extending our sensitivity studies to include quite large values for Ne. The resulting EIiceno under different conditions can approach those within the soot-rich regime.
In prior work we highlighted the potential increased significance of “cold” subvisible contrails integrated over their lifetime (Lewellen et al. 2014; Lewellen 2014). This was based on the potentially much longer lifetimes and significant reduction in ice crystal losses during the wake vortex regime for these contrails. The present simulations further support this potential: for fixed RHi and within the contrail limits, colder temperatures foster nucleation of more ice crystals whether they originate from ambient or exhaust aerosol. At the opposite limit, for temperatures approaching the contrail threshold for given conditions, our simulations show large sensitivity of EIiceno to many more factors, including details of mixing events. This makes reliable predictions or comparisons in this regime problematic: small errors in modeling or measurement could easily lead to spurious conclusions. Fortunately, from the perspective of potential climate impact, near-threshold contrails are unlikely to represent a significant fraction of persistent contrails, both because of the modest fraction of occurrence of conditions within 1 to 2 K of the threshold in heavy flight areas (see, e.g., Bier and Burkhardt 2019) and the expectation that the fewer and larger crystals produced will lead to shorter contrail lifetimes.
Finally, some caveats on the present work. We did not try to simulate the initial aerosol production (e.g., soot number), which clearly results from complex and sensitive processes within the engine and early jet leading to parameter dependencies we have not considered. The observations of Moore et al. (2017), for example, show significant increases in aerosol emission indices with an increase in thrust and accompanying increase in fuel flow rate.14 In addition to the uncertainty in ambient number densities of UFA, the treatment of microphysical processes involving nanometer-scale aerosol based on bulk parameters, though a standard approach (see, e.g., Zhang et al. 2012), is more uncertain than the corresponding situation for larger, better-studied CCN. This includes treatment of the Kelvin effect and the κ–Köhler parameterization, and values of ice deposition coefficient and surface tension—all affecting some of our simulation results. Finally, while we believe that our findings on when details of turbulent mixing affect EIiceno are robust, the actual values found for EIiceno may be affected by our use of a temporal representation, low-Mach-number approximation, and modest domain for the simulations. In particular, in some regimes the uncertainty range in EIiceno from turbulence realizations and differences between spatial and temporal simulations can be significant.
This work summarizes research supported by the Boeing Company over several years. We thank Boeing personnel, particularly Steve Baughcum, Mikhail Danilin, and Doug DuBois, for critical discussions, guidance, data, and feedback throughout the course of this work. We also thank Roberto Paoli and two anonymous reviewers for comments on the manuscript leading to significant improvements in the presentation.
a. Variable density implementation
In the expansion of the jet at constant pressure the specific enthalpy ht ≡ cpT + u2/2 − Lsqi is conserved (where cp is the isobaric specific heat, Ls is the latent heat of sublimation of ice and qi is the ice mixing ratio).A1 We define our density variable in terms of ht, (with Ra the gas constant for air), so that we can include the effects of kinetic energy dissipation and latent heating when solving for T,
We have found in simulation tests that including these terms can affect EIiceno when it is mainly produced early in the jet evolution. Using in place of ρ is accurate at the level of the low Mach number approximation and gives the correct density asymptotically as the jet dilutes. This is in the same spirit as the Boussinesq approximation—using a more accurate treatment of density change where it is needed (the microphysics is very sensitive to T) and a simpler one where it is not (the fluid dynamics).
In the LES the velocity is updated in two steps:
where u* includes everything other than the pressure update and the pressure perturbation p′ is determined to satisfy the condition ∇ ⋅ u = 0. By reexpressing the pressure gradient term in (A2) we can satisfy ∇ ⋅ u = 0 by solving a Poisson equation even with the variable density:
In the numerical implementation the last term in (A3) is evaluated at the prior time step, consistent with the leapfrog time scheme for the velocities. We have not attempted to include variable density effects in the subgrid model, leaving it unaltered other than ensuring that the diffusion terms properly deal with gradients of momentum flux (not velocity flux) to insure momentum conservation.
b. Spatial versus temporal LES
A spatial LES of an exhaust jet was performed by adapting an open-domain LES designed for tornado simulations (see, e.g., Lewellen and Lewellen 2007). The variable density approximation and suitable inlet conditions were implemented to otherwise match the temporal simulations considered here. Stretched grid spacings were used with finest spacing of 6 cm. A jet length of 120 m was simulated, corresponding to a flow-through time of ~0.5 s. Large lateral domain dimensions of 160 m × 160 m were used so that far-field secondary flows arising to feed the jet entrainment would not unphysically bias the near-jet flow. In addition to the exhaust tracer CE an “aging tracer” Ca was implemented specifically to be used in determining exhaust-averaged parcel age since emission. It is initialized to zero at the inlet and treated by the fluid mechanics as a conserved tracer but with CE as a local source term, ∂Ca/∂t = [∂Ca/∂t]FM + CE. In a temporal simulation starting with Ca = 0, the local ratio Ca/CE is spatially constant and equal to the elapsed time t (by construction). In the spatial simulation this ratio locally gives the mean exhaust-weighted age since emission of that parcel. For a given downstream station y the exhaust-weighted, cross-stream and ensemble averaged age of the exhaust jet is then given by .
The simulation was run without Ca for 2 s to reach a quasi-steady state, then for 0.5 s for Ca to flush through the domain, and finally for 1.5 s to form time averages. Basic results are shown in Fig. A1. As noted earlier there is an increase in the integrated amount of a conserved exhaust species per jet length over time (seen in Fig. A1a) due to an effective axial compression of the jet. For comparison corresponding temporal simulations were run in two ways: conventionally starting from near-field jet conditions (NF); and scaled so as to match “far field” normalizations (FF). The latter was done by scaling core and bypass cross-sectional areas in the initial temporal jet by Vc/Vflt and Vb/Vflt respectively while keeping initial core and bypass velocities (Vc, Vb) as well as densities and tracer concentrations unchanged. This imposes the accumulated axial compression of the spatial jet into the initial conditions of the temporal simulation to leading order.
The spatial simulation results in Fig. A1 generally fall in between the NF and FF temporal results. The geometry of plume dilution differs somewhat: in the temporal cases there is a one-to-one relation between dilution and lateral jet spread, while in the spatial case it is altered by the axial compression. For the microphysics of interest, however, what matters most is not the geometry but the relation between parcel lifetime and dilution. Fig. A1c shows quite respectable agreement in this measure between the spatial and temporal results (particularly the NF results, but moving toward the FF results with time). More detailed analysis also showed that the age differences in the spatial simulation at any fixed distance downstream were relatively modest (seen somewhat in the difference between the black and blue lines in Fig. A1b). Overall we concluded that the NF temporal approximation was adequate for our purposes. What error it produces acts to modestly increase dilution rates and hence potentially EIiceno (in the same manner as a decrease in engine size).
In Fig. A2 we compare the jet spread rate for the spatial LES with that from an axisymmetric RANS simulation of the jet. The latter simulation is that discussed earlier, provided to us by Boeing for use in initializing our temporal LES, continued here far downstream of the engine-exit plane. The agreement between the two provides a useful consistency check on our use of the variable-density, low-Mach-number approximation for the jet LES since the RANS was run fully compressible.
c. Initialization details
To initialize with environmental pe and Te differing from the base values and with which the jet RANS was performed, we followed common industry practice in estimating gas turbine parameters by preserving temperature and pressure ratios, Mach numbers, and η (see, e.g., Volponi 1998; Kurzke 2003). This requires scaling temperatures by , densities by , velocities by and heat release ffQ by .
To simulate the leading effects of scaling the size of an otherwise comparable engine (i.e., same bypass ratio, η and cruise conditions) by a factor sE we assume the core flow remains choked in the engine nozzle; Bernoulli’s law holds in the bypass flow; and Mach numbers and temperatures within the engine are limited by other factors (e.g., avoiding supersonic flows within the engine, and temperature limitations of engine materials). This implies that at least approximately jet dimensions scale by sE, but velocities, pressures and temperatures remain unchanged. An inspection of mass, momentum and total enthalpy conservation then implies that η will remain unchanged while thrust and fuel flow will both scale as (results that are at least broadly consistent with publicly available numbers for different sized engines with similar bypass ratios at “cruise” conditions). We have of course ignored here any more subtle differences between specific engines as well as differences during a flight as aircraft weight drops.
Default grids and domains used in the simulations are given in Table 1. The aircraft was chosen with a wingspan of 28.88 m. Using the results of Fig. A1b (modestly extrapolated) five spans downstream corresponds to an exhaust weighted mean age of 0.53 s, hence the choice in Table 1. The initial downstream domain length corresponds to ~9 initial jet core radii. Simulations with downstream domain length of only ~6 initial jet core radii [chosen as in Paoli et al. (2003) based on expectations for the most unstable modes in a round turbulent jet] were found to suppress the dilution rate after ~0.2 s by inhibiting larger-scale eddies; the choice in Table 1 successfully reproduced mean dilution histories from simulations with much greater domain lengths (24 radii) out to 0.7 s. As well as adding wake vortex fields when passing from grid 1 to grid 2, the jet results were axially compressed (with a corresponding dilation in the cross-stream directions to preserve volume) to move from a temporal approximation appropriate in the near field to one appropriate in the far field. Including this compression (NF to FF) all at the regrid rather than gradually as would be true in a spatial simulation can be expected to modestly overpredict mixing rates approaching the regrid and underpredict for a time afterward (i.e., likely exaggerating the inflection point seen in Fig. 1b somewhat). The compression factor, Vc/Vflt = 0.538, and periodic extension of the domain by 5 multiples leads to the somewhat odd domain length for grid 2. The Table 1 choices in the wake-vortex regime (for both fully 3D and Q3D simulation) are analogous to those given in Lewellen et al. (2014) for a 47.25-m-wingspan B-767 but scaled down by the ratio of wingspans.
The flow through a turbofan aircraft engine consists of a “core” flow where the combustion takes place (with its accompanying exhaust products) and a “bypass” flow ducted around the core flow.
Together with the assumption of ∇ ⋅ u = 0, a small expansion angle also implies small downstream gradients of conserved quantities.
We intentionally initialize the LES in the very near field where the shear layers are strong (and the jet not fully turbulent) in order to promote the rapid establishment of resolved turbulence in the LES. This efficiently produces LES fields that match the RANS mixing rates reasonably well (cf. Fig. A2 in the appendix) but via resolved turbulence satisfying the Navier–Stokes equations.
A pair of wake vortices forms downstream of the aircraft as a consequence of the lift imparted by the wings. For the cases considered here this significantly shears the jets by ~2 s, wrapping them around the vortex cores by ~6 s, and forcing a descent of the primary exhaust plume at an initial rate of ~1.5 m s−1 before the effects of vortex interactions and ambient stratification break up the organized vortex system by ~200 s.
Scaling time and plume cross-sectional area by sE and , respectively, in the former case, and time by in the latter. The sE scaling is approximately correct even including wake-vortex effects, assuming similar aircraft (e.g., engine placement, lift-to-drag ratio) and wingspan scaling approximately as sE.
It is only after ~10 s that the mean plume descent exceeds the plume dimension itself in the present cases, and only after ~30 s that the adiabatic heating in the primary wake approaches 0.5 K.
As a result, if a contrail forms predominantly on ambient aerosol (e.g., for a future aircraft with lower exhaust aerosol numbers) then for some conditions it can be expected to become optically visible only many wingspans behind the aircraft.
Dimensionally, for a free jet of radius rj and velocity Vj, the mixing time scales as τj ~ rj/Vj. Note that increasing Vj (e.g., by increasing thrust) has an opposite effect: reducing τj and hence potentially increasing EIiceno.
Here . For the accretion of a fixed ice mass the mean crystal radius will scale as giving then .
The latter result should perhaps not be surprising even from the box-model perspective: in parcels at the plume edges in the young jet where UFA are contributing significantly to EIiceno in the LES, all the soot is indeed forming ice, but there is not sufficient ice numbers so created to quench the subsequent ice formation on a fraction of UFA in those parcels, i.e., locally (as opposed to integrated over the full plume) conditions are not in this sense “soot rich.”
Kleine et al. (2018) also conclude that their analysis strongly suggests that ice nucleation is solely on soot in a soot-rich regime, conflicting with our findings here given large EIUFA. Their data analysis, however, includes a large correction to avoid double counting of soot particles that assumes that all the ice crystals form on soot cores [citing the box-model results of Kärcher and Yu (2009) as justification]. Accounting for this bias, the scatter in their observations is large enough that fractions of ice formation on cores other than soot at the level seen, e.g., in our Fig. 3 here are in fact not ruled out.
They also show differences between the two inboard engines on the NASA DC-8 by as much as a factor of 2 in the number of particles they produce despite being the same model engine run under nominally the same conditions.
Both ice sedimentation and latent heat from droplet condensation can be safely ignored for the small crystal/droplet sizes and short time scales in the jet.