Sensitivity of particle loss to the Kelvin effect in LES of young contrails

Different treatments of the Kelvin effect in LES modeling of early contrails are shown to cause variations in the survival rate of ice particles by up to a factor of 4 and in optical depth and mean particle size by up to 50%. The Kelvin effect which varies exponentially with particle size, can reduce or even suppress the impact of other important ambient parameters, such as ice supersaturation, on particle survival rate. Lowering or neglecting the Kelvin effect is shown to substantially alter the evolution of the ice particle size distribution and delay the onset of particle loss. A strongly Kelvin effect dependent exponential 5 relation between particle survival rate and particle size is shown for high EIsoot (O(1015)).

given as proof of Naiman et al. (2011)'s "outlier" behavior.This approach does not seem to provide sufficient information to determine whether an issue exists with the model.For example, what are the uncertainties in the parameter estimates in their Eqn.10(a-h)?What are the resulting z ∆ point-wise (or bin-wise) 95% confidence intervals to the curve fit (similar to the black dotted lines in our Fig.2(Left) in the revised manuscript, a global L-2 error is not sufficient to judge efficiency of the curve-fit over the entire range of z ∆ )?How do the models fare compared to the (presumably large) confidence intervals (77 data points to find minima of error and bias in a 6-D space)?
We have shown in our manuscript the major cause of discrepancy between models employing Lagrangian Particle Tracking(LPT) (Paoli and Shariff, 2016) without making claims of model errors or biases.

Comment:
Often such tests are added as appendices in more comprehensive studies.I agree that often details matter which require a somewhat longer description than possible in such appendices.Given the rather short length of the manuscript, I recommend to extend the study by a few aspects to make it a more substantial scientific contribution that is suitable to be published as independent research.
The presented results leave many questions open and I get sometimes the impression that the presented results are not fully understood.Some results seem implausible.At best, they are not explained well enough.Moreover, I think that the evaluation of some relationships or quantities is not overly useful.In general, the text touches many aspects only superficially.
I recommend publication, only after including a more careful analysis, a more thoughtful presentation and a better explanation of your results.

Response:
We intend for this manuscript to be concise and convey one major scientific finding -evaporative loss of particles in early persistent contrails is strongly sensitive to modeling of the Kelvin effect and this is the main source of discrepancy between some existing LES models.Many parametric models for contrails designed to be employed in larger GCMs use LES data from these discrepant LES models to infer parameter values (e.g.Unterstrasser (2016)).Some studies use LES data for young contrails to initialize simulations over longer time horizons (Unterstrasser et al., 2016;Jacobson et al., 2011).Such studies stand to directly benefit from our analysis.This is mentioned in the revised manuscript at p2.l16.We have included in the manuscript some of the most relevant explanations from the responses to the referee's comments given below so as to make the manuscript thorough in its analysis of the impact of the Kelvin effect.

Major Comments:
Comment: A1: Motivating your work with the low level of scientific understanding by citing a reference to Penner et al, 1999 (collecting scientific results about 20 years old) is awkward.Science has progressed since.Sausen et al, 2005, Lee et al, 2010, Burkhardt and Karcher, 2011, Bock and Burkhardt, 2016 Response: Newer studies have now been appropriately cited in the Introduction of the revised manuscript (at p1.l11-16).

Comment:
A2: The results section looks more like a technical report where hard facts are listed.Implications and interpretation of the simulation results are only touched in a few cases.I could live with it, if I understood all your model results.
However, this is not the case which is outlined in the following.

Response:
In the following responses, we have clarified the doubts raised by the referee.The implications of our results are concisely stated in the Conclusions section of the manuscript.Admittedly, we may have cut short on the discussion of our results in our attempt to keep the manuscript concise.Important clarifications sought by the referee are now included in the Results section of the revised manuscript.
A3: Ice crystal formation in contrails is completed after one second or so.Hence, the fraction of surviving particles as you show in Fig. 2 should be a monotonically decreasing function.It is not trustworthy when one simulation (red line) shows a late time increase after t=400s.So what's wrong, something in the model or in the post-processing tools?

