Contrails and Their Impact on Shortwave Radiation and Photovoltaic Power Production-A Regional Model Study

A high resolution regional-scale numerical model was extended by a parameterization that allows for both the generation and the life cycle of contrails and contrail cirrus to be calculated. The life cycle of contrails and contrail cirrus is described by a two-moment cloud microphysical scheme that was extended by a separate contrail ice class for a better representation of the high concentration of small ice crystals that occur in contrails. The basic input data set contains the spatially and temporally highly resolved flight trajectories over Central Europe derived from real time data. The parameterization pro5 vides aircraft-dependent source terms for contrail ice mass and number. A case study was performed to investigate the influence of contrails and contrail cirrus on the shortwave radiative fluxes at the earth’s surface. Accounting for contrails produced by aircraft enabled the model to simulate high clouds that were otherwise missing on this day. The effect of these extra clouds was to reduce the incoming shortwave radiation at the surface as well as the production of photovoltaic power by up to 10 %.


Abstract.
A high-resolution regional-scale numerical model was extended by a parameterization that allows for both the generation and the life cycle of contrails and contrail cirrus to be calculated. The life cycle of contrails and contrail cirrus is described by a two-moment cloud microphysical scheme that was extended by a separate contrail ice class for a better representation of the high concentration of small ice crystals that occur in contrails. The basic input data set contains the spatially and temporally highly resolved flight trajectories over Central Europe derived from real-time data. The parameterization provides aircraft-dependent source terms for contrail ice mass and number. A case study was performed to investigate the influence of contrails and contrail cirrus on the shortwave radiative fluxes at the earth's surface. Accounting for contrails produced by aircraft enabled the model to simulate high clouds that were otherwise missing on this day. The effect of these extra clouds was to reduce the incoming shortwave radiation at the surface as well as the production of photovoltaic power by up to 10 %.

Introduction
Contrails consist of ice crystals formed in the exhaust plume of aircraft due to mixing of the hot and humid exhaust with cold environmental air. Contrails form in case the Schmidt-Appleman criterion (Schmidt, 1941;Appleman, 1953;Schumann, 1996) is fulfilled; i.e., the ambient temperature is be-low a threshold of around −45 • C. With plume temperatures near −38 to −40 • C, contrail particles, which are formed in the liquid phase initially, freeze homogeneously and quickly to form ice crystals. Those ice crystals grow in air with relative humidity above ice saturation and contrails persist. With such a favorable state of the atmosphere, the originally lineshaped contrails undergo various physical processes at the microscale, spread by the influence of shear and sedimentation, and change their structure and microphysical properties. At some point, contrails are no longer distinguishable from natural cirrus in observations. This type of anthropogenic cloud is then called contrail cirrus (Heymsfield et al., 2010) and can have lifetimes of several hours. In a particular case, an 18 h old contrail cirrus could be tracked in a satellite imagery (Minnis et al., 1998). Other examples of longlifetime contrail observations are summarized by .
Contrails influence the radiative budget of the atmosphere in a way that is comparable to that of thin natural cirrus clouds (Sausen et al., 2005). Although a number of observations (e.g., Iwabuchi et al., 2012;Vázquez-Navarro et al., 2015; and other modeling studies have been performed, important properties, such as the optical depth or the spatial and temporal extent of occurrence, have not been sufficiently investigated and are also not quantified to a satisfactory extent (Boucher et al., 2013). The radiative forcing of aged contrails and contrail cirrus is of greater importance than that originating from young, line-shaped contrails (Stubenrauch and Schumann, 2005;Eleftheratos et al., 2007;Burkhardt and Kärcher, 2011).
Previous model studies primarily used GCMs (global circulation models) extended by parameterizations that are able to simulate line-shaped contrails. Here, the global radiative forcing due to contrails was quantified. The values range from roughly 5 to 20 mW m −2 depending on the simulated year and assumed contrail properties (Marquart et al., 2003;Stuber and Forster, 2007;Kärcher et al., 2010) and are consistent with satellite-based estimates by Spangenberg et al. (2013). Taking the larger effect of contrail cirrus into account, a global mean radiative forcing of approximately 38 mW m −2 (Burkhardt and Kärcher, 2011) (recently updated to 56 mW m −2 by Bock and Burkhardt, 2016b) to 40 to 80 mW m −2 (Schumann and Graf, 2013) is found. The major drawbacks of these methods are the coarse resolution of the models in both space and time.
Another class of studies simulates single contrails with high-resolution LES (large-eddy simulation) or RANS (Reynolds-averaged Navier-Stokes) models, hereby focusing on contrail formation (Paoli et al., 2013;Khou et al., 2015), young contrails with ages up to 5 min and their interaction with the descending wake vortices (Lewellen and Lewellen, 2001;Unterstrasser, 2014), and the transition into contrail cirrus over timescales of hours (Unterstrasser and Gierens, 2010a;Lewellen, 2014). Usually, the contrail evolution is studied for a variety of environmental scenarios, which helps to single out important process and ambient and aircraft parameters. Parameter studies allow for investigating the conditions under which contrails are persistent and the manner in which microphysical and optical properties change during transition and decay. Because this method is not applicable for a larger number of contrails and often limited to idealized environmental scenarios, it is not suited to quantify the impact of air traffic on the state of the atmosphere.
The spatial scale applied in this study lies between those typically used for large-eddy simulations and global climate models. We use the online coupled regional-scale model system COSMO-ART (Baldauf et al., 2011;Vogel et al., 2009). In this context, online coupled means that meteorology, chemistry, and contrail-related processes are simulated in one model at the same grid and one main time step for integration is used. For this study, the model is extended by introducing a new hydrometeor class for contrail ice crystals. The contrail initialization is based on a parameterization of early contrail properties by Unterstrasser (2016). The prognostic equations for contrail ice are similar to those of the natural cirrus ice class. Despite their similar treatment in terms of model equations, the evolution of the two cloud types can be quite different, as their formation mechanisms are different. Especially with respect to the spatial scale used in the regional-scale COSMO-ART model, the presented study is complementary to the aforementioned approaches and is, to our knowledge, the first study of its kind.
Ice crystals in young contrails are considerably smaller than the ones that form in natural cirrus clouds and occur with substantially higher ice crystal number concentrations (Febvre et al., 2009;Schröder et al., 2000;Voigt et al., 2010). Therefore, the original microphysical scheme was extended by a new hydrometeor class that allows for a separate treatment of the small contrail ice crystals separate from natural ice. This approach allows for the investigation of contrail microphysical properties and their changes during the various stages of development represented in a regional atmospheric model. In contrast to other studies using a two-moment microphysical scheme (Bock and Burkhardt, 2016a), the presented model configuration uses no fractional cloud coverage and no prognostic equations for contrail geometric properties like volume or area which describe the contrail spreading on the subgrid scale. Compared to GCM parameterizations, this omission seems acceptable in a regional model, as the spatial resolution is much higher (horizontal grid size of 2.8 km versus 50 km in a GCM).
Regarding the contrail microphysics and the interaction with the meteorological situation, the entire procedure is online coupled, thus allowing feedback processes between contrails and natural clouds in contrast to other models on a comparable grid scale (Schumann, 2012). There, a mixed Lagrangian-Eulerian approach is used instead of the usual Eulerian treatment. This approach allows covering the scale ranging from thousands of single contrails to multiyear global climate simulations (Schumann, 2012;Schumann et al., 2015;Caiazzo et al., 2017).
One of the key goals of this study is thus to quantify the influence of contrails and contrail cirrus on natural high-level cloudiness. Moreover, the radiative properties of contrails and their local influence on the shortwave (SW) radiative fluxes at the surface are examined.
The model uses a diagnostic radiation scheme (Ritter and Geleyn, 1992). Because the description of SW optical properties for ice clouds is optimized for various naturally occurring crystal habits (Fu et al., 1998;Key et al., 2002), a separate treatment of contrail ice crystals is introduced here as well.
Another feature of this study is the new and recently compiled data set of flight trajectories. Rather than statistical calculations for globally averaged fuel consumption, the basic data consist of real commercial aircraft waypoint data (flightradar24.com, 2015). Another approach using commercial flight data to study contrails on a regional scale is described in Duda et al. (2004). Here, a combination of commercial flight data and coincident meteorological satellite remote sensing data was used to perform a case study of a widespread contrail cluster. In the future, further case studies should be performed for which in situ observations of natural ice clouds and especially contrails are available.
The presented model configuration serves to study microphysical evolution of contrails and contrail cirrus, their influence on natural high-level cloudiness, and their impact on the radiative fluxes on a regional scale and short time periods. This gains importance, e.g., in predicting the energy yield from photovoltaic (PV) systems. During the past decades, the development of alternative, clean energy production was enhanced to counteract global warming and reduce air pollution. Within the scope of methods, one of the most promising sources is solar energy, gained by PV cells. To assure a sustainable supply, the demand of precise prediction of the energy yield from PV systems is desired (Lew and Richard, 2010). Several approaches exist to forecast PV power, such as statistical models, neural networks, remote sensing models, and numerical weather prediction models (Inman et al., 2013). Especially PV forecast using numerical weather prediction models is challenged by special weather situations or phenomena that are poorly represented in the models (Köhler et al., 2017). For example, Rieger et al. (2017) found a large impact of mineral dust due to Saharan dust outbreaks on the solar radiation over Germany. This becomes important, as mineral dust is currently not considered adequately in operational weather forecast. Another phenomenon that is not represented at all in numerical weather prediction is the influence of aviation.
In Sect. 2, the modification of the microphysical scheme and the radiation scheme are described. Section 3 presents a description of the parameterization providing the source terms for contrail ice. In Sect. 4, the results of a case study and comparison with satellite observations are presented.

Model description
In this section, the parameterizations to calculate the microphysical properties of ice crystals and the modifications to represent contrails are presented.
The model system COSMO-ART comprises a detailed treatment of aerosol dynamics and gas-phase chemistry (Vogel et al., 2009). Although most of these features are not used for the simulations shown in this study, large parts of the infrastructure contained in COSMO-ART, e. g. the tracer structure and modules for reading emission data, are adopted for the contrails parameterizations.
In this study, the COSMO-ART model is coupled with a comprehensive two-moment microphysical scheme following Seifert and Beheng (2006). Until now, the scheme contained one cloud ice class (besides other hydrometeor types in warm and mixed-phase clouds) that describes ice crystals in high-level ice clouds. Because ice crystals in freshly formed contrails are considerably smaller than those in natural cirrus, the basic microphysical processes in young contrails are treated in a separate, newly introduced contrail ice class. From now on, the original ice cloud class is called cirrus ice class. The separate treatment allows simulating local bimodal size spectra.
Consequently, also in the radiation scheme, contrails are treated separately from the cirrus ice using the new con-trail ice class. The applied radiation scheme is described in Sect. 2.2.
Unless otherwise indicated, the parameterized processes for the new contrail ice class are the same as those used for the cirrus ice class.
Similar approaches with a separate contrail ice class using climate models with coarser grid size are described by Burkhardt and Kärcher (2009) and Bock and Burkhardt (2016a).

The contrail ice class
In this section, we first describe shortly the treatment of ice crystals in the cirrus ice class. Hereby, we mainly focus on those aspects that are relevant for understanding contrailspecific modifications in the contrail ice class explained later on. A longer description including various formulae as well as a table showing the different coefficients characterizing both the contrail and the cirrus ice class (Table A1) is deferred to the Appendix.
In each grid box, the ice mass distribution is assumed to follow a generalized distribution f (m) (see Eq. A1). Prognostic equations are solved for the ice crystal number concentration n (zeroth moment of f (m)) and the ice crystal mass concentration q i (first moment of f (m)). Written concisely, the prognostic equations for the moments M k of order k = 0 or 1 (see Eq. A2 for a definition) read as where u is the grid-scale mean wind, K h is the turbulent diffusion coefficient, and v sed,k is the number or mass-weighted sedimentation velocity (see Eq. A7). S k comprises all source and sink terms like nucleation, deposition/sublimation, and aggregation. The deposition source term is derived from the growth equation of a single ice crystal, which is integrated over the whole ice crystal mass spectrum.
The ice crystals have a hexagonal shape and the mass m of a single ice crystal is related to its size L via the mass-size relation For m in kilograms and L in meters, the parameter values are a geo,nat = 1.59 and b geo,nat = 2.56 (Axel Seifert, personal communication, 1 June 2017) and are valid in the size range [L min = 17.5 µm, L max = 3800 µm]. The treatment of contrail ice is analogous to that of natural cirrus with only a few modifications. Nucleation is switched off in the contrail ice class. Instead, the generation of contrail ice depends on air traffic and the atmospheric state and is explained in Sect. 3.1.
As mentioned before, most ice crystals in contrails are smaller than in natural cirrus and different values for a geo and b geo (values are given in the Appendix) assure reasonable aspect ratios, also for small ice crystals down to sizes of L min = 1.24 µm. The upper size limit is set to a relatively small value of L max = 58 µm. If the mean ice crystal size in a grid box exceeds L max , then the total ice crystal mass and number from such a grid box are transferred from the contrail ice class to the cirrus ice class. This is reasonable, as contrails show distinct bimodal size spectra with many small ice crystals with sizes around 10 µm and fewer large ice crystals in the fall streaks (Unterstrasser et al., 2017a;Lewellen et al., 2014). The contrail ice class contains predominantly small ice crystals and the cirrus ice class allows for larger ice crystals that may also stem from aged contrails or contrail fall streaks. One drawback of this approach is that the anthropogenic contribution in the cirrus ice class cannot be directly determined. Instead, we analyze the differences in the cirrus ice class between simulations with and without air traffic. This indirect quantification of the aged contrail contribution could be circumvented by introducing further contrail ice classes and may be implemented in the future.
The sedimentation parameterization and other components are as in the cirrus ice class.
The introduction of a second ice cloud class leads to a more complex behavior as both ice cloud classes are coupled and interact with each other in several ways. They are directly coupled via the collection process (see Appendix) and the mass transfer of large contrail ice crystals as described above. Moreover, they interact indirectly via the competition for the available water vapor and possibly via dynamical changes through diabatic processes.
For a first illustration of our approach, Fig. 1 shows average size distributions (ASDs) of a simulation with and without air traffic. Each ASD is a superposition of local distributions. In this example, the contrails are at most 4 h old. More details on the simulation setup follow in Sect. 4.1. The solid blue line shows the contrail ice class, which has a peak at sizes around 20 µm. A less pronounced maximum is located at about 2 µm. The contribution of aged contrails becomes apparent by comparing the two dashed lines. Those show the cirrus ice class ASDs of a simulation with air traffic (blue) and without air traffic (green). Apparently, the anthropogenic contribution is substantial, in particular in terms of total number. The ice crystals in aged contrails are on average smaller than in natural cirrus and the peak size of the SD is shifted to a smaller value.

The radiation scheme
The atmospheric radiative fluxes in the COSMO-ART model are calculated using the GRAALS (General Radiative Algorithm Adapted to Linear-type Solutions) radiation scheme (Ritter and Geleyn, 1992). The algorithm needs as input several quantities such as temperature, pressure thickness of the model layers, trace gas concentrations, as well as the cloud cover and the mass mixing ratio for each hydrometeor class considered (Ritter and Geleyn, 1992). To include contrails and contrail cirrus in the radiative algorithm, we include a contrail ice cloud cover determined from the contrail ice class mass mixing ratio. Grid cells where the ice mass mixing ratio exceeds 10 −8 kg kg −1 are considered to be covered with contrails or contrail cirrus. The same threshold value is used in Seifert and Beheng (2006) for grid-scale natural ice clouds. As mentioned before, the aviation contribution to the natural cirrus ice class can only be determined by comparison with the reference simulation. Within the radiative algorithm, optical properties of hydrometeors are calculated. These are the mass specific extinction coefficient, single scattering albedo, asymmetry factor, and the forward peak of the phase function. For ice clouds, the parameterizations following Fu et al. (1998) and Key et al. (2002) are used. Because they are optimized for natural ice clouds, the scheme computes reliable values for effective radii r e between 5 and 60 µm (Fu et al., 1998), whereas for ice crystal populations with radii smaller than 5 µm the parameterization is not well defined. Ice crystals in young contrails often have effective radii smaller than 5 µm. To overcome this problem, we are using the parameterization of Fu et al. (1998) and Key et al. (2002), but we prescribe a lower limit of 5 µm for the calculation of the optical properties. For q i , no limit is prescribed; instead, the simulated q i is used to calculate the optical properties of the ice crystals. As the radiation scheme uses q i and r e for determining the radiative fluxes, implicitly fewer but larger crystals are assumed here. The extinction due to small ice crystals is expected to be larger than that for larger ones. Therefore, in our study, the radiative effect of young contrails may be underestimated.
The calculation of the contrail effective radii follows Fu et al. (1998). Here, the ice crystals are assumed to be randomly oriented, hexagonal columns.
D denotes the half width of an ice crystal, which is implicitly given by the mass size relation (see Eq. B2). The size distribution f L is related to the mass distribution f (m) via the transformation property f L (L)dL = f (m(L))dm. One can show that f L (L) (number concentration per size range) follows a generalized distribution, if f (m) (number concentration per mass range) does so and the mass size relation is a power law (see Eqs. B3 and B4). Other parameterizations exist that can compute reliable values for optical properties of small ice crystals with sizes down to 0.2 µm (Bi and Yang, 2017). For future studies, using such an approach clearly could overcome the necessity of the threshold described above.

Formation of contrails
For the description of contrails, the first step is to check whether the environmental conditions are favorable for the formation of contrails. Here, the Schmidt-Appleman criterion (Schumann, 1996) is used, which defines a threshold temperature below which contrail formation occurs. In the second step, the source term of contrail ice mass and number has to be calculated. The parameterization used to calculate those source terms is described in detail in Unterstrasser (2016).

Initial values for contrails
Microphysical properties of aged contrails depend a lot more on the number of ice crystals than on the ice mass after the vortex phase (Unterstrasser and Gierens, 2010b). The initial ice mass is of minor importance, as the later growth of contrail ice crystals and the related ice mass evolution in a persistent contrail is mainly controlled by the ambient water vapor supply. In contrast, the ice crystal number changes only slowly in a long-living contrail. Hence, its initial value determines the typical ice crystal sizes in the evolving contrail cirrus (for a given environmentally controlled ice mass), which affects the radiative properties and the sedimentationrelated life cycle. Therefore, it is appropriate to explicitly prescribe an initial ice crystal number concentration instead of an initial ice mass. The presented procedure is applied at each model time step and for each grid cell, given that both an aircraft is present and the Schmidt-Appleman criterion is simultaneously fulfilled. The parameterization provides ice crystal numbers for contrails that are about 5 min old. As meteorological input parameters, it requires the temperature T at cruise altitude, the ambient relative humidity with respect to ice RH i , and the Brunt-Väisälä frequency N BV . Furthermore, aircraft properties are characterized by the water vapor emission I 0 , an "emission" index for ice crystals EI iceno , and the wing span b span . We determine I 0 for medium fuel flow at cruise conditions as assumed in Unterstrasser and Görsch (2014). Here, a simple parabolic fit for I 0 depending on the wing span b span is used (Unterstrasser, 2016). Information on the wing span is available in the flight track data (see Sect. 3.2).
In this study only the most common JET-A fuel is assumed to be used; therefore EI iceno is set to 2.8 × 10 14 kg −1 following Unterstrasser (2014). The total number of ice crystals formed in the beginning, N 0 , is calculated using the following equation: with water vapor emission index EI H 2 0 = 1.25. Note that EI iceno is not reduced when the ambient temperature is only slightly below threshold temperature, even though in such situations fewer ice crystals would form . The descending movement of the primary wake of an aircraft causes adiabatic heating within the plume. Due to this, sublimation and loss of ice crystals occur, even in a supersaturated environment (e.g., Unterstrasser, 2016). As the spatial and temporal resolution of the model is too coarse to simulate these processes, the fraction of ice crystals surviving the vortex phase is parameterized by a loss factor λ N s . Details can be found in Unterstrasser (2016). The total number of surviving ice crystals per flight path N s is calculated with For the initial ice mass produced, the water vapor emission I 0 is used. The values for the produced ice mass and ice crystal number per flight distance are distributed equally on all grid cells, the aircraft passes within a time step: Here, V cell denotes the volume of the grid cells, and d is the flight distance within a grid box.
In global models, contrail parameterizations usually contain further prognostic equations that describe in some way the bulk contrail geometry (e.g., fractional cloud coverage or even some measure of contrail depth). In our approach, we assume that contrail ice crystals always populate the whole grid box as the horizontal scale is much smaller than in GCMs. In the vertical direction, it is assumed that ice crystals are distributed over the whole grid layer. Close to the ground, the vertical grid size is about 10 m and increases to 300 m at the tropopause. In supersaturated conditions, contrail depth varies between 100 and 500 m (mostly depending on aircraft type and stratification) and is similar to the depth of the grid layer. In the horizontal plane, the simplifications could have a larger effect. If few flight routes transect a grid box and the segments are in total d AT = 10 km long, then this implicitly results in an hypothetical initial contrail width of 2.8×2.8 10 km ≈ 800 m. This is larger than what is typically observed for 5 min old contrails and better fits to 15 min contrails (Freudenthaler et al., 1995). Hence, disregarding fractional coverage smears out the initialized contrails to some extent.

Determination of flight tracks
In contrast to most of the previously mentioned global modeling studies, this study uses a new and recently derived data set. Rather than statistical calculations for globally averaged fuel consumption, or radar data, the basic data consist of traffic waypoint information over a limited area recorded from ADS-B (Automatic Dependent Surveillance -Broadcast.) transponders on the plane (flightradar24.com, 2015). The ADS-B data are obtained mainly from flightradar24.com (2015). The DLR holds a historic data file that can be purchased from flightradar24.com (2015). These data are cleaned and combined with the input of the Official Airline Guide Database that is also held by DLR.
From this information, a data set is compiled containing spatially and temporally resolved information on geographical position, height, current velocity, and type of aircrafts. For reasons of efficiency, the information on geographical position is interpolated to fit onto the grid of the model. As a proxy for the aircraft emission parameters (see Sect. 3.1), the wing span b span is used.
The data set contains 8 h of air traffic, beginning at 08:00 UTC on 3 December 2013 and ending at 16:00 UTC on 3 December 2013. The trajectories of all flights during this period are displayed in Fig. 2.

Case study
To test the methods described in the previous sections, a case study was performed. For this purpose, a situation over Germany with a high density of contrails in an otherwise cloudfree environment was chosen as the simulation period.

Model setup
On 3 December 2013, the meteorological conditions over Central Europe were favorable for the formation of contrails. Additionally, the natural high-cloud coverage was relatively low, thus allowing for the identification of contrails on satellite images.
For the case study, two simulations were conducted, both running for 24 h, starting on 3 December 2013, 00:00 UTC. They use a horizontal 2.8 × 2.8 km grid and 60 vertical levels, resulting in a mean distance between the model layers of 300 m in the upper troposphere and 15 m in the lowest model layer. For boundary data, hourly COSMO reanalyses were used. In the reference run, air traffic is turned off and the cirrus ice class is the only ice cloud class to be active. For the run with air traffic, which we call aviation simulation, the previously explained configuration with two ice cloud classes is employed. Air traffic is switched on at 08:00 UTC and the two simulations evolve identically up to this point. Practically, a spin-up phase shorter than 8 h could have been used, but from an operational point of view it was simpler to start the simulations at 00:00 UTC. As the data set of flight trajectories contains no information about air traffic after 16:00 UTC, no new contrails form after this time.

Simulated contrail properties
The contrail treatment in the microphysics scheme is designed such that contrail-induced changes occur both in the newly introduced contrail ice class and in the existing cirrus ice class. The contribution of young contrails can be directly assessed by evaluating the contrail ice class. The contribution of aged contrails is found by comparing the cirrus ice class of the aviation simulation and of the reference simulation.     ice class of the aviation simulation, whereas the middle and right columns show the cirrus ice class of the aviation and the reference simulation. At this time, contrails are at most 2 h old and mostly consist of numerous very small ice crystals.
Independent of the considered quantity, the contrail ice class features mostly line-shaped structures. The IWC reaches values up to 5 mg m −3 , which is larger than in the simulated natural cirrus. This hints at an accumulation of emitted water vapor additional to the ambient supersaturation. Furthermore, the absolute humidity at the heights considered is relatively low. Therefore, the IWC of natural cirrus is small and cirrus clouds are very thin and almost invisible. The ice crystal number concentrations often lie between 1 and 100 cm −3 and can exceed those in natural cirrus by a factor of up to 1000. Consequently, the ice crystal effective radii in contrails are much smaller and lie below 10 µm. Figure 4e shows that the most dense contrails over Northern Germany already leave a mark in the cirrus ice class (see black box in panel b). Notably, each line-shaped structure in the left panel corresponds to a pair of lines in the middle  panel. This indicates growth of ice crystals particularly at the margins where the contrail ice mass is soon transferred to the cirrus ice class. Consistently, LES studies (Unterstrasser and Gierens, 2010a;Lewellen et al., 2014) and in situ measurements (Petzold et al., 1997;Heymsfield et al., 1998) indicate the strongest growth at the edges of a contrail. A closer inspection (not shown) reveals that each line-shaped structure in the aforementioned box consists of several contrails. Those were produced by several aircraft that fly along the same route with short time separations.
Nevertheless, the IWC, n, and r e values of the aged contrails and the surrounding natural cirrus shown in Fig. 3 are similar. Figure 4 show the situation at 12:00 UTC, analogous to Fig. 3 for 10:00 UTC. As time progresses, existing contrails continue to grow and further contrails are generated. Four hours after air traffic was activated in the model, a major part of the model layer is filled with contrails (left panel). Lineshaped patterns are still identifiable in some places, particularly over Northern Germany. South of 52 • N many contrails overlap and represent a huge contrail cluster. Based on their visual appearance, no clear distinction from natural cirrus would be possible.
Those contrail clusters still feature high IWC, high n, and low r e values comparable to those 2 h earlier and distinct to the surrounding natural cirrus. Moreover, we find lots of aged contrails over Germany (black box in middle panel), a region that is basically cirrus-free if there is no air traffic (black box in right panel). Here, local enhancements in both IWC and n occur in the cirrus ice class, where aged contrails are transferred to (Fig. 4e), compared to the reference simulation (Fig. 4f).
The results indicate that, in observations, microphysical criteria may help to separate at least young contrails from natural cirrus. In general, contrail fall streaks and aged contrails cannot be identified as such once their linear shape is lost or masked. Unterstrasser et al. (2017b) shows that natural cirrus that forms in high-updraft scenarios can have ice crystal numbers similar to those of young contrails, which renders even the separation between young contrails and natural cirrus impossible. Moreover, they show that contrails that be-come embedded in natural cirrus have large volumes where ice crystals of both origins co-exist. Hence, it is no longer meaningful to try to draw a strict separation line between natural cirrus and the anthropogenic cloud contribution.
Next, we analyze vertical distributions along the black line depicted in Fig. 4. Figure 5 displays the relative humidity (RH i ) and reveals a remarkably thick layer with strong supersaturation (maximum value: 1.34) that extends over the complete east-west extent of the model domain. The layer depth increases from 2.5 km in the west to more than 4 km in the east. These are generally very favorable conditions for the persistence and spreading of contrails.
The fidelity of such high values of RH i is corroborated by a comparison with observations. Here, vertical profiles obtained from radio soundings (UWYO, 2018) and simulated data evaluated at radiosonde stations are depicted in Fig. 6. In general, high values of RH i are observed by radio soundings, even supersaturation occurs over Lindenberg and Meiningen. The model clearly overestimates RH i for Idar-Oberstein, but Meiningen and Lindenberg agree quite well.
Additionally, the black curve shows the mean vertical profile of the cross section displayed in Fig. 5. From this it becomes apparent that the relative humidity in the displayed cross section is remarkably high compared to the radiosonde locations. The large layer of supersaturation is caused by lifting and radiative cooling. It can persist, as natural cirrus clouds are located mostly below this layer (Fig. 7b). Also, the cirrus clouds present in the layer are too thin, i.e., occurring number concentrations are too low (Fig. 7e) to effectively reduce supersaturation. Figure 7 shows the same ice cloud properties as Fig. 4 and additionally the extinction coefficient E at a wavelength of 1.115 µm. Again, the left column shows the contrail ice class of the aviation simulation, whereas the middle and right columns show the cirrus ice class of the aviation and the reference simulation.
In the reference simulation, natural ice is present over the entire supersaturated area. The number concentrations are mostly small, leading to optically thin cirrus clouds with extinction coefficients hardly exceeding 10 −2 to 10 −1 km −1 . Rather large ice crystals are present throughout the supersaturated area, with the largest values of r e found in the lower part of the cirrus around 10 to 14 • E.
In the aviation simulation, it becomes obvious from the left column plots that contrails form only at altitudes between 11 and 13 km. This is caused by the absence of air traffic below this layer. As already seen in Fig. 4a and b, the IWC is comparable to that of natural cirrus, whereas number concentrations reach much higher values. The effective radii in the aviation simulation are typically 1 order of magnitude smaller than in natural cirrus clouds and do not exceed 10 µm. These results are in agreement with large-eddy studies (Unterstrasser and Gierens, 2010a; Lewellen et al., 2014) and in situ observations (Poellot et al., 1999). The numerous small crystals lead to high values for the extinction coefficient (up to 2 km −1 ). Therefore, contrail ice is of great importance for the radiation budget.
In Fig. 8, the relative occurrence of ice crystal number concentrations and temperature for the cross section shown in Fig. 7 is depicted. The relative occurrence is normalized with the sum over all values. Both reference simulation (Fig. 8a) and aviation simulation (Fig. 8b) are similar for higher temperatures (i.e., lower heights) up to 220 K. For lower temperatures, high number concentrations up to 7 cm −3 occur in the aviation simulation, whereas number concentrations clearly decrease strongly with temperature in the reference simulation. Here, a rough comparison to measurement data can be made. In Voigt et al. (2017), mid-latitude cirrus clouds and contrails were probed in situ during an aircraft measurement campaign. Comparing their Fig. 6b to Fig. 8, a similar increase in n below temperatures of about 220 K is found. Therefore, most likely, the high values occurring in the aviation simulation and not in the reference simulation are due to aviation induced clouds.
In the aviation simulation, changes in the natural ice class can be found mainly at heights where contrails form and slightly below. The enhancement of ice number concentrations occurs mostly at flight levels, whereas an increase in IWC is found below. During the initialization, contrail ice crystals are vertically distributed over the whole grid layer and this indirectly accounts for the initial wake vortex induced vertical expansion of a contrail. Within our simulation, the vertical structure of the contrails is determined only due to the gravitational settling of the larger ice crystals. During this process, ice crystals number concentrations tend to decrease. In contrast, only a slight increase in r e is found. In areas where contrail ice enters the cirrus ice class, a large increase in extinction coefficient occurs. Values of the extinction coefficient are comparable to those of the contrail ice class and can be as high as 1 km −1 .
After 16:00 UTC, no new contrails are initialized in our simulation. Four hours after the end of new contrail formation, the remaining contrail ice has been advected to the southern part of domain (Fig. 9). Local patterns of increased number concentrations in the cirrus ice class are now limited to those regions, where contrail ice is still present. The line-shaped structures in the contrail ice class vanish, but relatively small values for r e are still found throughout the domain.
In the northern part of domain, the cirrus ice class is again undisturbed by aviation. Here, no remarkable differences to the reference simulation can be seen.

Comparison with satellite observation
In the following, satellite images (created with Global Imagery Browse Services (GIBS) NASA/GSFC/ESDIS, 2018) are shown in Fig. 10 for a qualitative assessment of the simulations. The panels a and b show the MODIS Terra Corrected Reflectance True Color at 10:00 and 12:00 UTC for the sim- ulated day, respectively, both with a resolution of 250 m. The True Color composition consists of MODIS bands 1, 4, and 3 (NASA/GSFC/ESDIS, 2018). Beside a cloud bank over the North and Baltic seas and fog over Southern and Western Germany, a considerable number of line-shaped contrails and diffuse cirrus clouds are present across both satellite images. Main contrail clusters are found over Central Germany for both situations. Contrails can also be identified over the Netherlands and Belgium, over the Czech Republic, and south of the Alps. At 12:00 UTC, the high-level cloud cover has increased and pervaded compared to 10:00 UTC.
Comparing Fig. 10a and e, the reference simulation clearly underestimates the coverage of high clouds in the center of the domain, which, in large parts, consists of contrails and contrail cirrus (Fig. 10c). Although a number of the patterns of high-level cloud cover detected by MODIS at 10:00 UTC can be identified in the aviation simulation, the amount of cloud cover seems to be slightly underestimated in Fig. 10c. This discrepancy is probably due to the fact that air traffic was switched on at 08:00 UTC and earlier flight movements are disregarded in our simulation. Hence, the simulation evaluation at 10:00 UTC neglects all contrails older than 2 h. The comparison at 12:00 UTC is more favorable than the one at 10:00 UTC. Observations match a lot better with the aviation simulation than with the reference simulation, as in the aviation simulation a much larger fraction of contrails could form since 10:00 UTC. Mostly over the center of the domain, areas with values of τ between 0.2 and 0.6 and peaks up to 1.0 are simulated. This is in good agreement with the very high values of the extinction coefficient of contrails compared to those of natural cirrus (Fig. 7j, k, l). The line-shaped patterns stem from the contrail ice class and the more patchy structures from aged contrails which have been transferred to the cirrus ice class. The case study for this particular day demonstrates that the inclusion of aircraft effects in a regional weather forecast model improves the representation of high-level clouds. Figure 11 shows the changes in the incoming SW radiation at two points in time at the surface. Clearly, additional ice crystals caused by air traffic reduce the incoming SW radiation. The line-shaped structure of young contrails at 10:00 UTC (Fig. 3a) is reflected in similar patterns of reduced SW radiation (Fig. 11a). As discussed earlier, the contrail ice class features numerous very small ice crystals with a high extinction coefficient. In the case investigated, this leads to a reduction in incoming SW radiation of 1 to 15 Wm −2 . As previously mentioned, our parameterizations are not able to calculate the optical properties of ice crystals with effective radii below the threshold of 5 µm. Rather, crystals with an effective radius of less than 5 µm are assumed to have optical properties of crystals as large as 5 µm. Therefore, the shading effect of the contrail ice crystals might be underestimated in our calculations.

Contrail impact on surface radiative fluxes
A strong spatial increase in the shading effect is found at 12:00 UTC (Fig. 11b). Here, the contrail coverage reaches its maximum. Still partly line-shaped contrails are found with a similar reducing impact as 2 h before. Additionally, clusters of contrail ice transformed into the cirrus ice class have evolved, particularly over the center of the domain (Fig. 4b).
The ice crystals herein are relatively small, but larger than those present in the contrail ice class. Therefore, they inhibit less SW radiation from reaching the ground than the ice crystals in the contrail ice class. Consequently, the reduction on SW radiation here is smaller, but still reaches 1 to 10 Wm −2 .
Especially in the north of the domain, small areas with the difference in incoming SW radiation attaining large negative values adjacent to large positive values are found. They occur when subtracting, for example, fields of radiative fluxes of the aviation simulation from those of the reference simulation. The introduction of contrail ice acts as source of disturbance for various processes like convection or turbulence. Those features are to be classified as noise as they do not influence the overall situation.
Non-negligible areas with an increase in SW radiation occur, e.g., at 12:00 UTC over the southeastern part of the domain (Fig. 11b). Besides only reducing the direct incoming SW radiation, contrail ice crystals also enlarge the flux of diffuse SW radiation (see below). On occasion, as the simulated contrails are still rather optically thin, this effect may be larger than the reduction of direct radiation.
The average change of incoming direct and diffuse SW for 3 December 2013 is shown in Fig. 12. Also here, the smallscale fluctuating values in the north are most likely noise.
The thin veil of contrail cirrus that spreads over most of the domain causes an average decrease of incoming direct and diffuse SW radiation of 1 to 5 %. The large and persistent contrail cluster over the northern and eastern part of the domain inhibits on average 5 to 10 % of SW radiation from reaching the ground during approximately 8 h of daylight. The reduction is strongest over the south and the west of the domain. Here, contrail ice is present for the longest time, as seen in Fig. 9. Figure 13 shows the temporal evolution of incoming SW radiation and normalized PV power. Panels a to d represent mean values for the entire simulation domain, panel e additionally represents mean values for an area of approximately 50 km × 50 km centered at a spot in the northeastern part of the simulated domain, indicated by the red circle (53 • N, 12 • E) in Fig. 12. This is the location of one of the largest solar parks in Germany (Solarpark Brandenburg-Briest). Also several other major solar parks are located around this area. For enhanced clearness, values for the sensitivity study (see Sect. 4.5) are depicted only as difference to the reference simulation. The differences between the aviation simulation and the reference simulation are in general more pronounced for the selected location that on average over the entire simulation domain.
During the whole time of daylight, contrails and contrail cirrus reduce up to 15 Wm −2 (about 7 %) of the total incoming SW radiation in the entire domain (Fig. 13a). The effect is largest during noon and ceases during the day. This corresponds to the size of contrail ice crystals. As in our simulation, contrails start to form at 08:00 UTC, the average contrail ice crystal size grows during the day. Accordingly, contrail ice effective radii also increase with time and lead to smaller values of the extinction coefficient.
Separating the total SW radiation into its direct (Fig. 13b) and diffuse (Fig. 13c) fraction, it is clear that especially the Figure 13. Temporal evolution for 3 December 2013 of incoming shortwave radiation: reference (black); aviation simulation as before (blue); "bio-fuel" scenario with EI iceno × 0.1 (orange), omission of crystal loss during contrail vortex phase with λ N s = 1 (red). Dashed lines (corresponding to right y axis) are difference to reference simulation. (a) total SW; (b) direct SW; (c) diffuse SW; (d), (e) normalized PV power. (a-d) Over entire domain; (e) mean for the location marked with the red circle in Fig. 12. Error bars represent mean values ± standard deviation. direct incoming SW radiation is strongly reduced by more than 20 Wm −2 due to the presence of additional ice crystals in the atmosphere. In contrast, the diffuse incoming SW radiation is increased by up to 10 Wm −2 . Enhanced scattering of SW radiation caused by the contrail ice crystals increases the diffuse SW radiation at the ground.
Notably, the peak in reduction of diffuse SW radiation occurs in the afternoon, around 13:00 UTC. Between 08:00 UTC and about 11:00 UTC, the amount of diffuse SW radiation reaching the ground is larger in the aviation simulation. During this time, as mentioned before, contrail ice crystals are smallest on average and forward scattering is less pronounced than for larger crystals, whereas contrail ice crystals grow on average during the day, resulting in enhanced forward scattering.
In the aviation simulation, young and aged contrails generally reduce the incoming SW radiation at the surface. This effect is currently neglected in operational weather forecast models. However, this effect is of relevance for the production of solar energy. The temporal evolution of the normalized PV power is depicted in Fig. 13d and e. The normalized PV is calculated using the open-source PV modeling environment PV_LIB for python (Andrews et al., 2014). For a specific combination of a PV module and a PV inverter combination, a nominal power and a reasonable tilt are assumed. These technical specifications are taken from Rieger et al. (2017). They assume panels consisting of a south-oriented PV module with a nominal power of 220 W and a size of 1.7 m 2 . Compared to the reference simulation, the normalized PV power is decreased in the aviation simulation most of the day. The largest losses of up to 10 % occur in the morning and diminish during the day; an increase even occurs for the selected location in the late afternoon (Fig. 13e). The normalized PV power is somewhat more strongly reduced than the total SW radiation. For production of PV power, the incoming direct radiation is of greater importance than the diffuse radiation; of the two, the direct experiences the larger reduction from contrails and contrail cirrus.
The error bars in Fig. 13d and e represent the mean values ± the standard deviation with respect to the entire simulation domain and the area of 50 km × 50 km around the selected location, respectively. The standard deviation in Fig. 13d reflects the large-scale variability of the impact of contrails and contrail cirrus on the incoming SW radiation, whereas Fig. 13e illustrates the small-scale variability.
Compared to the selected location, standard deviations are larger for the mean of the domain. Obviously, clouds modify the amount of SW radiation reaching the ground in a nonuniform manner. The magnitudes of the standard deviations for the aviation simulation are about the same magnitude as the ones for the reference simulation. Apparently, the impact of contrails and contrail cirrus on incoming SW radiation is as variable as the impact of natural clouds. For example, even at 12:00 UTC, when the impact of contrails and contrail cirrus is largest, confined areas exist that are unaffected by contrails and contrail cirrus (Fig. 11b).
However, the small-scale variability of the impact of contrails and contrail cirrus on SW radiation is rather small, reflected by much smaller standard deviations in Fig. 13e. One can therefore deduce that the exact location of contrails or contrail cirrus is not crucial for the strength of the impact on SW radiation.

Sensitivity to initial ice crystal number and early contrail ice crystal loss
In this last section, we briefly examine two sensitivities of our model setup.
For the first, we reduce the emission index for ice crystals EI iceno by a factor of 10 (orange lines in Fig. 13). This scenario explores lower engine soot emissions caused by either improved engine combustor technologies or fuel composition changes from, e.g., biofuel adoption (Moore et al., 2017). Due to this, the initial number concentration of contrail ice crystals is reduced, thus fewer but larger ice crystals are formed (Unterstrasser, 2014). In the simulation with reduced EI iceno , the reduction of total incoming SW radiation is only slightly weaker than for the simulation assuming standard fuel (Fig. 13a). As ice crystals are slightly larger, their extinction coefficient is lower and the reduction of direct SW radiation is smaller than for the standard setup (Fig. 13b). Also the increase in incoming diffuse SW radiation is less strong compared to the standard setup (Fig. 13c). Consequently, the reduction in normalized PV power is also less strong than for the standard setup, even though an enhancement occurs during afternoon (Fig. 13d, e).
Second, we set the surviving fraction of ice crystals λ N s in the parameterization of the initial ice crystal number to 1. This deliberately neglects the effects of crystal loss during the contrail vortex phase, as parameterized by Unterstrasser (2016). The influence of this parameter is large. During daytime, a reduction in total incoming SW radiation of up to 15 % is simulated (Fig. 13a). Both the reduction in direct SW radiation as well as the increase in diffuse SW radiation are much more pronounced for this case (Fig. 13b, c). Especially the reduction in direct SW radiation causes a strong reduction in production of PV power. Here, losses of nearly 15 % occur at about 10:00 UTC. Also concerning the temporal evolution, the reduction lasts much longer compared to the aviation simulation. When no early crystal loss is parameterized in the model, initial ice crystal number concentrations may be much higher than usual. As the initial IWC remains the same, the new crystals are smaller and, thus, the simulated contrails are optically thicker. The much stronger reduction in incoming SW radiation demonstrates that the early ice crystal number loss is non-negligible and an important aspect of contrail evolution as it has a long-lasting impact on contrail-cirrus radiative properties.

Conclusions
In this study, the regional atmospheric model COSMO-ART coupled with a two-moment microphysical scheme and a diagnostic treatment of radiation was extended by a parameterization describing contrails and the related physical processes.
Methods for a separate but consistent treatment of contrail ice were implemented to satisfy the special requirements describing the microphysics in young contrails and the transition phase to contrail cirrus. Performing a case study for a single winter day over Central Europe, it was shown how microphysical properties such as ice water content, ice crystal number concentration, and the mean ice crystal radius of ice crystals in contrails change over time and depend on the meteorological conditions. The ice water content in young contrails is comparable to that in thin cirrus clouds ranging from 0.2 to 5.0 mg m −3 , but with considerably higher ice crystal number concentrations between 1 and 100 cm −3 and effective radii below 10 µm. The numerous small ice crystals produce high values for the extinction coefficient and thus also for the optical depth.
The transition of contrail ice into the cirrus ice class causes increasing number concentrations. Here, effective radius of the ice crystals from the contrail ice class grows to an extent comparable to that of natural cirrus. Because of the still relatively high number concentrations, contrail cirrus still features high values for both the extinction coefficient and the optical depth.
Qualitative comparison with satellite data shows good agreement and proves advantages of considering contrails and contrail cirrus in a regional weather forecast model.
Contrail cirrus tends to cause changes in the microphysical appearance of high-level cloud coverage to a remarkable extent, which in turn influences the radiative effect in these regions.
In addition, the impact of contrails and contrail cirrus on SW radiation and the production of PV power were quantified. Although the case study was performed for 3 December 2013, when solar zenith angles are low and the length of days is short, a strong influence of contrails is still simulated. They inhibit up to 5 to 10 % of SW radiation from reaching the ground at noon. This results in a loss of PV power production of up to 10 %.
Moreover, it was demonstrated that ice crystal loss in young contrails is an important process which can significantly change the contrail-cirrus properties on a regional scale. This study is the first approach to simulate contrails and contrail cirrus using a numerical weather prediction model with high spatial and temporal resolution. Subsequently, the presented method can serve as a basis for improving the predictability of the solar radiation in regional weather forecasting by taking into account contrails and contrail cirrus. This section expands the description of the microphysical model from Sect. 2.1 and presents a collection of underlying equations. All formulae can be found in Seifert and Beheng (2006) as well. Values of constants used in this section are listed in Table A1. The generalized distribution is defined as Here, m is in units of kilograms. The parameters ν and µ are assumed constant, respectively. The parameter A is related to the total ice crystal number concentration n = M 0 and λ to the ice crystal mean mass m = q i /n = M 1 /M 0 . Expressions involving functions exist for the moments The growth equation of a single ice crystal is given by (Pruppacher and Klett, 1997) Here, T is the temperature, p is the pressure, and C = D/π denotes the capacity of hexagonal crystals (Harrington et al., 1995). S i is the supersaturation with respect to ice, L iv represents the latent heat of sublimation, p sat,i denotes the saturation vapor pressure over ice, R v is the specific gas constant for water vapor, K T is the conductivity of heat, and D v is the molecular diffusion coefficient of water vapor. The term F ven accounts for ventilation effects and G iv considers the diffusion of water vapor and the effect of latent heating: Integration of Eq. (A3) over the ice crystal mass spectrum yields the temporal derivative of the ice mass density q i : The mass-size relation given by Eq.
(2) uses the values a geo,nat = 1.59 and b geo,nat = 2.56 for the cirrus ice class. In order to avoid unreasonably small or large mean masses m = q i /n, a lower limit m min = 10 −12 kg and upper limit m max = 10 −6 kg are introduced. If m lies outside the interval [m min , m max ] in a grid box, the ice crystal number concentration is increased to q i /m min or reduced to q i /m max , respectively. The limits correspond to sizes L min = 17.5 µm and L max = 3800 µm. For the contrail ice class, a piecewise definition of a geo and b geo is employed following (Spichtinger and Gierens, 2009;Heymsfield and Iaquinta, 2000). For masses above m split = 2.15 × 10 −13 kg (corresponds to L split = 7.42 µm), a geo,con = 0.04142 and b geo,con = 2.2 are used. For masses below m split , a geo,con = 526.1 and b geo,con = 3.0 are prescribed, which define quasispherical hexagonal columns with aspect ratio 1 (see derivation in Spichtinger and Gierens, 2009). The latter constants are valid down to the prescribed lower limit m min = 10 −15 kg. For grid boxes with lower mean masses, the same bounding technique is used as in the cirrus ice class. The upper limit is set to a relatively small value of m max = 2 × 10 −11 kg and the treatment of grid boxes with too large mean masses is different compared to the cirrus ice class. Instead of bounding n, the total ice crystal mass and number from such a grid box are transferred from the contrail ice class to the cirrus ice class. The prescribed mass limits of contrail ice class correspond to the size limits L min = 1.24 µm and L max = 58 µm.
Similar to the mass-size relation, the terminal settling velocity v is approximated by a power law (Eq. A6) with coefficients a vel and b vel . v(m) = a vel m b vel (A6) Integrating v(m) and m v(m) over the ice crystal mass spectrum yields mean fall velocities v k for the ice crystal number (k = 0) and mass (k = 1), respectively:  (2) is beneficial. Then integrals over the mass distribution can often be expressed in terms of moments, which avoids employing more expensive numerical quadrature techniques. The contrail ice class uses a piecewise definition of the mass size relation. In this case, truncated moments have to be evaluated. In the model, both self-collection and collection between the individual classes of frozen hydrometeors are considered. Sink terms (not shown) are included in the prognostic equations of n. In case that two different hydrometeor classes are present, where class A has a larger mean mass than class B, the ice crystals of class B are collected by those of class A. In addition to the number loss in class B, mass is transferred from class B to class A.
The nucleation process for natural cirrus is not repeated here as ice crystal generation in the contrail ice class is quite different (see Sect. 3.1).  (2006)  f L (L) also represents a generalized distribution with parameters A L = A a ν+1 b; ν L = b(ν + 1) − 1; λ L = λa µ ; µ L = bµ.
Plugging Eq. (B2) into Eq. (3) and using the latter relations, the integrals in the effective radius definition can be reformulated in terms of moments.

Symbol
Definition Unit gamma function λ, λ L slope parameter of generalized distribution λ N s surviving fraction of ice crystals µ parameter of generalized distribution ν parameter of generalized distribution ρ, ρ ice air density, density of ice kg m −3 A, A L scaling parameter of generalized distribution m −3 kg −1 , m −4 a geo , a nat parameter in mass-size relation a vel parameter in fall speed relation b con , b nat exponent in mass-size relation b vel exponent in fall speed relation c p specific heat capacity of air collision efficiency for classes A and B -EI iceno emission index for ice crystals kg −1 f, f L number concentration size distribution m −3 kg −1 , m −3 m −1 F ven ventilation coefficient -I ice mass mixing ratio kg kg −1 q init "emitted" ice crystal mass concentration kg m −3 I 0 water vapor emission per flight distance kg m −1 q i ice crystal mass concentration kg m −3 K T conductivity of heat J m −1 s −1 K −1 L ice crystal size m L iv latent heat for sublimation J kg −1 L max maximum particle length m L min minimum particle length m M k kth power moment of f (x) kg k m −3 n ice crystal number concentration m −3 N number concentration kg −1 N BV Brunt-Väisälä frequency s −1 n init "emitted" ice crystal number concentration m −3 N s number of surviving ice crystals per flight distance m −1 N 0 number of produced ice crystals per flight distance m −1 p pressure hPa p sat saturation pressure hPa r radius µm r e effective radius µm R d specific gas constant for dry air J kg −1 K −1 R v specific gas constant for water vapor J kg −1 K −1 RH i relative humidity with respect to ice -S i supersaturation with respect to ice t time s T temperature kg v velocity m s −1 V volume m 3 b span wing span m