Box model trajectory studies of contrail formation using a particle-based cloud microphysics scheme

We investigate the microphysics of contrail formation behind commercial aircraft by means of the particle-based LCM (Lagrangian Cloud Module) box model. We extend the original LCM to cover the basic pathway of contrail formation of soot particles being activated into liquid droplets that soon after freeze into ice crystals. In our particle-based microphysical approach, simulation particles are used to represent different particle types (soot, droplets, ice crystals) and properties (mass/radius, num5 ber). The box model is applied in two frameworks. In the classical framework, we prescribe the dilution along one average trajectory in a single box model run. In the second framework, we perform a large ensemble of box model runs using 25000 different trajectories inside an expanding exhaust jet as simulated by the LES (large-eddy simulation) model FLUDILES. In the ensemble runs, we see a strong radial dependence of the temperature and relative humidity evolution. Droplet formation on soot particles happens first near the plume edge and a few tenths of seconds later in the plume centre. Averaging over the 10 ensemble runs, the number of formed droplets/ice crystals increases more smoothly over time than for the single box model run with the average dilution. Consistent with previous studies, contrail ice crystal number varies strongly with atmospheric parameters like temperature and relative humidity near the contrail formation threshold. Close to this threshold, the freezing fraction of soot particles depends strongly on the geometric-mean dry core radius and the hygroscopicity parameter of soot particles. This sensitivity is quite low 15 at ambient conditions far away from the formation threshold. Absolute ice crystal numbers, on the other hand, are controlled by the soot number emission index for all atmospheric conditions. The comparison with a recent contrail formation study by Lewellen (2020) (using similar microphysics) shows a later onset of our contrail formation due to a weaker prescribed plume dilution. If we use the same dilution data, our and Lewellen’s evolution in contrail ice nucleation show an excellent agreement cross-validating both microphysics implementations. This 20 means that differences in contrail properties mainly result from different representations of the plume mixing and not from the microphysical modelling. The presented aerosol and microphysics scheme describing contrail formation is of intermediate complexity and thus suited to be incorporated in an LES model for 3D contrail formation studies explicitly simulating the jet expansion. The presented box model results will help interpreting the upcoming, more complex 3D results. 25 1 https://doi.org/10.5194/acp-2021-361 Preprint. Discussion started: 21 May 2021 c © Author(s) 2021. CC BY 4.0 License.

in natural cirrus clouds. The present paper describes how the LCM model has been extended in order to cover the specifics of contrail formation and is a first step towards the goals of ConForm. Our main objective is to test the extended LCM in a simple dynamical framework by running it in a box model setup. The microphysical parameterisation is similar to that in Lewellen (2020) with the further improvement of prescribing soot particles by a lognormal size distribution. Moreover, the numerical approach of the microphysics differs as our study relies on a particle-based description and not on an Eulerian spectral bin 80 model. Paoli et al. (2008) and Vancassel et al. (2014) extended the box model framework to a multi-0D approach. Multi-0D here means that the box model was run for a large ensemble of trajectories that sample the expanding jet plume as provided by 3D LES and that represent many different plume dilution evolution. In such an "offline" approach, the microphysics can be more advanced than in a full 3D framework with "online" microphysics. In our study, we intend to highlight the importance 85 of the huge variability in thermodynamic and contrail properties resulting from the variability in the inhomogeneous turbulent mixed jet plume. Therefore, we will use an ensemble of many trajectories where the data have been taken from Vancassel et al. (2014). We expect that the temporal evolution of ensemble mean properties leads to improved scientific results compared to those using one single average trajectory. Paoli et al. (2008) mainly studied the temporal evolution of concentrations and (apparent) emission indices of volatile particles 90 and ice crystals depending on the fuel sulphur content. Vancassel et al. (2014) revealed the spatial variability in thermodynamic quantities and contrail properties. Moreover, they used an online method with less detailed and an offline method with more detailed contrail microphysics and highlighted differences between both methods. The two studies have the following improvements compared to the online 3D studies of Khou et al. (2015Khou et al. ( , 2017: They include different types of aerosol particles, the chemical processing of fuel sulphur and the formation and growth of volatile particles. Curvature and solution effect are 95 accounted for. A further improvement is prescribing lognormally size-distributed soot particles using 60 size bins and varying the geometric mean dry core radius between 10 and 30 nm. For conventional aircraft engines, soot particles are in general the major source for the formation of contrail ice crystals (e.g., Kärcher et al., 1996;Kärcher and Yu, 2009;Kleine et al., 2018). The number of emitted soot particles significantly influences the contrail cirrus climate impact . Since engine soot emissions are quite variable and uncertain particu-100 larly in terms of their number and size (e.g., Anderson et al., 1998;Petzold et al., 1999;Schumann et al., 2002;Agarwal et al., 2019), our main focus is to investigate the sensitivity of contrail properties (during their formation stage in the first seconds) on different soot properties. We will also analyse the variability of contrail ice nucleation with ambient conditions as already explored by Bier and Burkhardt (2019) using the analytical parameterisation of Kärcher et al. (2015).

105
This section gives a short overview of the basic "Lagrangian Cloud Module" (LCM)-model (Sect. 2.1), which has been extensively used for simulations of natural cirrus as well as young and aged contrails in the recent past. Section 2.2 then describes in more detail the modifications and extensions for simulating contrail formation and Sect. 2.3 gives an overview about the general plume evolution and our trajectory data. Finally, we will explain the box model framework, in which LCM is employed in the present study.

Particle-based LCM microphysics
The particle-based ice microphysics model LCM (Sölch and Kärcher, 2010) comprises explicit aerosol and ice microphysics and is employed for the simulation of pure ice clouds like natural cirrus (e.g., Sölch and Kärcher, 2011) or contrails (e.g., Unterstrasser, 2014;Unterstrasser and Görsch, 2014;Unterstrasser et al., 2017a, b). In a particle-based microphysical approach, hydrometeors (like ice crystals or water droplets) are described by simulation particles (SIPs) unlike traditional Eulerian ap- The SIP data structure is augmented and stores information about the particle type (soot, droplet, ice crystal) and properties (number, mass/radius, freezing temperature).

Microphysics of contrail formation 2.2.1 Exhaust particles
In this study, we consider soot as the only exhaust particle. Background particles will become relevant for aircraft with soot-130 poor or even soot-free emissions as for liquid hydrogen propulsion (Kärcher et al., 2015;Rojo et al., 2015). We neglect ultrafine volatile particles since this would require the inclusion of complex chemical processes like the recombination between ionclusters leading to huge computational costs. These ultrafine particles are in particular relevant at ambient temperatures more than 10 K below the contrail formation threshold and for engines with soot-poor emissions still containing sulphuric and/or organic compounds (Kärcher and Yu, 2009). The size spectrum of soot particles ranges from a few to hundreds of nanometer 135 with geometric-mean values of around 15 nm (Petzold et al., 1999). Soot particles are composed of several nanometer-sized primary spherules combined in aggregates. For simplicity, we assume a spherical shape. Recent studies indicate a fractal irregular shape of the larger (more than several 10 nm) soot particles (e.g., Wang et al., 2019;Distaso et al., 2020). Taking this into account would modify their effective surface area and, therefore, the activation into water droplets and condensational growth rates. But it would require additional complexity that may not be relevant in such an intermediate detailed microphysical We prescribe lognormally size distributed soot particles by an ensemble of N SIP simulation particles (SIPs) using a technique described in Appendix A of Unterstrasser and Sölch (2014). In a priory tests, we performed simulations with different values of N SIP ranging from 50 to 200. We find converged results of the analysed contrail properties for N SIP above around 100. Since this threshold value slightly varies with the geometric-mean dry core radius and particle size distribution width, we choose a 145 default value of around N SIP = 130.

Basic formation pathway
In general, the moist exhaust plume cools due to mixing with ambient air. Contrails can form if ambient temperature is below the so called Schmidt-Appleman (SA)-threshold temperature (Θ G ) and the plume becomes water-supersaturated in a certain temperature range. The threshold temperature varies with atmospheric parameters (ambient relative humidity over water, air 150 pressure) and engine/fuel properties (heat combustion, propulsion efficiency, water vapour emission index). The calculation of Θ G is described in Appendix 2 of Schumann (1996).
We use the Kappa-Köhler theory (Sect. 2.2.3) to calculate which soot particles are able to activate into water droplets. Thereby, we prescribe the hygroscopicity parameter of soot particles as a measure for their solubility. This is a more convenient approach than defining a fixed fuel sulphur content since other polar species may lead to active sites on soot particles as well. Smaller 155 soot particles require higher water-supersaturation to overcome the Köhler barrier than larger ones (Fig. B1). Therefore, the latter particles can activate earlier and may suppress the droplet formation on smaller particles due to depletion of water vapour.
The super-cooled droplets subsequently grow by condensation as long as the exhaust air remains water-supersaturated. They freeze into contrail ice crystals when the homogeneous freezing temperature (Sect. 2.2.5) is reached. Larger droplets can freeze earlier (i. e., at higher plume temperature) and grow further by deposition before smaller droplets manage to freeze into ice 160 crystals.

Koehler theory
To calculate the equilibrium saturation ratio over a droplet or ice crystal surface (S K ), we use the Kappa-Köhler equation (Petters and Kreidenweis, 2007)

165
consisting of the solution term, described by the activity of water a w , and the exponential Kelvin term. r d is the particle dry core radius, r the droplet radius and κ the hygroscopicity parameter characterising the solubility of the aerosol particle; σ is the surface tension of the solution droplet, ρ w the mass density of water and T temperature. M w and R denote the molar mass of water and the universal gas constant, respectively. For the activation of exhaust particles into water droplets, the plume saturation ratio needs to overcome the maximum of S K , denominated by the critical saturation ratio (S c ). This maximum 170 increases with decreasing r d due to the Kelvin effect ( Fig. B1 a)) and decreases with increasing κ due to the solution effect ( Fig. B1 b)). The calculation of S c is described in Appendix B.
The equilibrium saturation vapour pressures over the droplet (ice crystal) surface, used in Eqs.
(2) and (7), are the product of the saturation vapour pressure over a flat water (ice) surface and S K .

