## Abstract

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 (EI_{iceno}) 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 EI_{iceno} 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 EI_{iceno} and can underestimate the fraction from ultrafine aerosol. The effects of different parameters on EI_{iceno} 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 EI_{iceno}. We find EI_{iceno} 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.

## 1. Introduction

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 *N*_{ice} 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 *N*_{ice} 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 (RH_{i}) 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 EI_{soot} ~ 2 × 10^{14} to 2 × 10^{15} 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, EI_{UFA} ~ 10^{16} to more than 10^{17} (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 EI_{soot}, EI_{UFA} and ambient aerosol may impact the effective emission index for ice number EI_{iceno}.

It is known that the mean plume mixing rate is important in determining EI_{iceno}, 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 EI_{iceno}. This extends our prior LES treatment, where values of EI_{iceno} 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 EI_{soot} 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 *p*_{e} 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” *M*_{c}—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, *M*_{c} = 0.145 in the core-bypass^{1} shear layer, and *M*_{c} = 0.218 in the bypass-ambient shear layer initially, dropping with time as the exhaust jet mixes. At these levels of *M*_{c} 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 $Mc2$, 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 angle^{2} (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 ~*V*_{flt}/*V*_{a}(*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.

### b. Initialization

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 *f*_{f} = 0.32 kg s^{−1}, specific heat of combustion *Q* = 42.9 MJ kg^{−1} and emission index for water $EIH2O=1.25\u2009kg\u2009H2O\u2009per\u2009kg\u2009fuel$. 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 $\theta \u2261T\u2061(pe/p)Ra/cp$, where *c*_{p} is the isobaric specific heat and *R*_{a} the gas constant for air; “potential” density *ρ*_{s} ≡ *p*_{e}/*R*_{a}*θ* and “potential” velocity **u**_{s} ≡ (*ρ*/*ρ*_{s})**u**] are free from these fluctuations. We then use the RANS *θ*, *ρ*_{s}, and **u**_{s} fields 1 m downstream from the engine exit plane^{3} 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 $p*$, $T*$, 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 EI_{iceno} 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 *p*_{e} and *T*_{e} differing from the base values $p*$ and $T*$ or with an engine size scaled by a factor *s*_{E} are described in the appendix.

For longer simulation times the wake vortices influence the mixing^{4} 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).

### c. Microphysics

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

with

for droplet radius *r*, water activity *a*_{w}, Kelvin parameter for solution droplets $aks$, solution–air interface surface tension *σ*_{s/a}, water vapor gas constant *R*_{υ}, and liquid water density *ρ*_{w}. Because of competing solute and curvature effects *S*_{wd} will have a maximum at some critical radius *r**. The critical saturation ratio for aerosol activation, $Sw*$, is the value of *S*_{wd} at *r* = *r**. The aerosol properties are defined in this approximation solely by the dry aerosol core radius *r*_{d} 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 *r*_{d} = 20 nm; an ambient CCN with *r*_{d} = 20 nm and *κ* = 0.9 [consistent with the value for H_{2}SO_{4} given in Petters and Kreidenweis (2007)]; and an ultrafine aerosol (UFA) with *r*_{d} = 2 nm and *κ* = 0.9. For reference the $Sw*$ at 218.8 K for these three choices is respectively 1.065, 1.0099, and 1.347. For our base conditions quite large *S*_{w} are potentially reached, as shown in Fig. 1a; thus even quite small aerosol with $Sw*$ well above 1 can potentially activate, grow, and freeze. The $Sw*$ for the nominal “CCN” is sufficiently close to one that further increase in *r*_{d} 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 EI_{iceno} 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 EI_{iceno} 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 *V*_{w} to freeze in time interval *dt* is 1 − exp(−*J*_{hom}*V*_{w}*dt*). For *J*_{hom} we use the water activity dependent parameterization of Koop et al. (2000) with *a*_{w} given by (2). Riechers et al. (2013) argue that values for *J*_{hom} 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 *f*_{JH} to *J*_{hom} in the LES and use *f*_{JH} = 0.01 as a default value; extensive sets of runs with *f*_{JH} = 1 were also performed to assess sensitivities to homogeneous freezing rate in different regimes. Generally we have found that sensitivity to be low: *J*_{hom} increases extremely rapidly with decreasing *T* so even a large change in *f*_{JH} 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 EI_{iceno}.

