Contrail radiative dependence on ice particle number concentration

Recent studies on low aromatic fuels have shown that lower soot number emissions may reduce contrail ice particle number concentrations (N ice). Here we implemented, in a sophisticated radiative transfer model, two ice particle size distribution schemes in order to estimate the contrail radiative forcing’s (RFs) dependence on these prospective N ice reductions resulting from the introduction of sustainable aviation fuels. The results show that an 85% contrail N ice reduction produces a 35% smaller contrail RF, while neglecting all non-radiative effects. This estimate of an RF reduction only considers the effects of the N ice change assumed here, and neglects other potentially important microphysical mechanisms that may change the relationship between soot number emissions and N ice. A comparison of our results with previous published estimates from full climate model simulations, shows similar RF reductions to those which also take into account non-radiative mechanisms, evidencing the need for more studies in order to allocate the contribution from radiative and non-radiative changes, as this would guide possible mitigation implementations. Despite these modeled contrail RF reductions being largely independent of the assumed ice water content (IWC), it is only through simultaneous improvement of the IWC and N ice representation in models that contrail RF estimates can be better constrained. This is because our calculated RF varied by a factor of 3 when assuming a ±30% IWC range; and by a factor of 5 if a, still conservative, ±60% IWC range was prescribed, suggesting that the differences in the prescribed IWC and N ice values in different models may explain the large discrepancies amongst published RF estimates. Recent estimates of higher N ice values, and lower IWCs found in contrails even after several hours, compared to surrounding cirrus under similar atmospheric conditions, were assessed to conclude that it is mainly the differences in IWC that make young contrails have a smaller RF, and to reduce our previous estimate for linear contrail RF for year 2006 by 65%.


Introduction
The ice particle number concentration (N ice ) of contrails depends upon the initial soot aerosol number concentration emissions [1,2] for kerosene-based fuel. Given that, in turn, contrail radiative forcing (RF) is dependent on N ice , understanding the dependency of contrails' RF on N ice becomes relevant in the ongoing discussion around the usage of lower-carbon footprint fuels (often referred to as 'Sustainable Aviation Fuels' , or 'SAF'). SAF have low to zero aromatic content, and many measurements at ground and at altitude show that the soot number emissions depend upon this (Moore et al. [3] and references therein). Understanding how SAF might impact the optical properties and otherwise of contrails [4,5] is of critical importance to better-quantifying the potential climate co-benefits of low-carbon footprint fuels in reducing contrails' RF.
Burkhardt et al [6] estimated, from global climate model (GCM) simulations, that an 80% reduction in N ice resulted in a 50% reduction in contrail cirrus RF, and a 50% reduction would lead to a significant decrease in contrail cirrus optical depth and coverage, leading to a decrease in RF by approximately 20%. Similarly, Bock and Burkhardt [7] found a 14% RF reduction from an N ice reduction of 50%. Both studies account for changes in contrails coverage, lifetime and optical depth, while the aim of the present study is to estimate the sensitivity of the radiative-only properties of contrails with respect to differences in N ice in order to isolate their dependence on the prospective introduction of SAF, assuming a relationship between soot number and N ice.
Contrails are ice clouds formed by aircraft under sufficiently cold and humid atmospheric conditions [8]. They may persist for several hours in ice supersaturated regions, eventually losing their linear shape, and becoming contrail cirrus as they increase their horizontal and vertical extent [9]. The confidence in the understanding of the RF of linear contrails remains constrained by the difficulties in the measurement, both in situ and from remote sensing, of their complex microphysical and optical properties, as well as by the complexity involved in parameterizing their broad variability [10] in large-scale models.
The present study is not only aimed at assessing the differences in the radiative properties of contrails linked to the introduction of SAF, but also at estimating the radiative differences between natural and contrail cirrus, as well as the sensitivity of their calculated RF to differences in the N ice and ice water content (IWC) amongst different models.
Most ice radiative contrail and cirrus parameterizations make use of two variables, commonly an effective size, like the equivalent diameter (D e ), and IWC. D e is defined as the ratio of the surface area to the volume of the ice particle size distribution (PSD), while the IWC is simply the mass concentration for a given volume. The fact that the smallest ice crystals largely determine N ice without significantly influencing D e , nor IWC, implies that with this approach two PSDs may show similar D e and IWC values despite large N ice differences, which may significantly influence the calculated radiative properties. This has already been shown [11][12][13] for an atmospheric single column, as well as for ice single scattering properties, but not for global coverage RF calculations, prompting the need for the present study. Equally important is the use of a sophisticated efficient radiative transfer model that allows a realistic representation of the particle size dependence on IWC and temperature. This has not always been the case in previous contrail sensitivity studies (e.g. Markowicz and Witek [14]) in which the contrails' optical depth, or altitude, are prescribed using fixed values. This approach neglects the dependence that contrail IWC and D e have on temperature and therefore on altitude. Furthermore, these simplifications neglect the fact that the results may be different if a different fixed value is assumed.
Given the lack of climatological contrail data, some linear contrail studies have prescribed natural cirrus properties (e.g. Ponater et al [15] and Rodríguez De León et al [16]), concluding that, like in the case of natural cirrus, the calculated RF of persistent contrails is mainly driven by their IWC, while being marginally dependent upon particle size, but evidence has grown about contrails showing higher N ice than natural cirrus, as reported in several studies [17][18][19][20], suggesting the need to reassess this conclusion.
Despite the currently limited information about contrail N ice and IWC properties and their temporal evolution, young contrails have been shown to have higher N ice values than natural cirrus by an order of magnitude [19,20], with enhanced N ice values being found even after several hours compared to surrounding cirrus under similar atmospheric conditions [18]. A more recent study [21] found more modest N ice differences of around 50% between natural and contrail cirrus. Based on these findings here we also here estimate the calculated contrail net RF sensitivity to N ice and IWC differences.