Condensational droplet growth 175
The droplet growth of activated soot particles is calculated according to Barrett and Clement (1988) and Kulmala (1993). In the steady state (small changes in the vapour flux) and for spherical droplets, the single droplet mass growth equation is given by where e v is the partial vapour pressure and e K,w equilibrium saturation vapour pressure over the droplet surface. L c denotes the 180 specific latent heat for condensation/evaporation, D v the binary diffusion coefficient of air and water vapour,K the conductivity of air and R v the specific gas constant of vapour. The terms F M and F H represent the mass and heat diffusion term with β m and β t denoting the transitional correction factors calculated according to Fuchs and Sutugin (1971).

Homogeneous freezing
We calculate the homogeneous freezing temperature (T frz ) of a super-cooled droplet with radius r according to the method of 185 Kärcher et al. (2015). This method is suitable for a strong cooling situation typically occurring in a jet plume. We introduce the nucleation rate as the number of droplets freezing per unit time where J is the nucleation rate coefficient and LW V = 4 3 π(r −r d 3 ) is the liquid volume available for freezing. The temperature dependent nucleation rate coefficient J is parametrised according to Riechers et al. (2013) J/(m 3 s) = 10 6 · exp(a 1 T + a 2 ), where a 1 = 3.574 K −1 and a 2 = 858.72 are empirical constants. Assuming that the whole droplet freezes immediately if ice nucleates within its volume, the pulse-like nature of j may be expressed by a characteristic freezing time scale withṪ denoting the plume cooling rate. Using Eq. (4) we obtain τ frz (a 1Ṫ ) −1 .

195
The homogeneous freezing requirement will be fulfilled if the product of j = LW V ·J and τ frz is unity. Using this requirement, inserting the relations from above and rearranging to T := T frz yields The homogeneous freezing temperature in contrails declines with decreasing droplet radius and increasing cooling rate. As shown in Fig. 8 of Kärcher et al. (2015), T frz ranges between around 230 and 232 K for typical droplet sizes (100-500 nm) at 200 cruise altitude conditions.

Depositional ice crystal growth
The mass change of a single ice crystal (Mason, 1971) is given by where C is the shape factor, r the equivalent volume radius and L d is the specific latent heat for deposition/sublimation. In our 205 study, we set C to unity assuming spherical ice crystals as they are typically observed for young contrails (e.g., Schröder et al., 2000). Since sedimentation is of low relevance during contrail formation, we here neglect ventilation effects. The calculation of the transitional correction factors β v and β k is described in Appendix A.