Response:
Both the model and the post-processing tools are thoroughly vetted and we were unable to find any errors.Indeed, we observe an increase in the survival rate after t ≈ 400s for high EI soot and no Kelvin correction.Monotonicity is not an intrinsic property of the survival rate but a mere consequence of modeling assumptions.All particles reduced to the bare soot core (r = 15nm) are considered lost, but not deactivated for vapor deposition.In the absence of the Kelvin effect, crystal loss during the vortex phase is primarily due to plume being warmer than the ambient with additional aid from turbulent fluctuations (Kärcher et al., 2014).After the vortex phase, as the plume sloshes due to buoyancy and mixes with the ambient, its temperature relaxes to that of the local ambient (see Fig. 1(Left) in the revised manuscript) and RHi begins to rise.Now, these still activated soot cores can resume uptake of the available excess vapor.This is what is observed for the red curve in Fig. 3 (a, Top) of the revised manuscript.Cases with Kelvin effect are not seen to have this behavior as the lost particles never overcome the Kelvin barrier.Lewellen (2012) theoretically shows monotonicity of the particle survival rate only for cases with Kelvin effect.Schumann (1996) notes that for high temperatures as seen in the jet phase, most soot activation occurs through "condensation freezing mode" (Pruppacher and Klett, 2010).However, for low temperatures (as seen during the buoyant sloshing and dispersion phases), ice nuclei can be activated directly for vapor deposition (Pruppacher and Klett, 2010).Bailey andHallett (2004, 2002) observe that RHi < 125% is sufficient to activate particles for temperatures less than −42 o C.Even at the high temperatures observed during the jet phase, Kärcher et al. (2015) notes that in-situ measurements do not rule out direct vapor activation.In our model we assume that for low temperatures, soot is activated directly for vapor deposition (Paoli et al., 2004) while Picot et al. (2015) sets the higher bar of condensation freezing even during the dispersion phase.Noting that p water sat (T = 219K) ≈ 3.5P a and p ice sat (T = 219K) ≈ 2.1P a, for the RHi values considered in Picot et al. (2015), this means that a particle that is once reduced to its soot core is never reactivated -even in their no Kelvin cases.
It is pertinent to note that this difference in the model assumptions causes a mere ∼4% increase in the survival rate of our red curve in Fig. 3(a, Top) of the revised manuscript -the no Kelvin case that is anyway unphysical.We have noted this in the revised manuscript at p3.l21-24, p5.l20-26

Comment:
A4: Several aspects of the evolution of the size distribution (SD) in Fig. 3 look peculiar.
1.First of all, your chosen colours for B0 and B2 are hard to distinguish.Please replace one.

Response:
B2, B3 curves are now green while L0, L1 curves are now blue.The contrasting colors in each panel are now clear and easily distinguishable.
2.Why do you change the y-range for t ?60s, when there is no need to?

Response:
The purpose of changing the y-range is simply to make sure that the figure shows only the relevant range and reduce unnecessary white space.
3. The large droplet mode is controlled by RHi.Hence for the right end of the SD, all solid lines are basically identical and similarly all dashed lines.For t = 15, 30 and 45 s it looks like this may not be the case.Can you clarify this.

Response:
The plume RHi remains the same irrespective of the Kelvin effect modeling.So, the total uptake of available vapor for ice growth is the same, however its distribution between small and large particles is different based on how strong the Kelvin effect is.The near 100% plume RHi (Naiman et al., 2011;Paoli and Shariff, 2016) suggests that almost all available excess vapor in the plume volume is used up by the ice particles.We see the ice mass in the plume is indifferent to the Kelvin effect in Fig. 1(Right) of the revised manuscript.So, if a k = 0, some of the total excess vapor budget is accounted for by the abundant small particles leaving less for the large particles and hence reducing their prevalence, i.e. when Kelvin effect is lowered or suppressed, the vapor uptake by small particles results in the reduction in availability of vapor for, and thus, the prevalence of large particles, albeit marginally.This is what is seen in the SDs at 90s,120s and 300s (we presume these are the times the referee meant as the SD is almost identical for large sizes at 15s, 30s, and 45s).This has been clarified in the revised manuscript at p6.l16-21.
4. What really irritates me is the fact the left tail (the so called sublimation tail) develops faster for higher RHi.How can this be?There's three times more water vapour available for deposition in the RHi = 130%-cases.So why should ice crystals be more susceptible to sublimation in this case?This result is counter-intuitive, implausible and different to all other models.Hence, it must be explained in detail, why your models behaves like this.As a consequence, the ice crystal number in Fig. 2 drops faster for higher RHi.For the red case, even the final survival fraction is lower.This result is hard to believe.Did you carefully check all your model components?