Model setup
The model setup is described in the supplementary information section, and is similar to that of Rodríguez De León et al [16] (RL18, hereafter), from which the reference case is here included for comparison purposes. The calculated linear contrail cover is prescribed as monthly mean 6-hourly vertical profiles, based on the REACT4C [22] air traffic inventory for year 2006. Given that the vertical layers are variable and predetermined in the model, their IWC was reduced in order to approximate a thickness of approximately 500 m. Also, given that only one ice cloud type can be prescribed in the model, we run the simulations under clear sky conditions.
In order to assess the sensitivity to N ice , two radiative parameterizations were here used to estimate the global linear contrail instantaneous RF. These parameterizations are based on two PSD schemes with different N ice values. The PSDs are defined as functions of ambient temperature and IWC.
The first scheme was developed by Field et al [23] (F07 hereafter) and the second one, which was previously implemented in RL18, was developed by McFarquhar and Heymsfield [24] (MH97 hereafter). The MH97 scheme is based on PSD measurements of deep convection cirrus from the Central Equatorial Pacific Experiment (CEPEX), while the F07 parameterization is based on mid-latitude ice clouds. The MH97's PSD parameterization provides N ice and its partition into the small and large ice particle PSD modes for any specified IWC, which is here prescribed, together with its dependence on temperature, from the IWC climatology developed by Schiller et al [25] (S08, hereafter). The contrail radiative properties are then calculated by integrating their single scattering properties over the MH97 PSDs in order to produce a parameterization in terms of the PSD's IWC and D e . For the F07 scheme, the radiative parameterization described in Baran et al [26] was applied, which is based on an ensemble of ice crystal shapes, described in Vidot et al [27]. This radiative parameterization does not include the dependence on the effective ice particle size, because the particle shape ensemble implicitly incorporates the size dependence from the combined IWC and ambient temperature values.
It was shown in RL18 that N ice values from the MH97 scheme were slightly larger than the climatological values reported in Krämer et al [28]. In contrast, the F07 scheme values are smaller than those from the MH97 parameterization by around 85% at representative persistent contrail temperatures, and for this reason these schemes provide a wide range of modeled N ice values to analyze their impact on the calculated RFs.
For the MH97 PSDs, we assumed a hexagonal cylinder particle habit with a variable aspect ratio, which reproduced the observed mass and area contributions from small and large particles reported by MH97. The D e values prescribed in the simulations are also obtained from the MH97 scheme, ensuring that a realistic dependence of D e and of the IWC on the local ambient temperature was always assumed in the sensitivity cases. As mentioned in the introduction section, the accuracy of the results from sensitivity studies relies on prescribing physically realistic values, which take into account their interdependencies.
The IWC was prescribed as a function of local ambient temperature in all the cases, but in order to assess the RF dependence on the IWC, we scaled it by 30% and 60% increments/decrements with respect to the mean IWC compiled by S08 for middle latitudes, assumed in the base case. Prescribing ±60% IWC differences results in an IWC range covering a factor of 4. This can be contrasted to the ice water path differences between GCMs used in the IPCC's 4th Assessment Report, which showed a factor of 20 difference between the largest and smallest values, and a factor of 6 after removing outliers [29], while reporting regional differences of nearly 2 orders of magnitude. Given that the ice water path would also include an uncertainty linked to the physical thickness of the cloud layer, it seems reasonable to assume that the ±60% IWC uncertainty implemented here is similar to that between GCMs. Also, contrail cirrus are reported to have a 60% smaller IWC compared to cirrus under the same atmospheric conditions [21].