Plume evolution and trajectory data
In this section, we first explain the plume evolution assuming a mean state. Subsequently, we describe which kind of trajectory 210 data we will use and how dilution, a measure for the air-to-fuel ratio and the volume of the plume, is derived from these data.

General plume dilution equations
The temporal evolution of the mass mixing ratio of a species χ(t) is given by (Kärcher, 1995;Kärcher, 1999) where the first term describes the mixing of the hot exhaust with ambient air and the second term (ξ) represents sink and source 215 terms. The entrainment rate ω χ characterises the dilution and χ a is the ambient mass mixing ratio of the species. Applying the equation above to water vapour and setting up an analogous equation for the evolution of the plume temperature (T ) yieldṡ where T a and x v,a are ambient temperature and vapour mass mixing ratio, respectively. The subscript LH relates to latent 220 heating and P C stands for phase changes from water vapour to liquid water or ice and vice versa.
From the evolution of a passive tracer χ p (t) (defined by Eq. (8) with ξ = 0), we can derive the entrainment rate ω p and set ω v = ω p since water vapour is supposed to be diluted like a passive tracer. The entrainment rate characterising the temperature change, ω T , differs from ω v by a factor called the Lewis number Le = ω v /ω T , which is the ratio of the thermal 225 and mass diffusivity. We set Le to one assuming that heat and mass diffusion proceed equally fast since turbulent exchange processes are dominating over molecular exchange processes. Due to ω v = ω T , we skip the index from now on and simply use ω. The entrainment rate ω describes the instantaneous mixing strength. Integrating it over time, we obtain the plume dilution factor (Kärcher, 1995) 230 where the subscript "0" refers to conditions at the engine exit at time t 0 := 0 s.
At the engine exit, the emitted air is already a mixture of ambient air and combustion products. The initial plume dilution N 0 can be derived from the temperature excess where Q is the specific heat of combustion, η the propulsion efficiency andc p = 1020 J(kg K) −1 an average heat capacity for 235 dry air over a temperature range between around 200 and 600 K (see derivation in Schumann, 1996).
The overall dilutionN generally increases with plume age due to continuous entrainment of ambient air and is related to the dilution factor (this quantity decreases with plume age) viâ implying D 0 = 1.

240
The initial plume area A 0 results from mass conservation (Schumann et al., 1998) where ρ 0 = pa Rd T0 is the air density with p a denoting air pressure, R d the specific gas constant of dry air and m F the fuel consumption (in kg m −1 ). The evolution of the plume area is given by

245
In the numerical implementation, the ordinary differential equations (ODE) like Eqs. (9) and (10) are discretised by a simple Euler forward scheme. Due to the linearity of the ODEs, the implementation of an implicit Euler backward or a second-order trapezoidal scheme is straight-forward, but we find no improvement over the forward scheme. The source terms in the ODEs are provided by the LCM. In the LCM, the computation involves summations over all SIPs with phase changes or water vapour uptake/release during the time step. All model components (plume evolution ODE and LCM microphysics) use a time step of 250 0.001 s. This time step is sufficiently small to account for the fast diffusional growth at high plume supersaturations. Vancassel et al. (2014) used the LES solver FLUDILES to simulate the plume evolution in a 3D domain over 10 s starting behind the engine nozzle of an A340-300 aircraft. They sampled the plume with 25000 trajectories that were randomly positioned inside the initial plume area A 0 = πr 0 2 with radius r 0 = 0.5 m. We assume that all trajectories represent an equal share 255 of the plume. For each trajectory, the temperature evolution T 3D,k (t) has been tracked and we use these data to infer the dilution factor (or equivalently the entrainment rate) along each trajectory k by assuming that temperature is a passive tracer:

Trajectory ensemble data
T 3D,0 = 580 K and T 3D,a = 220 K are the plume exit and ambient temperature of the FLUDILES simulation. These equations are only meaningful, if T 3D,k > T 3D,a is fulfilled. However, the LES temperature occasionally falls below T 3D,a due to a pressure 260 drop resulting from the interaction of the plume with the wake vortex. In the computation of ω k (t) or D k (t), we replace all values of T 3D,k being below T 3D,a + with this lower limit. The threshold is chosen such that the maximum dilution is still reasonable for a few seconds of plume age (D min ≈ 2.5 · 10 −3 ).
In the FLUDILES simulation, the initial plume temperature profile is based on the Crocco-Busemann relation (White, 1974) and introduces a smooth radial transition between the plume and the environment for numerical reasons. This means that 265 T 3D,k (t 0 ) equals T 3D,0 only directly behind the nozzle exit centre and decreases down to around 400 K towards the plume edge (see later in Fig. 1 a)). These trajectories start with D k,0 < 1 and an accordingly higher initial dilutionN k,0 . For simplicity, we mostly neglect this fact in the following description, but the initialisation in this transition region is consistently treated in our model.
FLUDILES temperature is higher than a "passive tracer temperature". Therefore, we underestimate dilution causing a later onset of contrail formation. This effect is partially compensated by the inhomogeneous initial temperature profile implying that exhaust air parcels are already diluted near the plume edge.
Basically, we only infer the dilution from the trajectory data according to Eq. (17) and plug the corresponding entrainment rate values into Eqs. (9) and (10) for each trajectory. For the radial plots presented later on, the coordinate data of these 275 trajectories are used, but the actual computations are independent of this type of information.
Note that the ambient temperature T a and plume temperature at engine exit T 0 in Eqs. (9) and (10) are allowed to differ from T 3D,a and T 3D,0 used in Eq. (17). Similarly, the ambient and initial plume water vapour can be chosen independently of the 3D data. In such a case, we basically use only the information of the dilution and its variability.
Moreover, a single mean state trajectory can be derived from the trajectory ensemble. We determine the corresponding mean 280 dilution D AT by averaging the pre-processed temperature data over all 25000 trajectories and plugging the average temperature evolutionT 3D into Eq. (17).