### 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 *C*_{E} 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 $VE$ follows by demanding $\u222bdxCE\u2061(x,t)=\u2329CE\u2061(t)\u232aEVE\u2061(t)$, 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 *C*_{E}), ⟨*ϕ*⟩_{E} and $VE$ 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 (*s*_{E}) or *T*_{e}; 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 (*T*_{e}, RH_{i}, *p*_{e}, stratification *dθ*/*dz*), aerosol (*r*_{d}, *κ*, EI_{soot}, EI_{UFA}, or combination EI_{aero} ≡ EI_{soot} + EI_{UFA}, environmental number density *N*_{e}), engine size scaling (*s*_{E}), 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 EI_{iceno}. Where unspecified, base conditions were used of *s*_{E} = 1, *p*_{e} = 238.4 hPa, *T*_{e} = 218.8 K, RH_{i} = 110%, *dθ*/*dz* = 2.5 K km^{−1}, and *f*_{JH} = 0.01.

### a. Exhaust aerosol

Figure 2 shows sample EI_{iceno} evolution beginning from soot aerosol, varying the dry-size specifications, turbulence realization, and model treatment. A relatively large EI_{soot} 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 EI_{iceno} 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 EI_{iceno} 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 EI_{soot} to EI_{iceno}. 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 EI_{iceno} 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 EI_{iceno} production than the competition between different sized cores in the bimodal dry-core-size distribution provides. Figure 3 shows a related example. The eventual EI_{iceno} 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 EI_{aero} 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 *T*_{e} = 225 K EI_{iceno} 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 RH_{i} = 110% or 130%, respectively). Second, away from threshold conditions, all EI_{aero} is converted to EI_{iceno} if EI_{aero} is small enough, irrespective of the details of aerosol type or ambient conditions. As EI_{aero} 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 *T*_{e} (at fixed RH_{i}), consistent with the greater supersaturation levels in the plume that are present (cf. Fig. 1a). The effect of ambient RH_{i} 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 EI_{iceno} from the box model relative to the corresponding LES.