Results
The RF results, shown in table 1, comprise longwave (LW), shortwave (SW) and net (LW + SW) annual mean global instantaneous contrail RF values, as well as their percentage differences with respect to the MH97's base case. The cases combine scaled IWC values with the two described PSD schemes, amongst which N ice values differ by around 85% at representative persistent contrail temperatures, emulating plausible differences between contrails and natural cirrus [30], as well as possible soot reductions linked to a prospective introduction of SAF. Figure 1 shows the percentage differences with respect to the MH97 base case, from which it can be seen that the N ice differences resulted in the F07 net RFs being around 35% smaller than those calculated with the MH97 scheme. This reduction resulted from F07 SW RF magnitudes being around 33% larger than those from MH97, partially counterbalanced by a comparatively modest strengthening of around 11% in the LW. As a result, all F07 net RFs showed smaller values than the MH97 base case, even for case 5 which had a 60% larger IWC prescribed value.
The combined N ice and IWC sensitivity in table 1 covers a range of net RF values that spans a factor of 3 when prescribing IWC differences of ±30%, and a factor of 5 when IWC is varied by ±60%, compared to a more modest factor of 3 for the ±60% IWC scaling if only one PSD scheme is prescribed. This shows that despite the fact that the results are more sensitive to the prescribed IWC differences compared to those from N ice when taken separately, their contribution to the combined RF sensitivity is not largely dissimilar, and therefore it is only by simultaneously improving their model representation and validation that the uncertainty of linear contrail and contrail cirrus RF estimates can be constrained.