Box model framework
We use the LCM box model, which is based on the 3D-LCM (introduced in Sect. 2.1) and has been extended for the contrail formation microphysics. Similar to Paoli et al. (2008), we either perform a single box model run with an average dilution 285 (labelled "average traj") or an ensemble of box model runs for the trajectory dilution data. In this case, physical quantities are suitably averaged in the end (labelled "ensemble mean"). The simulated time t sim is 2 s. In the following, we describe the initial conditions (Sect. 2.4.1) and introduce important diagnostics (Sect. 2.4.2). Table 1 summarises initial conditions of our baseline scenario, which are described in the following. The plume exit temperature 290 (near the plume centre) is set to 580 K. We prescribe an ambient temperature of either 220 K (baseline case) or 225 K (nearthreshold case). The baseline SA-threshold temperature (Θ G ) is 226.2 K for the given atmospheric conditions and engine/fuel properties. The prescribed ambient relative humidity over ice, RH ice,a , determines the ambient water vapour mixing ratio x v,a being proportional to RH ice,a · e s,ice (T a ). We use RH ice,a = 120% for the baseline case.

Initial conditions
Engine and fuel parameters like plume exit temperature, propulsion efficiency and water vapour mass emission index have 295 values according to an A340-300 aircraft at cruise speed (Table 1). This has been chosen consistently with the trajectory data setup. Only for the fuel consumption we prescribe a lower value than typical for this aircraft type in order to be consistent with the initial plume area given by the trajectory data. For given T 0 , combustion heat and propulsion efficiency, we obtainN 0 ≈ 75 according to Eq. (13). The initial plume water vapour mixing ratio is controlled by the water vapour emission index EI v

300
where the second term is the contribution of ambient humidity entering and exiting the engine. This term is less than 1% of the total plume water vapour and, therefore, neglected in our initialisation.
Based on laboratory measurements, we prescribe lognormally distributed soot particles with a geometric-mean dry core radius (r d ) of 15 nm and geometric width of 1.6 (Hagen et al., 1992;Petzold et al., 1999) and set the engine soot number emission index (EI s ) to 10 15 kg −1 as a typical value for commercial aircraft (e.g., Schumann et al., 2002Schumann et al., , 2013Kleine et al., 305 2018). Furthermore, we assume the same hygroscopicity parameter value (κ s = 0.005) as in Kärcher et al. (2015). The overall number of emitted engine soot particles per flight distance is given by shown that a 50:50 mixing of conventional kerosene with a biofuel causes a decrease in EI s by around 50% (Moore et al., 2017). Additionally, ground measurements with synthetic Fischer-Tropsch fuels have found a decrease of the average soot particle size besides the decreased soot particle number (Beyersdorf et al., 2014). This is because the switch to alternative fuels leads to a reduced number of primary soot spherules. The subsequent coagulation between those spherules decreases and the final soot aggregates are on average smaller. Therefore, in PV3, we vary both EI s andr d and combine a 50% reduction 315 in EI s with a reduction inr d to 12.5 nm. In PV4, we change the hygroscopicity parameter of soot particles in a reasonable range between 0.001 and 0.01 (pers. communication with M. Petters). Finally, we study the variation of contrail properties with ambient temperature T a (PV5) and relative humidity over ice RH ice,a (PV6).

320
Here, we introduce diagnostics that are used to evaluate the simulations.
The fraction of soot particles activated into water droplets (φ act ) or frozen to ice crystals (φ frz ) is given by where ζ i characterises, whether soot particles (ζ i,act = ζ i,frz = 0) represented by a SIP i have turned into water droplets (ζ i,act = 1 and ζ i,frz = 0) or ice crystals (ζ i,act = ζ i,frz = 1). Note that we classify ice crystals as activated, hence we get 0 ≤ φ frz ≤ φ act ≤ 1. 325 We introduce the apparent ice number emission index, AEI ice , as a measure of the number of formed ice crystals per fuel mass burnt. Since we consider only soot as exhaust particles in our study, AEI ice is simply given by the product of φ frz and EI s .
Furthermore, we define the effective radius of the soot/droplet/ice number distribution 330 Typically, r eff is a quantity that is relevant in radiation problems. Here, we prefer r eff over the ordinary mean radius since the time evolution is smoother and there is no need to specify the total particle number. ambient conditions fuel/engine properties engine exit conditions soot particle properties numerical parameters

Results
In this section, we first describe the spatio-temporal evolution of thermodynamic and contrail properties during the formation stage (Sect. 3.1). We then analyse the temporal evolution of contrail properties depending on different soot properties and 335 atmospheric parameters in Sect. 3.2. Finally, we summarise our results by studying the sensitivity of final contrail ice crystal numbers to ambient temperature and relevant soot properties. Our simulations are based on the FLUDILES dilution data and we particularly compare simulations using the full trajectory ensemble and a single average trajectory.

Spatial variation of contrail properties
The first two rows of Fig. 1 depict the spatio-temporal evolution of temperature and relative humidity. The panels on the right 340 show contour plots of RH liq in (r, t)-space. The size of the exhaust plume is immediately apparent from this representation.
The jet plume expands from an initial radius of 0.5 m towards around 7.5 m after t = 1 s. The black line indicates the plume radius r p as a function of time. All other panels show line plots, where the horizontal axis uses the normalised radial distancẽ r = r/r p as coordinate. Panel a) shows the radial temperature distribution for different plume ages, as indicated in the inserted legend. The black line shows the initial temperature profile with the smoothing in the outer areas according to the Crocco-Busemann profile. In the very beginning, the plume centre is not affected by entrainment and the core temperature remains at its initial value. The radial temperature gradient is initially large (nearly a temperature decrease of 300 K at t = 0.05 s) and decreases with increasing plume age as the plume centre cools. Later on, the absolute changes in temperature, spatially and temporally, decrease continuously.
Panel b) shows the radial RH liq distribution of a simulation with switched-off contrail microphysics, where water vapour 350 behaves like a passive tracer. We call the displayed quantity hypothetical RH liq ; it increases with increasing plume age and decreasing plume temperature due to the non-linearity of the saturation vapour pressure. Since plume cooling starts earlier and proceeds faster near the plume edge, water-saturation is surpassed earlier in this region (after 0.25 s) than in the plume centre (after 0.5 s). For fixed plume age, the hypothetical RH liq increases with increasing radial distance. This becomes also apparent from panel c), where the radial dependence of the hypothetical RH liq is depicted for several points in time.

355
Panels d), e) and f) show results of a simulation with switched-on contrail microphysics. As soon as droplet formation on soot particles is initiated (i.e., after plume relative humidity exceeds at least water-saturation), the further increase in RH liq is inhibited due to condensational loss of water vapour. The subsequent rapid decrease in RH liq (first at the plume edge and later on in the plume centre) is due to ice crystal formation and subsequent depositional growth. After around 1 s, the relative humidity with respect to ice reaches values of 100%, which translates into water-subsaturated plume conditions. Panel e) shows 360 the fraction of activated soot particles (φ act ). According to the evolution in plume relative humidity, soot particles activate into water droplets first near the plume edge and afterwards in the plume centre. Water saturation in the plume is first reached for T sat ≈ 240 K so that below those temperatures soot particles may form water droplets and grow further by condensation. The homogeneous freezing of the super-cooled water droplets typically occurs at T frz ≈ 230 − 232 K depending on the droplet size and the cooling rate. Panel f) shows the fraction of soot particles that turned into ice crystals (φ frz ). The radial profiles of φ frz 365 behave similar to the corresponding profiles of φ act but with a time lag of around 0.2 s. In the outer part withr > 0.5, basically all soot particles have become ice crystals after around 0.6 s. At t = 0.6 − 0.8 s, the profiles exhibit local minima atr ≈ 0.2, which results from the minimum in RH liq after around 0.7 s (Fig. 1 b).   Regarding the ensemble mean evolution, the fraction of soot particles forming droplets/ice crystals increases smoothly over a 380 certain time window. This is explained as follows: The spatial variability in RH liq is quite large (Fig. 1 d)) due to a systematic radial gradient and superposed turbulent deviations. Therefore, droplet and ice crystal formation occur for different trajectories at a different plume age mostly depending on the radial distance to the plume centre (Fig. 1 e) and f)). Although the ensemble mean features a smaller RH peak value than the average traj evolution, a larger fraction of soot particles forms ice crystals (96% vs. 88%).

385
In the near-threshold case (right column), T sat ≈ 236 K is lower and closer to T a . Hence, water-saturation is reached at a higher plume age when the cooling rates are already lower. Additionally, the peak RH values are clearly reduced. According to the evolution of the plume relative humidity, droplet and ice crystal formation on soot particles are initiated a few tenths of seconds later compared to T a = 220 K (panel d)). The time between the activation of the first droplets and the freezing of the last droplets is also longer (1.6 s versus 1.2 s). The final φ act and φ frz -values are clearly reduced relative to the baseline case since only the 390 larger soot particles can activate into water droplets due to the decreased peak RH liq . Again, the ensemble mean evolution displays slightly higher final φ frz -values than the average traj.
A closer look at the average traj evolution of the near threshold case reveals that φ act is not a monotonically increasing function and that the final φ frz -value is lower than the φ act -value at some intermediate time (in this example at t = 0.9 s). This is explained as follows: The largest soot particles are activated first and over time smaller soot particles grow large enough to be considered 395 activated as well. However soon after, RH liq drops below S K (for the smaller soot particles) as condensational growth removes water vapour. Those particles shrink, such that their radius drops below the activation radius. Once the largest water droplets freeze, depositional growth leads to a fast depletion of water vapour. Again, the smallest of the currently activated droplets face water-subsaturated conditions and also become de-activated (at around t = 1.25 s). We will encounter this "deactivation phenomenon" later again. Even though it is not visible in the ensemble mean evolution, this deactivation phenomenon also 400 occurs for most trajectories of the ensemble runs.

Variation with soot properties
In this section, we investigate the sensitivity of contrail formation to various soot properties. Fig. 3 presents results of the sensitivity studies PV1 to PV4. As in the previous section, the baseline and near-threshold temperature cases are considered.
We first analyse the baseline temperature case (T a = 220 K; left column in Fig. 3). The ensemble mean results reveal that 405 the temporal evolution of the activation fraction φ act and freezing fraction φ frz hardly changes with any variation of the soot properties. In the average traj framework, the sensitivity to soot properties is slightly higher. There, activation and freezing fraction moderately increase with decreasing EI s (panel a)) and increasingr d (panel b)). A higher number of soot particles leads to a faster depletion of plume relative humidity prohibiting the activation of smaller particles. Likewise, a decrease ofr d reduces the number of soot particles that manage to be activated. Changing the two parameters, EI s andr d , simultaneously 410 (panel c)), the effects on φ act and φ frz cancel out more or less resulting in an almost negligible effect. The variation of φ frz with the hygroscopicity of soot particles (panel d) displays the lowest impact within the considered parameter range. Note that the average traj simulations with higher prescribed EI s (dotted cyan lines) exhibit the aforementioned deactivation phenomenon.
In the near-threshold scenarion, the φ act and φ frz -values are clearly lower due to decreased plume water-supersaturation (Sect. 3.2.1). While the variation with EI s (panel b)) is still low, both average traj and the ensemble mean evolution are 415 more sensitive to the variation ofr d and κ than for T a = 220 K. This is because of two basic reasons: First, the critical saturation ratio for the activation of soot particles varies non-linearly with the soot core size (due to the non-linear Kelvin term) displaying a strong change for r d < 10 m (Fig. B1 a)). This means that far away from the formation threshold, larger changes of r d and κ are necessary to affect the activation of the several nanometer sized soot particles. Second, these small soot particles are at the tail of the size-distribution and, therefore, their activation makes only a low contribution to the overall φ act . In the 420 near-threshold case, only the larger soot particles (r d 20 nm) can form water droplets. This activation threshold is close tō r d which simultaneously defines the maximum of the lognormal size distribution. Hence, changingr d or κ has a significantly higher impact on φ act for T a = 225 K than for lower ambient temperature.
Again, we find a stronger sensitivity to r d (panel b)) than for κ (panel g)). Interestingly, the ensemble mean simulations show in the end larger r d and κ-induced variations in φ frz than the average traj simulations, which is opposite to what we saw in panel 425 g). This is mainly due to the deactivation phenomenon. In the combined variation of EI s andr d (panel f)), the impact ofr d on φ act and φ frz dominates over that of EI s . Halving the EI s -value and decreasingr d to 12.5 nm leads to a decrease in φ frz from 0.41 to 0.31 and translates into a decrease of AEI ice even by 62% (instead of around 50% if only the EI s -value is halved).
Furthermore, we have investigated the sensitivity of contrail formation to the geometric width of the soot particle size distribution σ (not shown): In contrast to the variation withr d and κ, there is a negligible impact near the contrail formation threshold 430 (T a = 225 K). This is because the activation of soot particles with dry core radii belowr d is already inhibited due to the Kelvin effect. For our baseline case, in general a higher fraction of soot particles forms ice crystals if the geometric width is reduced.
This is due to the fact that the smallest soot particles within the size distribution become larger and hence can be activated more easily. Decreasing (increasing) geometric width from 1.6 to 0.8 (2.4) causes a decrease (increase) in φ frz from 0.88 to 0.83 (0.98) for the average traj and from 0.96 to 0.93 (0.99) for the ensemble mean.

Sensitivity to ambient conditions
In this subsection, we study the dependence of contrail formation on ambient conditions. We vary the ambient temperature T a and background relative humidity over ice RH ice,a (see PV5 and PV6 in Table 2).
The top row of Fig. 4 displays the variation of ambient temperature T a at default RH ice,a = 120%. As already explained above, a change of T a affects the contrail formation time and the number of soot particles turning into ice crystals. The lower T a , 440 the earlier water saturation is first reached in the plume and the higher is the attained peak water-supersaturation. Hence, reducing T a to 215 K, droplet and ice crystal formation are initiated earlier and more ice crystals form (panel a)). The results for T a = 220 K and 225 K have already been described in Sect. 3.2.2. Thereby, the difference in φ frz is significantly larger between 225 K (dotted) and 220 K (solid) than that between 220 and 215 K (dashed).
The evolution of the effective radius r eff in panel b) shows a steep increase from the moment on when the first droplets are 445 activated. The increase slows down in most cases and is accelerated again, once the droplets freeze. The curves level off as soon as ice-saturation is reached. This happens after around 1 s for the lower T a and after around 2 s for T a = 225 K. The final r eff values hardly change between T a = 215 K and 220 K and between the average traj and the ensemble mean. Only for T a = 225 K higher values are attained. This is because the absolute amount of ambient water vapour is higher for higher T a and gets distributed among fewer ice crystals.

450
The bottom row of Fig. 4 shows the variation of RH ice,a . For T a = 220 K (panel c)), contrail ice nucleation is not at all affected by the variation of RH ice,a . In contrast, the sensitivity of φ frz to RH ice,a is very high for T a = 225 K. At ice-saturation (dotted line), only few ice crystals form (φ frz < 10%). This is because a decrease in RH ice,a causes a decrease in Θ G from 226.2 to 225.5 K so that those contrails form only 0.5 K below the SA-threshold temperature and maximum plume water-supersaturation is only a few percent. Analogously, more ice crystals form for RH ice,a = 140% as Θ G and hence |T a − Θ G | ≈ 3 K is higher.

455
Our results highlight that the sensitivity of contrail ice nucleation to atmospheric parameters is only large near the contrail formation threshold which is consistent with previous studies (e.g., Kärcher and Yu, 2009;Kärcher et al., 2015;Bier and Burkhardt, 2019). Therefore, it is meaningful to display the variation of the nucleated contrail ice crystal number with ∆T := T a − Θ G as well (see next section).

460
In the following, we aim at summarising the results from the previous sensitivity studies. We analyse the sensitivity of final contrail ice crystal numbers to ambient temperature and relevant soot properties. Thereby, we display both the freezing fraction (left column of Fig. 5) and the absolute ice crystal number in terms of AEI i,f (right column of Fig. 5). The latter is shown only for sensitivities with varying EI s since otherwise there is simply a linear relationship between φ frz,f and AEI i,f .
The freezing fraction in general decreases with increasing EI s (panel a)) because the plume humidity is depleted faster at 465 higher exhaust particle number concentrations and also distributed among more ice crystals. Hence, the smaller soot particles do not manage to form contrail ice crystals. The ensemble mean φ frz,f only slightly changes with EI s , which confirms the low sensitivity pointed out in Sect. 3.2.2. Accordingly, the ensemble mean AEI i,f increases nearly linearly with increasing EI s . The average traj φ frz,f -evolution (blue curve) features a stronger decrease than the ensemble mean, in particular for the higher soot number emissions (EI s > 10 15 kg −1 ). This is due to the deactivation phenomenon explained in Sect. 3.2.2. Therefore, the 470 increase in AEI i,f with rising EI s flattens in the average traj case.
Panel c) displays the variation of φ frz,f with ambient temperature T a (lower label) and the difference between T a and Θ G (∆T , upper label). Note that the relation between T a and ∆T is not purely linear since the ambient relative humidity over water varies with T a at fixed RH ice,a and hence, Θ G slightly changes as well. Near the contrail formation threshold (i.e., |∆T | 3 − 4 K), φ frz,f strongly increases with rising |∆T |. This is because the maximum plume water-supersaturation increases and more and more soot particles manage to activate into water droplets. Below T a = 220 K (|∆T | > 6 K), more than 90% of the soot particles turn into ice crystals so that AEI i,f approaches EI s . The average traj simulation produces somewhat lower φ frz,f values than the ensemble mean simulation. Furthermore, low soot simulations (with EI s reduced by 80% relative to the baseline case) feature a similar evolution in φ frz,f and a smaller difference between the average traj and the ensemble mean (not shown).
Finally, we investigate the contrail formation dependence on the soot particle size, as specified by the geometric-mean soot dry 480 core radiusr d . Panel d) shows a variation withr d for fixed baseline EI s . At T a = 220 K (solid lines), nearly all soot particles can form ice crystals as long asr d is larger than 20 nm. For smaller mean dry core radii, φ frz,f starts to decrease. This decrease is stronger for the average traj (blue vs. red curve). Near the formation threshold (dashed), φ frz,f values are clearly lower. Again, the ensemble mean values are larger than the average traj values. Moreover, we see a strong dependence onr d over the full parameter range.

485
In the third row of Fig. 5,r d is varied together with EI s such that the overall soot mass is fixed. We find a qualitatively similar evolution in φ frz,f compared to the previous sensitivity study (panel d)) so that the impact ofr d dominates. Yet, AEI i,f (panel f)) is controlled by the change in EI s (EI s decreases with increasingr d ) and the trend is opposite to that in panel d).
The following evaluation of our box model computations relies on the comparison with other contrail formation simulations.
Generally, such model simulations can differ in various aspects, which we group into four categories: the initial and ambi-490 ent conditions, the prescribed dilution, the model (micro)physics and the model framework (average traj, ensemble mean or fully 3D LES). If simulations have different settings or properties in several or all four categories, a comparison can hardly be conclusive as differences in the simulation results can have more than one origin. Hence, we made efforts to reproduce the simulation setup of the recent study by Lewellen (2020, abbr. as Lew20 hereinafter) by using his baseline conditions (more specifically those used in his Fig. 2). They differ from our baseline conditions as summarised in our Table 1    Lew20's dilution data are compared to the box model simulations of Lew20 (green vs. orange lines). This aims at validating the microphysical approach and implementation in both models. Results are presented for a monodisperse and bidisperse (as in Lew20) as well as for our regular lognormal soot size distribution. All simulations feature a pulse-like increase in AEI ice once freezing starts. However, ice crystal formation is initiated significantly earlier for Lew20's dilution, namely after 0.25 s instead of 0.6 s for the FLUDILES dilution. As explained in Sect. 2.3, the FLUDILES dilution is too weak and the onset of 510 contrail formation is retarded.
Using the bidisperse soot initialisation (blue dashed) we find AEI ice = 0.5EI s for the FLUDILES dilution as the 10 nm soot particles are not activated into water droplets. Yet, all soot particles are activated for Lew20's dilution. This is remarkable as for both dilution data the maximum plume water-supersaturation would be the same if the microphysical processes were suppressed. The reason for the different activation fraction is that the time period between droplet formation and freezing is 515 longer and correspondingly also the cooling rates are lower within this period in the FLUDILES case. This leads to a reduced peak saturation ratio, which is smaller than the critical saturation ratio for the activation of the 10 nm particles.
For the lognormal size distribution (blue solid), some smaller soot particles are only temporarily activated, but then shrink again as the larger droplets originating from the larger soot particles deplete water vapour. With Lew20's dilution (green solid), there is not enough time between activation and freezing to deactivate the smaller droplets. This difference originating from 520 the different dilution data implies that the deactivation phenomenon, which we frequently encountered in our results, may not occur that frequently if plume cooling proceeds faster. Further implications concerning this aspect will be discussed in Sect. 5.
As a next step, we compare our box model simulations using Lew20's dilution (green) with Lew20's box model simulation (orange). From the latter model only results for the monodisperse soot initialisation are available, yet the bidisperse case gives basically identical results. We find a perfect agreement between both models with an only some hundredth of seconds earlier 525 increase in AEI ice for the LCM. (Note that the green dashed and dotted lines are mostly covered by the solid line). This is because our freezing temperature (T frz ) is slightly higher than that of Lew20 (pers. communication with D. Lewellen) leading to an earlier onset of ice crystal formation.
The right panel displays FLUDILES ensemble mean simulations (red lines) and 3D LES results with online coupled microphysics of Lew20. The latter are shown for three different turbulence realisations demonstrating irreducible uncertainty 530 involved in predicting contrail formation. The ice crystal number from the 3D LES is not only considerably reduced relative to our AEI ice but also relative to the box model framework of Lew20 using the same microphysics as in his 3D LES. D. Lewellen explains this as follows: "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. . . . Thus, the much narrower time window for ice number production seen in the box 535 model leads to greater EI iceno [AEI ice ] than in the LES." Originally, we supposed that switching from the average traj to the ensemble mean framework leads to more realistic quantitative results as already assumed in Paoli et al. (2008). However, we see an opposite trend compared to the 3D LES, namely the ice crystal number rises when we switch from the average traj (blue lines in left panel) to the ensemble mean (red lines in right panel). We believe this comes from the fact that we do not consider the inter-trajectory mixing mentioned above. This means 540 that offline ensemble simulations miss an important process and their quantitative prediction of AEI ice is not more reliable than from an average traj simulation, even though they account for plume heterogeneity and are in several aspects more plausible.
Note that the reasoning above is based on comparing only a few simulations. For instance, Lew20 shows that for EI s ≤ 10 15 , all soot particles form ice crystals in the LES (see his Figs. 5 and 11) and hence there is no longer a difference between Lew20's average traj and LES predictions of AEI ice . Moreover, the deactivation phenomenon leading to the large difference between 545 FLUDILES average traj and ensemble mean results may be exaggerated with the FLUDILES dilution. This implies that the difference between the average traj and the ensemble mean could be smaller with a stronger dilution. We test this hypothesis in the following discussion section.
Similar to our study, Paoli et al. (2008) compared average traj and ensemble mean results for two specific parameter settings.
They also used a set of 25000 trajectories even though it is not the same data set provided by Vancassel et al. (2014). Several 550 aspects are similar to ours, e. g., the quasi pulse-like activation versus a formation period. Interestingly, they find that for the average traj the final ice crystal number is even three times lower than for the ensemble mean.
In the previous section, we have shown that ice crystal formation is initiated significantly earlier and the number of formed ice crystals may be even larger if we use Lew20's instead of the FLUDILES dilution data. Furthermore, we stressed that some 555 differences found between the ensemble mean and the average traj results may be particularly due to a too weak dilution. In the following, we will estimate the time lag of the plume mixing and discuss potential improvements in our results by "accelerated" trajectory data.
While 38% of the total combustion energy E tot is used for the propulsion of the aircraft (according to the A340-300 propulsion efficiency), the remaining energy (62%) is transferred into the jet plume in the form of kinetic energy E kin and ther-560 mal energy E therm . For the initial partitioning, we assume that E therm is in the range 51% -31% of E tot and, accordingly, refers to the smallest/largest assumed initial jet kinetic energy.
In the following, we introduce a small correction to "speed-up" our existing trajectory ensemble. More specifically: we artificially increase the dilution rate of each FLUDILES trajectory such that in the end the average dilution of this cor- simulations (blue and red lines, as already depicted in Fig. 5 a) and c)) and for the speed-up scenario (green and magenta) ac-575 cording to Lew20. We are particularly interested in comparing the evolution between the ensemble mean and the average traj.
In terms of the EI s -sensitivity, the difference between the two frameworks is apparently smaller in the speed-up sceanrio than in the original data set. Considering the variation with T a , the original average traj displays the lowest φ frz,f (blue) of all other simulations since the deactivation phenomenon leads to fewer ice crystals. For the speed-up scenario, the difference between the ensemble mean and the average traj is again significantly lower than for the original data set. Both sensitivity scenarios 580 confirm our hypothesis that the weak dilution of the original data is conducive to excessive deactivation in the average traj set-up. Furthermore, the ice crystal number is in general higher for the speed-up scenario (magenta) than for the original data set (red).
Additionally, we show the results of the parameterisation by Kärcher et al. (2015) in an adapted version. This parameterisation was recently implemented in a global climate model (Bier and Burkhardt, 2019), where it refines the contrail initialisation procedure. Those results (brown lines in Fig. 7) lie well in the range of the other simulations results. Overall, the extensive evaluation may help interpreting previous studies based on an average traj or a trajectory ensemble approach. (red) with the speed-up ensemble (magenta), whose average dilution matches that of Lewellen (2020). The blue and green lines show the corresponding average traj evolution of both ensembles. Additionally, the results from the parameterisation of Kärcher et al. (2015) with an adaption regarding the calculation of the activation dry core radius of soot particles (brown) are shown.