In Fig. 5a as EI_{aero} increases, EI_{iceno} appears to saturate toward a *T*_{e} 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 EI_{iceno} 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 EI_{iceno} 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 EI_{iceno} at the end of the jet-dominated regime (~10 s) regardless of how large EI_{iceno} 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 EI_{iceno} 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 *T*_{e}, RH_{i}, 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 EI_{iceno} evolution for cases with only ambient aerosol. As in the cases with exhaust aerosol (and for the same reasons), colder *T*_{e} leads to greater EI_{iceno}. 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 *S*_{w} > 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 EI_{iceno} 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 EI_{iceno} on *N*_{e} is more clearly illustrated by plotting the ratio EI_{iceno}/*N*_{e} as in Fig. 8. Values of *N*_{e} above and below typical expectations in the upper troposphere are included to clarify the behavior in these limits. EI_{iceno}(*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 EI_{iceno}/*N*_{e} that is just a measure of evolving plume volume, set purely by plume dilution. With increasing *N*_{e}, EI_{iceno}/*N*_{e} breaks off from this limit at progressively earlier times as moisture within the plume is more rapidly consumed by ice growth. For small enough *N*_{e}, EI_{iceno}/*N*_{e} breaks from this limit purely due to dilution reducing *S*_{w} within the plume (cf. Fig. 1) rather than ice growth; below this point the ratio becomes independent of *N*_{e}.

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 EI_{iceno}. Given the success of the box model it is possible to more thoroughly explore EI_{iceno} production from ambient aerosol over a large parameter space. Figure 10 provides examples varying *N*_{e}, *T*_{e}, altitude and aerosol size. Some of the dependencies seen, particularly in the low and high *N*_{e} limits, can be understood from simpler models considered in section 4 below.

### c. Exhaust plus ambient aerosol

If EI_{aero} 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 EI_{iceno}. On the other hand, if EI_{aero} 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 EI_{iceno} contribution from the exhaust aerosol for simple reasons (for low-enough EI_{aero}, 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 *N*_{e} and modest EI_{aero} present the resulting EI_{iceno} will lie below the sum of the EI_{iceno} produced by each if present in isolation, but generally only moderately so. There can be exceptions, however, such as for *T*_{e} = 205 K and *N*_{e} ~ 10^{8} m^{−3} in Fig. 12. There the addition of EI_{aero} of 10^{13} kg^{−1} leads to less EI_{iceno} than the ambient aerosol would have produced alone. At this low *T*_{e} and modest *N*_{e} the EI_{iceno} production is spread over ~10 s (cf. Fig. 7); over this long a time the early ice production from the addition of EI_{aero} 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 *s*_{E} as discussed in the appendix. For low-enough EI_{aero} all the aerosol converts to ice and EI_{iceno} is independent of *s*_{E}. For larger EI_{aero}, the EI_{iceno} after formation falls with increasing engine size (*s*_{E}). 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 EI_{iceno}. Figure 14a shows related results from box-model simulations of EI_{iceno} production from ambient aerosol (where the box-model results have proved trustworthy). Again for larger *N*_{e} a larger engine produces a smaller EI_{iceno} than a smaller (but otherwise comparable) engine. Further there is a clear similarity in the effects on the falloff of EI_{iceno}/*N*_{e} in Fig. 14a due to a decrease in engine size versus due to a decrease in *N*_{e} (curves for different *s*_{E} have similar shapes, differing by a shift in *N*_{e}). 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 *s*_{E} and *N*_{e} while keeping other factors (*T*, RH_{i}, *κ*, etc.) constant, *τ*_{mix} scales as *s*_{E} and *τ*_{fc} approximately as $Ne\u22122/3$.^{9} If the dependence of EI_{iceno}/*N*_{e} on *τ*_{mix} and *τ*_{fc} is only through their ratio (scaling as $sENe2/3$) then plotting the results against any function of this ratio [e.g., $\u2061(\tau mix/\tau fc)3/2~sE3/2Ne$] should collapse the *s*_{E} 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 EI_{iceno}, 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 EI_{iceno} production as follows. We assume that aerosol nucleates ice if and only if it activates (i.e., encounters $Sw>Sw*$). Consider an infinitesimal plume parcel following the products from burning a mass of fuel *f*_{p}, diluting in an environment with total water mixing ratio *q*_{te} (assuming no condensate there). The evolving mass of water in the parcel (*M*_{wp}) is related to the growing parcel air mass (*M*_{ap}) as

Combustion provides an increase in enthalpy in the parcel that is conserved thereafter, *δH*_{p} ≡ *Q*(1 − *η*)*f*_{p} 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 *ϵ* ≡ *R*_{a}/*R*_{υ} and *e*_{s}(*T*) the saturation vapor pressure over liquid water at temperature *T*. Setting $Swp=Sw*$ 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 $\u2061(T\u02dc)$,

where

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 $T\u02dc$, corresponding to when the parcel first reaches $Swp=Sw*$ and when it last does so, which we denote *T*_{p1} and *T*_{p2}, respectively. In the LND limit under our assumptions the number of ice crystals produced from ambient aerosol in the parcel will be given by *N*_{e}*M*_{ap}/*ρ*_{e} with *M*_{ap} given by (7) with *T*_{p} = *T*_{p2}. 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 *T*_{p2} replaced by *T*_{p1} is in principle relevant in the opposite limit, i.e., *N*_{e} large enough and mixing time scales long enough that moisture depletion sufficient to cut off activation occurs before *T*_{p} dilutes much below *T*_{p1}. 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 *N*_{e} considered; nonetheless, as discussed below, it explains some behavior encountered for large *N*_{e}.

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 *e*_{s}(*T*) given by Buck (1981) about *T* = *T*_{e} and solving; after some algebra,

where *S*_{we} is the environmental water saturation ratio. In the term in curly brackets *T*_{e} 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 *r*_{d} and *κ* is solely through $Sw*$. 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 *p*_{e}.

For cases with exhaust aerosol alone the Eq. (11) LND prediction is simply EI_{iceno} = EI_{aero}, already seen in the LES results (e.g., the green lines in Figs. 5 and 13). More nontrivially, much of the results for EI_{iceno} 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 EI_{iceno}/*N*_{e} for low *N*_{e} and their dependence on *T*_{e} and *r*_{d} 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 “*T*_{p1}” branch of the solution in Fig. 15 illuminates the behavior of the simulation results for large *N*_{e}. In particular the switchover in Fig. 10a from EI_{iceno}/*N*_{e} increasing with decreasing *T*_{e} at low *N*_{e} to increasing with increasing *T*_{e} at high *N*_{e} is predicted by the corresponding switch in behaviors seen in the “*T*_{p2}” and “*T*_{p1}” branches in Fig. 15. The LND limit model is not in exact agreement with the simulations even for low *N*_{e}. At colder temperatures for larger aerosol the LND limit model tends to slightly underpredict EI_{iceno}/*N*_{e} from the simulations (up to 15% below has been seen), and the simulations show a small dependence of EI_{iceno}/*N*_{e} on plume mixing time scale even for low *N*_{e} 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 $Sw>Sw*$ 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 RH_{i} 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 *S*_{w} 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 EI_{iceno} to only a fraction of that implied by the activatable aerosol present. We will not try to formulate an analytic estimate for EI_{iceno} 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 *t*_{0}, 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 $aki$ the analog for ice of $aks$ in (2), *q*_{si} the ice saturation mixing ratio, *D*_{υ} the water vapor diffusion constant, *M*_{E}(*t*) the equilibrium ice mass for unit mass of fuel burn, and Λ a parameter depending on the growth rate of *M*_{E}(*t*). For *t* ≫ *t*_{0} the plume *T* will have relaxed to good approximation to *T*_{e} and *M*_{E}(*t*) will be given by

with *A*(*t*) the plume cross-sectional area and *f*_{m} the fuel mass per length of exhaust plume. To get an upper limit we assume arbitrarily large initial EI_{iceno}(*t*_{0}), so that $m\xaf\u2061(t0)$ 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*)/*f*_{m} 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 *T*_{e} 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 *M*_{E}(*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 EI_{iceno} 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 EI_{soot} ~ 10^{15} or less away from near-threshold conditions; the only partial conversion of soot to ice for much larger EI_{soot} with the fraction decreasing with increasing temperature; the rising (and significant) contribution of ambient aerosol to EI_{iceno} with decreasing temperature and EI_{soot}), but not all (e.g., their conclusion that UFA do not contribute to EI_{iceno} in the soot-rich regime of EI_{soot} ~ 10^{15} or above). We have shown that well below the contrail threshold temperature and for sufficiently low numbers of activatable exhaust or ambient aerosol, EI_{iceno} 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 EI_{iceno} 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 EI_{iceno}. Even for larger aerosol numbers we have found that LES results for EI_{iceno} are approximated well by a box-model representation if EI_{iceno} 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 EI_{iceno} relative to the LES results and an underrepresentation of the fraction of EI_{iceno} 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 EI_{soot} alone may not significantly change EI_{iceno}: 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 EI_{iceno} 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 EI_{iceno} determinations of a few times 10^{15} kg^{−1} (their Fig. 8), consistent with expectations here for aged contrails in Fig. 5b. More recent measurements of Kleine et al. (2018) (for *T*_{e} = 215.5 K and a comparable aircraft to what we consider here) suggest an EI_{soot} ~ 5 × 10^{15} and EI_{iceno} of ~3 × 10^{15} initially after formation, dropping to ~1 × 10^{15} with age as the wake descends, consistent with expected ranges here in Fig. 5 at *T*_{e} = 215.5 K and EI_{aero} = 5 × 10^{15}.^{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 EI_{iceno}(*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 EI_{iceno} 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 EI_{soot} and EI_{iceno} 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 EI_{soot} in the range measured for current jetliners (~2 × 10^{14} to 2 × 10^{15} kg^{−1} fuel) EI_{iceno} 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 EI_{soot} and EI_{UFA} has the potential to reduce EI_{iceno} from typical current levels (~10^{14}–10^{15} 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 EI_{iceno} 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 × 10^{4} cm^{−3} (cf. their Fig. 5). This motivated extending our sensitivity studies to include quite large values for *N*_{e}. The resulting EI_{iceno} 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 RH_{i} 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 EI_{iceno} 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 EI_{iceno} are robust, the actual values found for EI_{iceno} 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 EI_{iceno} from turbulence realizations and differences between spatial and temporal simulations can be significant.

## Acknowledgments

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.

### APPENDIX

#### Model Details

##### a. Variable density implementation

In the expansion of the jet at constant pressure the specific enthalpy *h*_{t} ≡ *c*_{p}*T* + **u**^{2}/2 − *L*_{s}*q*_{i} is conserved (where *c*_{p} is the isobaric specific heat, *L*_{s} is the latent heat of sublimation of ice and *q*_{i} is the ice mixing ratio).^{A1} We define our density variable in terms of *h*_{t}, $\rho ^\u2261cpp/Raht$ (with *R*_{a} 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 EI_{iceno} when it is mainly produced early in the jet evolution. Using $\rho ^$ 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 *C*_{E} an “aging tracer” *C*_{a} 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 *C*_{E} as a local source term, ∂*C*_{a}/∂*t* = [∂*C*_{a}/∂*t*]_{FM} + *C*_{E}. In a temporal simulation starting with *C*_{a} = 0, the local ratio *C*_{a}/*C*_{E} 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 $t\xafage\u2061(y)=\u222b\u200bCa\u2061(x,t)dxdzdt/\u222b\u200bCE\u2061(x,t)dxdzdt$.

The simulation was run without *C*_{a} for 2 s to reach a quasi-steady state, then for 0.5 s for *C*_{a} 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 *V*_{c}/*V*_{flt} and *V*_{b}/*V*_{flt} respectively while keeping initial core and bypass velocities (*V*_{c}, *V*_{b}) 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 EI_{iceno} (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 *p*_{e} and *T*_{e} differing from the base values $p*$ and $T*$ 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 $\u2061(Te/T*)$, densities by $\u2061(pe/p*)\u2061(T*/Te)$, velocities by $Te/T*$ and heat release *f*_{f}*Q* by $\u2061(pe/p*)Te/T*$.

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 *s*_{E} 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 *s*_{E}, 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 $sE2$ (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, *V*_{c}/*V*_{flt} = 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.

## REFERENCES

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*J. Comput. Phys.*

*J. Geophys. Res. Atmos.*

*J. Appl. Meteor.*

*Atmos. Environ.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Nat. Commun.*

*Geophys. Res. Lett.*

*J. Atmos. Sci.*

*J. Geophys. Res. Atmos.*

*Meteor. Z.*

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*Nature*

*Annu. Rev. Phys. Chem.*

*Proc. Turbo Expo*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*AIAA J.*

*J. Atmos. Sci.*

_{x}and HO

_{x}wake chemistry

*J. Geophys. Res.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Atmos. Chem. Phys.*

*Nature*

*Annu. Rev. Fluid Mech.*

*Phys. Fluids*

*J. Fluid Mech.*

*Phys. Fluids*

*J. Fluid Mech.*

*Atmos. Chem. Phys.*

*Phys. Chem. Chem. Phys.*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Meteor. Z.*

*Ice Formation and Evolution in Clouds and Precipitation: Measurement and Modeling Challenges*,

*Meteor. Monogr*.

*Atmos. Chem. Phys.*

*J. Geophys. Res.*

*Atmos. Chem. Phys.*

*Proc. Int. Gas Turbine and Aeroengine Congress*

*Atmos. Chem. Phys.*

*Geophys. Res. Lett.*

*Chem. Rev.*

## Footnotes

^{1}

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.

^{2}

Together with the assumption of ∇ ⋅ **u** = 0, a small expansion angle also implies small downstream gradients of conserved quantities.

^{3}

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.

^{4}

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.

^{5}

Scaling time and plume cross-sectional area by *s*_{E} and $sE2$, respectively, in the former case, and time by $T*/Te$ in the latter. The *s*_{E} 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 *s*_{E}.

^{6}

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.

^{7}

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.

^{8}

Dimensionally, for a free jet of radius *r*_{j} and velocity *V*_{j}, the mixing time scales as *τ*_{j} ~ *r*_{j}/*V*_{j}. Note that increasing *V*_{j} (e.g., by increasing thrust) has an opposite effect: reducing *τ*_{j} and hence potentially increasing EI_{iceno}.

^{9}

Here $\tau fc~(r\xafNe)\u22121$. For the accretion of a fixed ice mass the mean crystal radius will scale as $r\xaf~Ne\u22121/3$ giving then $\tau fc~Ne\u22122/3$.

^{10}

The “slow growth” regime generally requires the crystal population to be in a quasi-equilibrium with local conditions (which may be changing) and a total ice mass in the parcel growing no faster than linearly in time (Lewellen 2012).

^{11}

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 EI_{iceno} 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.”

^{12}

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 EI_{UFA}. 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.

^{13}

See, e.g., Lewellen and Lewellen (2001b) for a related example of the difficulties of interpreting in situ sampling in this plume regime.

^{14}

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.

^{A1}

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.