Response:
We believe all model components have been implemented accurately.The slower "spectral ripening" of the size PDF for low RHi is explained below.
Let N tot (t) be the total number of ice crystals at time t.Until the onset of particle loss, the survival rate is near 100%.So, N tot (t) = N tot = const.. Let N (r, t) represent the number of ice particles of size ≤ r at time t.Thus, the size PDF is defined as n(r, t) = 1 Ntot ∂N (r,t) ∂r .
It follows that, where σ 2 (t) is the variance of the size PDF and a measure of its width.Until onset of crystal loss, the mean particle size grows very slowly (see Fig. 3(a, Center) of revised manuscript) so we have dr dt ≈ 0. Hence, we expect dσ 2 dt ≈ dr 2 dt .Using Naiman et al. ( 2011) (their Eqn.9), we thus have, The smallest particles are large enough to be indifferent to the Kelvin effect (99.98% of the particles are larger than 2×10 −7 m up to t ≤ 20s for RHi = 130% and t ≤ 35s for RHi = 110%).So, ∆ρ vapor excess is only a weak function of r and ∆ρ vapor excess is driven by the temperature experienced by the particles (up to the aforementioned times).This is sufficient to reason the faster spreading of the size PDF in high RHi cases -higher RHi ⇒ higher ∆ρ excess vapor ⇒ faster growth of σ 2 (t).Thus, the size PDF of cases with higher RHi starts having small particles (susceptible to Kelvin loss) earlier (∼ 20s for RHi = 130% and ∼ 35s for RHi = 110%).
In cases with Kelvin effect, once particles appear in size bins affected by the Kelvin effect (i.e.r bin ≤ 2×10 −7 ), the assumptions of the above analysis are no longer valid.Now, cases with higher a k loose particles faster due to apparent subsaturation experienced by them and consequently, the sublimation tail spreads faster for higher a k values.Simultaneously, we see a rise in mean particle size r = as N tot begins to fall while the total mass uptake is the same (given an RHi) for all a k values.Eventually the size PDF evolves in a self-similar manner (Lewellen, 2012) (their "slow growth regime").However, for a k = 0, the above analysis continues to hold until the smallest size bin (or the bin with lost particles, r bin = 15nm) gets populated (this happens earlier for RHi = 130%).Thus, case B4 sees early onset of particle loss than B5.As explained earlier, this "loss" in the case of no Kelvin correction is only temporary as after buoyant sloshing of the warm plume and mixing with the cool ambient, these activated soot cores regrow.
In Fig. AR 1 we have shown the evolution of σ 2 for high EI soot cases.We see that for a given RHi, the evolution of σ 2 is identical up to t = 20s for RHi = 130% and t = 35s for RHi = 110% and that σ 2 grows faster for higher RHi.Beyond these times, the cases with higher values of a k see faster growth of σ 2 .This discussion is now included in the revised manuscript in Appendix B. As to why other models do not observe this behavior, we will refrain from speculating as we are not privy to the specifics of their implementation.

Comment:
A5 -1: It is discouraging, if you mix up things and reviewers have to disentangle them: The terminology of your simulations is misleading.Runs B0, B1, L0 and L1 use a temperature-dependent a k which are compared to runs with constant a k = 1 × 10 −9 m.Your presentation implies that including the temperature dependence makes a large difference.However, the temperature dependence itself is not the reason for the observed differences (the temperature dependence is anyway weak, as you shows in Table 1 and Fig. 1).It is simply that the a k -value is about 2.3 times higher if you use your temperature-dependent expression.So you basically compare cases with a k = 1 × 10 −9 m and a k = 2.3 × 10 −9 m.

Response:
Using a k = 10 −9 m is an order-of-magnitude accounting of the Kelvin effect .For a k to be 10 −9 m as used in Picot et al. (2015) , the σ would have to be ∼ 0.047Jm −2 which is unphysically low for the temperatures relevant to development of persistent contrails (T = 218K, say) (Pruppacher and Klett, 2010).Using a physically plausible, temperature-dependent value of σ ≈ 0.107Jm −2 the a k value is obtained is ≈ 2.3 × 10 −9 m.Indeed, this means that we fix for B0, B1, L0, L1, U0 and U1, a k = 2.3×10 −9 m.The usage of the term "temperature-dependent" may be misleading and this has been removed in the revised manuscript.

Comment:
A5 -2: By the way, why do list a k -values for three different σ-values, if only one case is used