Discussion
The impact that differences in N ice play on the radiative properties of contrails acquires particular relevance when considering theoretical calculations that have shown a dependence of N ice on soot number emissions [2], as well as limited in situ measurements [31] reporting soot and N ice reductions between 50 and 70%. This discussion becomes especially pertinent to current efforts to produce low-carbon footprint fuels (bio, synthetic), since they will have low to zero aromatic content, and much lower soot emissions [3].
Contrasting the 50% reduction in contrail cirrus RF reported by Burkhardt et al [6] from an 80% N ice reduction, which took into account changes both in contrail cirrus optical depth and coverage, with the 35%  [25] assuming relatively large and small N ice values from McFarquhar and Heymsfield [24] and Field et al [23] schemes, respectively. RF reduction found in the present study for a similar reduction (85%) in N ice , we conclude that, despite the fact that our model only considers the change in radiative properties linked to N ice , both results show an encouraging similarity. It becomes clear that, given the relevance of the topic, more studies from different models and approaches will be needed in order to understand the implications of introducing SAF. It should be stressed here that although a reduction in RF is found with lower N ice , that in the context of low-aromatic SAF such a response is not assumed. Theoretical modeling has shown that at low soot number emissions, enhanced activation of ultrafine aqueous particles may occur at very cold temperatures, resulting in an increase in N ice for low soot numbers, past a critical threshold of 10 14 soot particles kg −1 fuel.
The effect of N ice differences had already been explored by Zhang et al [11] for single column simulations for spherical ice crystals hexagonal columns and polycristals, but not for global simulations nor for ice particle shape ensembles. The results presented here showed that the effect of an 85% N ice reduction is a reduction of around 35% in the calculated net RF, irrespective of the assumed IWC.
IWC differences between contrails and natural cirrus have been also measured. Li et al [21] reported that countrail cirrus present smaller IWCs by around 60%, while Chauvigné et al [19] reported smaller contrail and contrail cirrus IWC values in all stages of contrails compared to natural cirrus, and larger N ice for young contrails suggesting the MH97 scheme better represents linear contrail N ice values. Based on these results, we reduced our previous estimate from 9 mW m −2 , reported in RL18 (represented here as case 3 for MH97), to 4 mW m −2 , based on case 1 of the MH97 scheme, which assumes an IWC reduction of 60% as reported by Li et al [21]. This significant reduction evidences the importance that contrail measurement campaigns have on reducing global RF uncertainties by improving combined IWC and N ice climatological data compared to surrounding natural cirrus.
Important discrepancies amongst published studies on the linear contrail net RF's dependence on the prescribed particle size include RL18 reporting a RF sensitivity of −9% when implementing the MH97 parameterization used here, while Yi et al [32] reported differences of +46%. In both studies the differences emulated a 50% increment in the effective particle size, suggesting that, as pointed out by Zhang et al [11] and Rodríguez De León and Haigh [12], the commonly adopted approach of parameterizing the radiative properties of ice crystals using their IWC and an effective size produces RFs which are not univocally estimated unless some information about the PSD's shape is also included.
Here we suggest that information about the N ice values of the PSDs used to parameterize the radiative properties of persistent contrails can provide useful information to explain differences amongst RF estimates. Similarly, the sensitivity to particle size obtained by Yi et al [32] would also be affected by the fact that in their simulations they assumed fixed contrail altitude and optical depth values. As a result of this, the RF differences would be biased in one direction compared to a setup in which ambient conditions would more realistically determine contrail properties. It is therefore suggested, that despite the advantages of implementing simplified setup cases in sensitivity and inter-comparison studies, these simplifications often produce results not only linked to the assessed variable.
The radiative differences from the PSD schemes used here are also related to the fact that they assume two different particle shapes, hexagonal cylinders and an ice crystal ensemble. Similarly to the present study, Yi et al [32] also compared schemes based on a hexagonal cylinder habit and on an ice habit mixture. They obtained a 11% larger net RF from hexagonal cylinders, which they attributed to the fact that roughened ice particles were assumed in their multi-habit mixture, which would present smaller SW asymmetry factors compared to smooth hexagonal cylinders, and therefore would reflect more SW radiation.
This, in turn, would result in a reduction of the net RF for the multiple habit mixture. The discrepancy between the 11% difference obtained by Yi et al [32] and the 35% obtained in the present study may be explained by the use of fixed optical properties and contrail altitude assumed by Yi et al [32], which would exclude the N ice variability expected at different temperatures if the altitude and the optical depth were allowed to vary. Again, these results suggest that simplified setups may render their sensitivity estimates unrealistic when fixing properties.
The results in table 1 show −SW/LW values of consistently around 67%, and 83% for the MH97 and the F07 scheme, respectively, showing that, with all other variables being fixed, this ratio shows little variability for each PSD scheme. On the other hand, this ratio is expected to be linked to the contrails' altitude, as De León et al [33] reported −SW/LW percentages between 64% to 81% when fixing contrail altitudes between 30 000 and 38 000 feet, while allowing IWC and the effective size to vary. Schumann and Graf [34] compiled −SW/LW values from various published studies, and reported a similar range size, from 32% to 52%, not only linked to N ice differences, but in that case comprising differences in the contrail fields, parameterizations, radiative transfer schemes, and the background meteorology differences amongst models.
The results in table 1 and figure 1 show that the F07 scheme, which comprises smaller N ice values by 85% compared to the MH97 scheme, produces net RF values that are always smaller than the reference case; even cancelling out the effect of an IWC 60% increase. In contrast, by excluding the IWC variability, the net RF differences linked to the N ice differences produce a net RF reduction of around 35%.
A model intercomparison study by Myhre et al [35], reported contrail RF differences amongst models within a factor of 2 when prescribing a homogeneous 1% global contrail cover at a fixed altitude. They also reported global net RF values between 9.3 and 15 mW m −2 when prescribing a realistic 3D linear contrail cover. Their contrail setup consisted of a fixed PSD with fixed optical properties, which is equivalent to fixing N ice , meaning that their setup excluded the N ice sensitivity analyzed here, instead quantifying RF differences linked to the radiative schemes and background meteorologies amongst the compared models. Their 3D results implied a relatively modest (22%) standard uncertainty compared to 60% differences from the combined N ice and IWC obtained in the present study if we assume a 30% IWC uncertainty, and 80% differences if a 60% IWC uncertainty is assumed, suggesting that the discrepancies between published contrail RF estimates are more likely explained by the uncertainties in the ice properties fed into the radiative schemes than by other differences in the radiative treatment in large scale models, including the background meteorology.