Conclusions and outlook
In this paper, the particle-based Lagrangian Cloud Module (LCM) (Sölch and Kärcher, 2010) has been extended to cover the specific microphysics of contrail formation. One of our objectives was to test the newly extended LCM model in a simple box 590 model framework prior to its coupling to the 3D LES model EULAG. Furthermore, our aim was to analyse the most relevant sensitivities of contrail ice crystal numbers to various soot properties and ambient conditions and to compare our results with previous studies.
We employ a data set with dilution histories along 25000 trajectories sampling the expanding jet. These data are provided by 3D FLUDILES simulations behind the engine nozzle of an A340-300 aircraft (Vancassel et al., 2014). We compare a single 595 box model run with an average dilution with an ensemble mean of box model runs for all dilution histories.
Plume mixing leads to a strong cooling at the plume edge where water-saturation is reached first. Later on, also the plume centre cools. Therefore, soot particles form water droplets (and subsequently freeze into ice crystals) first near the plume edge and afterwards in the plume centre.
Consistent with previous studies (e.g., Kärcher and Yu, 2009;Kärcher et al., 2015;Bier and Burkhardt, 2019), contrail ice 600 crystal number varies strongly with atmospheric parameters (like ambient temperature and relative humidity) only near the contrail formation threshold. Thereby, the difference between the ambient and Schmidt-Appleman threshold temperature is of relevance since it controls the maximum water-supersaturation in the plume.
Investigating the impact of soot properties on contrail evolution, we find that the fraction of soot particles freezing to ice crystals significantly changes with mean size and solubility of the soot particles only near the contrail formation threshold. This is 605 due to the combined effect of the non-linear Kelvin term and the lognormal soot particle size distribution. Absolute ice crystal numbers are, on the other hand, controlled by the soot number emission index for all ambient conditions. To investigate the impact of current biofuel blends, we analysed a combined variation of the soot number emission index and the average soot core radius. We find that the decrease in contrail ice crystal number does not only result from the reduced soot particle number but also from the smaller soot particle sizes.

610
We evaluate our study with results from a recent study of Lewellen (2020) who uses a similar contrail formation microphysics with a transient liquid phase. Prescribing the same initial and ambient conditions, we particularly find a later onset of contrail formation in our model. This is because our dilution is too weak as discussed in the previous section. When employing the dilution data provided by Lewellen (2020), we find a perfect agreement in the evolution of contrail ice crystal number suggesting correct implementations of the microphysics in both models.

615
Switching from a single trajectory to the ensemble mean framework shows only an improvement in terms of a more continuously (instead of pulse-like) evolution in contrail ice nucleation. On the other hand, we see an increase in final ice crystal numbers which is, at least for high soot number emissions, an opposite trend compared to the recent 3D LES results of Lewellen (2020). This is likely due to the fact that we do not consider the interaction (mixing and depletion of water vapour) between the different trajectories and, therefore, overestimate ice crystal formation. Hence, using a large spatially resolving trajectory 620 ensemble does not necessarily lead to improved scientific results contrary to what we expected in the beginning.
The presented aerosol and microphysics scheme describing contrail formation is of intermediate complexity and thus suited to be incorporated within 3D simulations (LES) of contrail formation explicitly simulating the jet expansion. Or next task is to extend the EULAG-LCM model system, which has been extensively used for high-resolution contrail simulations during the vortex and dispersion phase, by the specific contrail formation physics.
If droplet radii are initially small, i. e., having a magnitude comparable with the molecular mean free path, the condensational growth in the transition regime has to be considered. We use the correction terms for the mass and heat term in Eq.
(2) according to Fuchs and Sutugin (1971): 640 β t = Kn a + 1 (b 1 /α t,c + b 2 ) · Kn a + b 1 /α t,c · Kn a 2 , where b 1 = 4/3 and b 2 = 0.377 are empirical constants. The Knudsen number (Kn) describes the ratio between the mean free path of the surrounding gas molecules and the droplet radius so that Kn a = λ a r −1 .

645
The mean free paths of vapour (λ v ) and air molecules (λ a ) are given by Williams and Loyalka (1991) λ with R d and R v denoting the specific gas constants for dry air and vapour, respectively.
α m,c and α t,c are the mass and thermal accommodation coefficients for condensation. Here, we set both parameters to unity since 650 droplet growth rates from cloud model studies agree best with experimental studies when choosing these values (Laaksonen et al., 2005).
The transitional correction factors that appear in Eq. (7) are given by Pruppacher and Klett (1997) β where ρ a and c p are mass density and specific heat constant of dry air and α m,d and α t,d represent the mass and thermal accommodation coefficients for deposition, respectively. In this study, we set α m,d to 0.5 according to Kärcher (2003) To obtain the critical saturation ratio for activation of soot particles into water droplets (S c ), we search for the local maximum of 665 S K which requires ∂SK ∂r to be zero. The root r c is called the critical radius for activation. Substituting α = 1−κ and Ke = 2 σ Mw R T ρw yields ∂S K ∂r rc = −Ke r c 6 − 3 r d 3 (α − 1) r c 4 + r d 3 r c 3 Ke (α + 1) − α Ke r d 6 (r c 4 − r d 3 α r c ) 2 · exp Ke r c = 0.
The determination of r c requires the nominator of the first term f (r c ) = −Ke r c 6 − 3 r d 3 (α − 1) r c 4 + r d 3 Ke (α + 1) − αKe r d 6 (B3) 670 to be zero. We solve this numerically and apply a Newtonian iteration until the root r c is determined with an accuracy of < 0.01 nm. Figure B1 a) shows that S c increases with decreasing r d due to the Kelvin effect. The variation of S c with temperature T is quite low for the range where activation in exhaust plumes typically occurs (see legend). Using the simplified Eq. (10) of 675 Petters and Kreidenweis (2007) (green line) causes a significant overestimation of S c for r d < 15 nm. The critical saturation ratio increases with decreasing hygroscopicity parameter κ due to the solution effect (Fig. B1 b)).