Response:
The purpose of listing a k values for 9 different (σ, T ) combinations, is to show that for a realistic range of σ and T , the value of a k is still ∼ 2 × 10 −9 .The dotted black lines in Fig. 2(Left) of the revised manuscript (that use extreme values of a k from our Table 1) then stress that a realistic variation in a k due to changes in (σ, T ) does not affect the Kelvin correction dramatically.However, an order-of-magnitude value of a k is seen to significantly change the correction factor. Comment: A5 -3: p4.l.9-10: "results in substantially lower Kelvin correction for smaller particles".This is misleading as the correction factor is constant over the whole radius range, only the absolute change is larger for smaller particles.

Response:
The Kelvin correction factor Φ = exp ( a k r ) is indeed a function of radius!So, the value of the parameter a k may be constant, but the correction factor experienced by different radii is certainly different -which is shown in Fig. 2(Left) of the revised manuscript.

Major to Minor Comment:
Comment: B1: The paragraph (p.3 l.18-l.20)sounds like a perfect motivation to carry out sensitivity simulations varying the initial size distribution.If the true initial size distribution is not known, a model offers the unique opportunity to vary this parameters.This is particularly interesting in this study.The Kelvin effect has a prominent effect on the shape of the size distribution as you show in Fig. 3.So a variation of the initial size distribution is directly relevant to the main aspect of your paper.This may be also a reason for discrepancies between models.

Response:
Thank you for this comment.The major discrepancies -low RHi sensitivity and correspondingly heightened EI soot sensitivity -are already adequately accounted for by the different Kelvin effect treatment (Please see p5.l9-12, p6.l28-34 in the revised manuscript).Besides marginally reducing the time to onset of particle loss, changing the initial size PDF is expected to not have any significant impact.In our above analysis, σ 2 (0) = 0, so evolution of a thicker initial size PDF (i.e.σ 2 (0) = σ 2 0 > 0) will put particles in r ≤ 2 × 10 −7 m range earlier to start particle loss earlier (also see Fig. 6 of Unterstrasser and Sölch (2010)).
This has been noted in Appendix B of the revised manuscript.

Comment:
B2: Similarly, I recommend to carry out the L4 and L5 simulations.You say, those simulations are not necessary, as Picot et al, 2015 showed that no crystal loss occurs.One main motivation of your work was the discrepancy between the various models.So in this sense, referring to another model is a bit contradictory.It would be interesting to know, if your model behaves similarly.

Response:
There is negligible discrepancy between our low Kelvin effect case for low EI soot and Picot et al. (2015).So, at least ex-post, Fig. 3(b) in the revised manuscript vindicate our choice of not choosing to run cases L4, L5.

Comment:
B3: To me it is unclear, what you want to demonstrate with the bottom row of Fig. 2. Mean particle size is mainly controlled by growth of detrained ice crystals being outside of the vortex system.The crystal loss, on the other hand, occurs inside the vortex system.For me it makes no sense to link those to quantities, as they are not really physically connected.I recommend to remove the paragraph from p.4 l.32 to p.5 l.5 and the sentence in the abstract/conclusions.

Response:
The argument about exponential relationship between size and survival rate is being made only for the vortex phase(our p5.l33 in the revised manuscript ...we observe an exponential relation between survival rate and size during the vortex phase).
As seen in Inamdar et al. (2016) (their Fig. 7) and also Unterstrasser (2014)(first panel of their Fig.1), there is hardly any secondary curtain when the coherent vortical structures are descending through stratified ambient -thus size and survival rate are certainly physically linked up to the end of the vortex phase and it is logical to seek a relation between them.This observation will help in modeling the size and survival rate in any model for contrail evolution up to the vortex phase.As we argue in our manuscript (p5.l35 of revised manuscript But, as the plume temperature...), this becomes progressively invalid after the coherent vortical structures are destroyed and the plume sloshes and mixes with the ambient while its temperature falls.In Fig. 3(a, Bottom) of the revised manuscript, we see that for a given a k , the plot for both high and low RHi is linear and has similar slopes up to the destruction of the coherent vortex system.After that, the curves for different RHi diverge and we may no longer seek an exponential relation beyond this point.In the revised manuscript we stress that this exponential relationship is valid only during the vortex phase at p6.l1.

Comment:
B4: Personally, I think that analyses of optical depth of young contrails are not overly useful, as this quantity is linked to radiation and climate aspects.LES of young contrails are not directly relevant to this.For this, contrails must be simulated over a much longer time (at least several hours).Optical depth decreases substantially over the first half an hour, as the contrail gets usually tilted by vertical wind shear (a process absent in your simulations).So the given optical depth values are only a snapshot.The differences you find may not be long-lasting.Indeed, Unterstrasser et al, 2016 presents contrail-cirrus simulations over 6 hours and switching off the Kelvin effect had barely an effect on contrail properties (all simulations were initialised with the same 5 min old contrail, though)