Conclusions
The relevance of the dependence of the contrails' RF on their assumed N ice is manifold. The results obtained here showed that smaller N ice values produce a SW RF enhancement that is 3 times that of the LW, resulting in a NET RF reduction. Given that our calculations only cover the radiative aspects of an N ice reduction, our results of a contrail net RF reduction potentially linked to an N ice reduction by introducing SAF can be expected from the corresponding change in optical properties alone, while a similar RF reduction has been reported from GCM calculations in the literature. This is an encouraging result, but until further measurements are made, a simple relationship between lower soot numbers from SAF and N ice cannot be assumed since theoretical modeling shows that at very low soot numbers (<10 14 soot kg −1 fuel) and cold conditions (⩾12 K below contrail formation threshold temperatures), ultrafine aqueous particles may act to increase N ice again [2]. Nonetheless, reducing soot, and therefore N ice , at temperatures above contrail formation values minus 12 K, to this threshold point of 10 14 soot kg −1 fuel, could represent a mitigation option for contrails by reducing N ice for fossil-based fuel, with no added potential CO 2 penalties, as opposed to other strategies like contrail avoidance by flight rerouting.
Based on the differences in N ice and IWC between natural cirrus and contrails reported in the literature, our results show that a 60% IWC reduction, similar to the differences found by Li et al [21], produces a larger net RF reduction (∼55%), compared to that (∼35%) from the reported larger N ice values also reported for contrails. It is nevertheless concluded that the accuracy and representation of these variables combined needs to be improved in order to better estimate contrail RF and its dependence on the introduction of SAF.

Data availability statement
The data cannot be made publicly available upon publication because they contain commercially sensitive information. The data that support the findings of this study are available upon reasonable request from the authors. cover from ERA-Interim was used by COMA to calculate the critical relative humidity for cloud formation. The specific humidity and temperature were then used to determine the potential contrail coverage. This coverage was then multiplied by the travelled distance in each cell to produce the actual contrail coverage, assuming an aircraft fleet propulsion efficiency of 0.36. This resulted in an annual global contrail cover of 0.10%, based on traffic for the year 2006, which incorporated a representative diurnal variability.
The instantaneous RF, defined as the difference at the top of the atmosphere of the radiative fluxes with and without contrails, was calculated using the UK's Met Office's GCM's two-stream radiative transfer scheme (SOCRATES) [38] with a 144 (longitude) × 72 (latitude) × 23 (altitude) resolution. Solar zenith angles and day lengths at the middle of the month were used to represent the monthly mean SW fluxes with a five points Gaussian integration. The contrail radiative properties were represented offline in SOCRATES as functions of two prognostic variables, IWC and ambient temperature. The latter is prescribed using ERA40 [39] ECMWF reanalysis data, from which all background meteorology is prescribed as monthly climatological averages for water vapor, temperature, and surface albedo.
The parameterizations used in the simulations are schematically represented in table A2. The contrail's IWC is determined using its climatological dependence on ambient temperature developed by Schiller et al [25] and used as input in McFarquhar and Heymsfield's [24] (MH97 hereafter) PSD scheme by choosing six contrail representative temperature values, covering a range from 195 to 245 K. Three fits from Schiller et al [25] for minimum, mean and maximum IWCs for each temperature are then used to produce 18 PSDs. These PSDs are then used to integrate the ice crystal single scattering properties (extinction coefficient, single scattering albedo, and asymmetry parameter) of randomly oriented hexagonal cylinders to parameterize these properties in terms of the PSD's IWC and generalized diameter (D ge ) [40] which ranged from 11 to 75 µm.
The methodology presented in Fu et al [40] is used to integrate the LW optical properties calculated by Baran et al [41]. In the SW, the parameterization described by Fu [42] and the single scattering properties calculated by Yang et al [43] are applied, integrated over 6 and 9 bands in the SW and in the LW spectral regions, respectively. The aspect ratio (ratio of the ice crystal length and width) of the prescribed hexagonal cylinders varies with particle size in a way that matches both the total volume and projected areas retrieved by MH97, making the micro-physical and the optical properties in the model consistent.