Response:
The referee points out that Unterstrasser et al. (2016) observe little sensitivity of contrail optical properties to the Kelvin effect past 5 mins., given a particular initialization (that was inferred from early contrail simulation which assumed one value of a k ).However, it is pertinent to ask -"How sensitive are these long time horizon simulations to different initializations?".So, if the same simulations were initialized using our B2 case vs. our B0 case, how different would the long-time predictions be?
Extinction is proportional to number of ice particles and we see by 5min.our B2 case to have ∼ 200% higher number of ice particles as compared to our B0 case.This has been noted at p5.l13-19 in the revised manuscript.Indeed, the OD values for young contrails are not directly useful for studying radiation and climate aspects.However the objective here is to see if we produce reasonable values of OD and extinction compared to other models employed in similar simulations for similar time horizons (Picot et al., 2015).

Comment:
B5: A point-to-point comparison between various models as done on p5 l21 is not leading anywhere.Contrail optical depth depends on a multitude of parameters.So you always find simulations with similar, yet not identical parameters which leave enough room for arguing that for this or that reason the optical depth is similar or smaller/larger in the one case.Unterstrasser, 2016 presents a more rigorous evaluation exercise that accounts for the multi-parameter nature of the problem and that is also able to disclose systematic model differences as mentioned in the introduction of this review.

Response:
As stated before, the objective here is merely to show that the model produces OD values that are "reasonable" and not "outliers".This is not uncommon in literature (Picot et al., 2015;Naiman et al., 2011;Paugam et al., 2010).We address the comment about Unterstrasser (2016) in our response to the referee's first comment.

Comment:
B6: Naiman speculated that they might have used too few computational particles and that this could have led to the discrepancy with other LES models.How many particles did you use?Do your results depend on this numerical parameter?Did you check convergence of your results?

Response:
We use a total of 8 × 10 6 simulation particles(SIPs)/numerical particles(NPs).Please see Appendix A of Naiman (2011) for a detailed analysis of sensitivity to this parameter.We see that the results are converged for this SIP resolution.

Minor Comment:
Comment: C1: I don't want to downplay the possible effect the early temperature difference by including/excluding the exhaust enthalpy has on contrail properties.Nevertheless, it is noteworthy that after 100s the excess plume temperature is not affected at all by this model aspect.

Response:
This is exactly our finding in Conclusion # 4.

Comment:
C2: I recommend to split Fig. 1 for clarity reasons.The left column shows LES results and the right column shows simple physical relations without a connection to your LES results.

Response
This has been done in the revised manuscript.

Comment:
C3: You cite several Inamdar papers from the recent past.I am not sure, if all those are peer-reviewed contributions.
If not, I recommend to reduce references to them and instead repeat the results in this study.For example, p.6, l.9/10 cites an important result of your recent work.Has it gone through peer review? Response: These are conference papers and are publicly available.Results from these papers are being compiled into Inamdar et al. (in progress, 2016b, i).Citing conference papers is not uncommon in the field, e.g Huebsch and Lewellen (June 5 to 8, 2006, San Francisco, CA, USA) is referenced widely in contrail LES literature.

Comment:
C4: Can you add the expression for σ? Do you vary it independently of T? The legend of Fig. 1 right alludes to this.

Response:
As stated in the manuscript, the value of σ is held constant for "temperature dependent" cases at 0.107Jm −2 .

Comment:
C5: p.2 l.27:What do you want to say here?Can you make a clearer connection between the availability of measurement data and what you say in the rest of the sentence.

Response:
Febvre et al. ( 2009) reports measurements for ∼ 2.5min.old contrails, Gayet et al. (2012) reports measurements for 70s, 105s, 205s old contrails.We wished to see if the impact of exhaust enthalpy subsides by ∼ 100s or else the exclusion of it may have affected our ability to compare simulation data with observations made at those time horizons.As it turns out, exhaust enthalpy will not affect this comparison.This has been elaborated in the revised manuscript at p3.l5-9.

Comment:
C6: p.3 l.33:The plume temperature is constant!?I do not understand this statement.The plume temperature increases due to adiabatic heating.It may help if you describe in more detail how you compute the excess temperature.
How is your reference temperature determined? Response: