Process‐Based Simulation of Aerosol‐Cloud Interactions in a One‐Dimensional Cirrus Model

A new microphysical cirrus model to simulate ice crystal nucleation, depositional growth, and gravitational settling is described. The model tracks individual simulation ice particles in a vertical column of air and allows moisture and heat profiles to be affected by turbulent diffusion. Ice crystal size‐ and supersaturation‐dependent deposition coefficients are employed in a one‐dimensional model framework. This enables the detailed simulation of microphysical feedbacks influencing the outcome of ice nucleation processes in cirrus. The use of spheroidal water vapor fluxes enables the prediction of primary ice crystal shapes once microscopic models describing the vapor uptake on the surfaces of cirrus ice crystals are better constrained. Two applications addressing contrail evolution and cirrus formation demonstrate the potential of the model for advanced studies of aerosol‐cirrus interactions. It is shown that supersaturation in, and microphysical and optical properties of, cirrus are affected by variable deposition coefficients. Vertical variability in ice supersaturation, ice crystal sedimentation, and high turbulent diffusivity all tend to decrease homogeneously nucleated ice number mixing ratios over time, but low ice growth efficiencies counteract this tendency. Vertical mixing induces a tendency to delay the onset of homogeneous freezing. In situations of sustained large‐scale cooling, natural cirrus clouds may often form in air surrounding persistent contrails.


Introduction
Unlike other clouds, cirrus does not reflect much incident sunlight back to space but traps heat originating from the Earth's surface and atmosphere efficiently. Cirrus also influences stratospheric water vapor concentrations by freeze-drying air at entry to the stratosphere. Therefore, cirrus clouds are important regulators of our climate. Exactly how strongly cirrus influence Earth's energy balance is uncertain, in part because their radiative effect-a subtle balance between life cycle integrated, net shortwave, and long-wave radiative flux changes-is difficult to capture in global models. This means that predictions of how cirrus responds to global warming remain uncertain (Kärcher, 2017a). Without suitable countermeasures, the next three decades will witness an increasing atmospheric impact of aviation-induced contrail cirrus (Bock & Burkhardt, 2019), the most easily perceived evidence of high clouds perturbed by anthropogenic activities. Lack of understanding of how those human-caused ice clouds evolve in the atmosphere creates uncertainty in quantifying the overall climate impact of aviation (Kärcher, 2018).
Small-scale microphysical processes within cirrus must be understood well to predict their atmospheric effects with confidence. As demonstrated in this study, key processes during cirrus formation unfold on vertical scales of tens of meters and smaller. Those scales are too small to be resolved in global climate models with vertical resolutions of 100 m or larger. Therefore, aerosol-cirrus interactions must be parameterized in such models. Parameterizations representing ice crystal formation (Kärcher et al., 2006;Liu & Penner, 2005;Phillips et al., 2008) contain simplifications regarding heterogeneous ice-nucleating particles (INP) due mainly to the lack of constraints on INP composition and abundance at cirrus levels (Jensen et al., 2018) and the difficulty to measure number concentrations of small ice crystals accurately in cirrus (O'Shea et al., 2016). Importantly, the updraft speeds driving these parameterizations are only partly resolved in regional and global models (Kienast-Sjögren et al., 2015) and need to be parameterized as well (Kärcher & Podglajen, 2019).
The tight interaction between realistic dynamical forcing and the competition of aqueous aerosol particles and INP for water vapor during ice nucleation events has not been studied systematically on the process level, leaving an important gap in simulating indirect effects of INP on cirrus with confidence. The number 10.1029/2019JD031847 concentration of cloud ice crystals depends on the type of aerosol available for nucleation, the number of INP, the mode of ice nucleation, and the evolution of ice supersaturation in the vicinity of ice-forming particles. The latter is influenced by vertical air motions. Homogeneous freezing of aqueous aerosol particles in cirrus conditions as the most fundamental ice formation process has been well researched. Recent years have seen a major advance in quantifying vertical air motions, or cooling rates, due to mesoscale gravity waves in the tropopause region (Podglajen et al., 2016). At the same time, advances in airborne instrumentation have paved the way for improved characterization of INP and cloud ice properties (Baumgardner et al., 2012;DeMott et al., 2011) so that there is hope for better constraints on INP concentrations and properties in future cirrus simulations.
Many simulations of ice crystal formation in cirrus have been carried out with parcel models. At the same time, studies have pointed out the importance of ice crystal sedimentation in simulations of cirrus (Jensen & Pfister, 2004;Murphy, 2014;Spichtinger & Gierens, 2009). Sedimentation can be included in a parcel model framework only with assumptions. The growth of ice crystals by uptake of water vapor (H 2 O) is a fundamental microphysical process in any study of ice-containing clouds. In associated ice crystal growth theories (Lamb & Verlinde, 2011), the H 2 O deposition coefficient is introduced as a measure of growth efficiencies. It is well known that deposition coefficients are variable and low due to the difficulty of H 2 O to enter and irreversibly attach to an ice crystal lattice. Yet most cloud models employ constant (high) deposition coefficients. Laboratory studies move towards quantifying the limitations imposed by the kinetics of H 2 O attachment on ice crystal surfaces at cirrus temperatures (<240 K). This heightens the prospect for more reliable simulations of ice crystal evolution during and after nucleation events.
While knowledge of the shapes (habits) of ice crystals is key in retrieving cirrus properties from remote -sensing instruments, to date, there is no consensus on how cirrus ice crystal habits can be predicted. This issue is strongly tied to the inherent complexity of the physics of vapor growth-the prediction of ice crystal morphology requires sound knowledge of a plethora of processes controlling H 2 O surface attachment kinetics (Libbrecht, 2005). At the same time, in situ measurements reveal a rich suite of ice crystal habits in cirrus, depending on cirrus type and meteorological conditions in which these clouds form and mature (Lawson et al., 2019), calling for better understanding and explanation.
Arguably, it is not necessary to leap from a parcel model framework to the full complexity of three -dimensional cloud-resolving models, as the task of better understanding and quantifying aerosol-cirrus interactions does not immediately demand this expansion. This work introduces a novel one-dimensional (1-D) cirrus model augmenting existing 1-D models (Cirisan et al., 2014;Jensen & Pfister, 2004;Kärcher, 2005;Lin et al., 2005;Murphy, 2014;Reichardt et al., 2008). We make use of a kinetic (non-equilibrium) treatment of the interaction between water vapor, aerosol particles, and ice crystals and a particle tracking approach for the ice phase. The latter offers a number of advantages over Eulerian microphysics schemes in cloud simulations (Grabowski et al., 2019). We predict deposition coefficients that depend on ice crystal size and respond to changing environmental conditions, most notably, ice supersaturation. While keeping the spherical approximation for small ice crystals, we employ spheroidal vapor fluxes towards their surfaces to allow for primary habit prediction once advanced models for H 2 O attachment kinetics valid in cirrus conditions become available. Previous studies have established that the spherical assumption for ice crystals is capable of reproducing the overall H 2 O mass uptake rate fairly well (Zhang & Harrington, 2014) but is clearly poor as far as the prediction of radiative properties is concerned. We account for the effect of turbulent diffusion on vertical profiles of water vapor and heat, which may be particularly relevant for young aircraft contrails.
We opt for a 1-D model framework because capturing vertical variations is more important for short-term process studies of aerosol-cirrus interactions than the simulation of horizontal variability. In such studies, it is possible to ignore effects of vertical shear of the horizontal wind field. For long simulation times, column simulations may still be useful representing horizontal averages of cloud variables. We ignore radiation-cirrus interactions that constitute an important physical component in, for example, tropical anvil cirrus (Gasparini et al., 2019), as our simulations address brief episodes of ice formation and growth processes in shallow, mostly optically thin cirrus columns. For the same reason, we work with spatially homogeneous wind speed profiles since it is more important to accurately capture the temporal wind speed forcing than to exactly balance divergent and nondivergent wind components. Moreover, we work with Figure 1. Salient features of, and physical processes simulated within, the one-dimensional cirrus model primeice1D, developed to simulate the formation of pristine (primary) cirrus ice crystals from aerosol precursor particles. The latter are characterized by a size distribution with total number concentration (n aer ) and, in each size category, volume mixing ratios of their (insoluble) core material, liquid water content (q lw ), and hygroscopicity parameter ( ). In a vertical column of air, the prognostic meteorological variables, potential temperature ( ), and water vapor mixing ratio (q wv ) are driven by a prescribed vertical wind field (w) and optionally turbulent diffusivity (). Air pressure (p) and temperature (T) are inferred based on hydrostatic equilibrium. The vertical coordinate and time (their numerical increments) are denoted by z (Δz) and t (Δt), respectively. Besides a slowly varying mean component (w), the wind forcing may include stochastic gravity wave-driven perturbations (w ′ ). Water substance is conserved locally, at each vertical station, between gas, aerosol, and ice phases. Growth and evaporation of supercooled liquid aerosol particles are calculated in a Lagrangian manner. They may freeze homogeneously depending on T and q lw . The ice phase is simulated by tracking the history of a large number of simulation ice particles (NSIP). The latter contain a number of real ice crystals; they are represented by spheroids as a surrogate for planar or columnar crystal shapes characterized by semiaxis lengths (a, c), number concentration (n), and ice water mixing ratios (q iw ). Deposition growth allows for axis-dependent water vapor fluxes and supersaturation-dependent deposition coefficients. Ice crystal deposition/sublimation rates and fall speeds (V) are size-and shape-dependent. No equilibrium assumptions are made to simulate the water content in aerosol and ice particles and their interaction with the vapor. microphysical assumptions appropriate for sub-100-m ice crystal sizes, as mean maximum dimensions of ice crystals during and shortly after in situ formation events do not exceed several 10 m. After describing the model physics (section 2), we present first applications of the short-term evolution of a persistent aircraft contrail and of a homogeneous freezing event during cirrus formation (section 3). The study concludes with a summary and outlook (section 4).
The microphysics formulation in primeice1D allows for a size-resolved, fully kinetic treatment of interactions between H 2 O, aqueous aerosol particles, and ice crystals. No thermodynamic equilibrium assumptions are made to simulate the supercooled liquid water or ice water content. The dry cores of aerosol particles contain soluble and optionally insoluble matter. The aerosol liquid water content is affected by the solubility parameter, . Soluble particles expand or shrink in size by adjusting their water content to variable relative humidity conditions (hygroscopic growth/evaporation and condensational growth/evaporation below and above liquid water saturation, respectively). In this way, aerosol particles can be activated into cloud water droplets when the relative humidity surpasses liquid water saturation. Water vapor uptake and the reverse process are treated individually for ice crystals (depositional growth/sublimation). The physics of water vapor uptake at ice crystal surfaces; hence, the effect on crystal dimensions and habits is embodied in molecular deposition coefficients, . Ice crystal sizes along with their number concentration, n, determine the ice water volume mixing ratio, q iw . When ice crystals are produced by nucleation, aerosol core material is returned to the respective size distribution upon full sublimation of ice crystals. This core return makes recycled aerosol particles available for the simulation of subsequent nucleation processes.
The ice phase is modeled by tracking simulation ice particles (SIPs), each representing a certain number of real ice crystals resulting, for example, from nucleation. SIP positions change according to w and their terminal fall speeds (sedimentation velocities), V. Each SIP is tracked individually by integrating its equation of motion over time. Ice crystals grow by uptake of H 2 O, and when their fall speeds are larger than the local updraft speed, they begin to settle out of the growth regions. In each SIP, all real ice crystals have the same physical properties. The total number of SIPs in a simulation, NSIP, is affected by nucleation and sublimation processes. SIPs are not allowed to fall across contiguous altitude bins within a time step when their interaction with the gas phase is fast. Individual SIPs may also be initialized with predefined properties at specified altitudes and can be associated with attributes recording their spatiotemporal evolution.
When particles fall in air, they quickly achieve a constant (terminal) fall speed when their gravitational acceleration is balanced by a hydrodynamic drag force. We model the fall speeds of pristine ice crystals employing a capacitance-based formulation (Westbrook, 2008). The capacitance, C, captures both size and habit of ice crystals. As C also affects H 2 O fluxes, depositional growth, habit evolution, and sedimentation are treated in a physically self-consistent manner. The restriction to sub-100-m sizes removes the necessity to include ice crystal aggregation, which is an important growth process in anvil cirrus (Gasparini et al., 2019) and some midlatitude cirrus (Sölch & Kärcher, 2010). For the same reason, it is not required to account for ventilation effects or to introduce a size-dependent ice crystal mass density, both modifying growth rates. Aerosol particles and the cloud droplets evolving from them are not allowed to sediment in primeice1D owing to their small sizes and rapid homogeneous freezing rates, respectively.
Quasi-spherical particles together with hexagonal plates and columns are frequently observed habits of sub-100-m ice crystals in many cirrus types (Lawson et al., 2019). Spheroids are therefore used as a first-order characterization of the habits of model ice crystals for simplicity and generality. Two semiaxes, a and c, define the aspect ratio, = c∕a, and therefore model ice crystal habits as oblate ( < 1) and prolate ( > 1) spheroids, approximating plates and columns of real ice crystals, respectively. We use the theory developed by Zhang and Harrington (2014) to simulate kinetically limited H 2 O fluxes towards ice crystals along a (prism faces) and c axes (basal faces), potentially enabling habit prediction. The latest version of this method, Diffusion Surface Kinetics Ice Crystal Evolution (DiSKICE), is described in Harrington et al. (2019) and allows to incorporate crucial physical mechanisms and dependencies that are entirely absent in conventional uptake models. However, other shapes, most notably bullet rosettes, can be found dominating the habit distribution below maximum crystal dimensions 50-100 m in some types of cirrus (Lawson 10.1029/2019JD031847 et al., 2019. In such cases, habits may at this point be diagnosed from single ice crystal mass-dimension relationships. A microscopic model for axis-dependent H 2 O deposition coefficients must be specified, determining ice crystal growth rates. At cirrus temperatures, uptake models are associated with uncertainty. Normal growth of ice crystal facets by dislocations (Burton et al., 1951) and stacking faults or two-dimensional ledge nucleation (Nelson & Baker, 1996) is taken into account in DiSKICE, approximately explaining the generation of primary ice crystal habits and other growth features. A ledge defines a region on an ice crystal face abundant of sites for preferred molecular attachment. Nelson and Swanson (2019) discuss less well-explored lateral growth modes leading to complex (secondary) habits, but a theory considering these growth modes is not yet available. While laboratory studies hint at columns and rosettes as prevalent habits depending on supersaturation at temperatures below −40 • C (Bailey & Hallett, 2009), exactly which ice crystal growth mode prevails in cirrus clouds remains uncertain, and associated laboratory results constraining relevant model parameters (critical surface supersaturation) are not available.

Initialization
We set up an altitude grid by specifying column base and top altitudes, z bot and z top , respectively, along with the number of vertical grid points, NZ, using a constant grid spacing, Δz = (z top − z bot )∕NZ. Initial vertical profiles, T(z), (z), and q wv (z), or alternatively, ice saturation ratio (fractional relative humidity over ice), S(z), may be taken from observations or numerical models. Air pressure is diagnosed from the adiabatic relationship p = p o (T∕ ) c p ∕R air , with the surface reference pressure p o = 1, 000 hPa, the specific gas constant for air, R air , and the specific heat of air at constant pressure, c p . Given the H 2 O saturation vapor pressure over ice, p sat (Murphy & Koop, 2005), we set q wv = S p sat ∕p. While overall air mass conservation places a restriction on the spatial distribution of the vertical wind, initial profiles, w =w + w ′ , are assumed to be spatially homogeneous. We prescribe a constant mean updraft speed,w > 0, and may add spatially uncorrelated fluctuations, w ′ (t), due to gravity waves.
In the simulations analyzed in this study, we specify initial values, T(z bot ) and p(z bot ), along with the environmental lapse rate, . Then, T(z) and (z) follow from integrating dT∕dz = − and d ∕dz = (Γ − ) ∕T, respectively, with the dry adiabatic lapse rate, Γ = g∕c p ≥ ; g is the acceleration due to gravity. Moreover, we construct a profile, S(z), with well-defined ice crystal nucleation, growth, and sublimation zones by means of a Gaussian distribution. While this idealized, single-mode S(z) profile suffices for the purpose of this study, upper tropospheric water vapor profiles, for example, at midlatitudes, can be highly variable (Flentje et al., 2005), introducing more vertical structure in S(z).
In simulations accounting for vertical turbulent diffusion, we prescribe a constant H 2 O diffusivity in air, . The Lewis number, Le, defined as the ratio of thermal to mass diffusivity, is used to estimate the associated thermal diffusivity affecting and T. Here, we use Le = 1. Particles are not affected by diffusion, assuming that they do not completely follow turbulent velocity fluctuations due to their inertia. As turbulence levels are typically low in the upper troposphere, we focus in the case of ice crystals on tracking sedimentation and capturing effects of water vapor and heat diffusion on nucleation.
During nucleation events in cirrus, ice crystals are created from aerosol precursor particles. We have implemented aqueous aerosol particles containing supercooled liquid water, solutes, and optionally insoluble cores as internal mixtures. The soluble core fraction dissolves in the presence of condensed water creating hygroscopic solution droplets, which may activate into nearly pure water droplets by condensing large amounts of supersaturated vapor. Aerosols do not coagulate among themselves and are not scavenged by ice crystals. Instead, they may freeze homogeneously at high ice supersaturation, typically S > 1.5 but below liquid water saturation. In future simulations, we may associate internally mixed insoluble cores with immersion freezing nuclei and further add externally mixed INP.
Aqueous aerosol particles and their dry cores are assumed to be spherical. To initialize mixed aerosol particles, we prescribe a size distribution of the number concentration and volume mixing ratio of the dry cores at each vertical station. The distribution is divided into a number of size categories around the modal particle size. The number-based mixing ratio in a given size category is q aer = n aer (Δr)∕n air , where is the fractional number concentration of particles initially residing in a given radius range, Δr, and n aer is the total number concentration of aerosol particles. In this study, follows a log-normal distribution. The dry core volume mixing ratio follows from q core = V core q aer , with the spherical core volume V core = 4 r 3 core ∕3. During KÄRCHER 5 of 21 10.1029/2019JD031847 a simulation, particles from each size category are allowed to grow to their exact total radii, r w , determined by predicting the supercooled liquid water content mixing ratio, q lw , for each particle size category. The latter is defined via q lw = (V w − V core )∕ · q aer , with the total liquid particle volume, V w = 4 r 3 w ∕3, and the volume of an H 2 O molecule in liquid water, .

Vertical Advection and Turbulent Diffusion
The sequence of processes determining the overall fluid-dynamical and thermodynamic state of the air column is implemented as follows. After updating w(t), we advect and diffuse the potential temperature along the column by solving The diabatic source term,  , specified in section 2.6.4, accounts for the transfer of latent heat during water phase changes.
We employ a semi-Lagrangian solver (Xiao et al., 1999) to integrate the linear hyperbolic advection equation (equation (1) without the diffusion and latent heat terms), supplemented by Neumann boundary conditions. Its ability to suppress oscillations near steep gradients is especially valuable to model ice supersaturation profiles with narrow nucleation and sublimation zones.
The assumption of spatially constant  leads to a linear parabolic diffusion equation (equation (1) without the advection and latent heat terms). After discretization using centered differences on the fixed z grid, the resulting fully implicit difference equations are cast into a system of linear equations in tridiagonal matrix form, which can be solved very efficiently (Press et al., 1986).
Changes in due to latent heat exchange are calculated explicitly after the associated microphysical source and sink terms have been estimated (equation (1) without the advection and diffusion terms). The latter is accomplished after the new pressure profile is diagnosed by integrating with the Exner function, Π. Equation (2) is based on hydrostatic balance, z p = −g air , with the mass density of air, air = p∕(R air T). After solving equation (2), both and p are known, and the temperature profile is The gas phase water vapor mixing ratio is also subject to advection and diffusion and evolves according to with the microphysical source and sink term  q , specified in section 2.6.4. For the water vapor tracer,q wv , we solve equation (3) The Courant-Friedrichs-Levy criterion is used to limit the time step, Δt, for advection. Turbulent diffusivities up to 0.6 m 2 /s have been reported to characterize conditions in aging aircraft exhaust plumes (Schumann et al., 1995). Values representative for the background upper troposphere are considerably smaller. In high-resolution simulations, the vertical diffusion time scale, ( z) 2 ∕, becomes very short for small spatial scales, for example, about 1 min for  = 0.3 m 2 /s and z = 4 m. Time steps are not allowed to exceed the above diffusion time scale and are further constrained such that  q does not change q wv more than a small fraction of its previous value at any vertical station. Figure 2 shows idealized initial profiles based on p bot = 300 mb, T bot = 230 K, and = 7 K/km. The altitude domain spans 1.5 km, from z bot = 9 km to z top = 10.5 km; the domain is discretized into NZ = 1, 000 grid points, resulting in a spatial resolution of Δz ≃ 1.5 m. The width of the Gaussian ice saturation ratio profile, centered at 10 km, is 250 m. The entire domain is lifted withw = 5 cm/s, and we set w ′ = 0.
The profiles at t = 1.5 hr are also shown in Figure 2. Advective cooling causes the ice saturation ratio to increase over time, creating a ≈ 400-m thick, strongly ice supersaturated layer. The increase in S near the peak is larger than in the region below, owing to the lower temperature in the upper layers. The advected and S profiles are not significantly affected by vertical diffusion (not included) owing to their shallow gradients. Perturbations, w ′ , due to mesoscale gravity waves and resulting adiabatic temperature fluctuations are not included in this work, but the methodology to generate them is discussed in Kärcher and Podglajen (2019).

Contrail Simulations
primeice1D may be applied to simulate the evolution of persistent aircraft contrails and contrail cirrus clouds. This ice cloud type is only poorly researched at the process level. A number of atmospheric processes control the microphysical and optical properties of these ice clouds in the spreading stage that commences about 10 min past formation behind the jet engines of cruising aircraft (Kärcher, 2018). Over several hours of lifetime, contrails contain more smaller ice crystals than natural cirrus (Voigt et al., 2017). Peak ice crystal number concentrations in fresh (age < 10 min) contrails are at least two orders of magnitude larger, and associated mean ice crystal diameters are about tenfold smaller than those typically found in mature midlatitude cirrus clouds.
At the beginning of the spreading stage, fluid-dynamical perturbations caused by dissipating aircraft wakes have ceased, and further evolution of the contrail is governed by atmospheric processes. A combination of shear and turbulent diffusion leads to mixing of environmental air into the contrail volume, resulting in a continuous reduction (dilution) of the number concentrations of contrail constituents. We parameterize the effect of entrainment by multiplying at each time step the contrail ice crystal number concentrations by the incremental dilution factor, (t) = [(t 0 + t)∕(t 0 + t + Δt)] b , with the initial contrail age t 0 and t ≥ 0. The power law exponent b ≈ 1 may vary slightly depending primarily on wind shear and turbulence levels.
Consistent with the above mixing law, environmental air is entrained at a rate = b∕(t 0 + t), with t ≥ 0. This process is taken into account at every vertical station for the contrail variables = {T, q wv } according to d ∕dt = − ( − env ), where env (z, t) is calculated based on equations (1) and (3) setting  =  q = 0, respectively. This entrainment law is similar to a traditional nudging equation with a time-dependent nudging time scale, 1∕ , and assumes that the environmental air undergoes the same meteorological evolution as the contrail but without microphysical processing. Entrainment becomes very slow in old contrails.
There is little consistent information regarding which ice crystal habits prevail as contrails evolve. Shapes of ice crystals in a young contrail were observed to become aspherical few minutes after formation, and the degree of deviation from initial isometry varies (Gayet et al., 2012). Some aircraft measurements reported plate-like and columnar crystal with sizes up to about ≈ 20 m in the core region of a contrail and larger 10.1029/2019JD031847 crystals (> 100 m) at its periphery (Lawson et al., 1998). One study found that 75% of all ice crystals in a cold contrail at 212 K were hexagonal plates (Goodman et al., 1998). Another study found a predominance of small (≈ 10 m) plate-like ice crystals in a young contrail (Jones et al., 2012); occasionally, larger crystals (≈ 100 m) in aging persistent contrails were detected in these aircraft measurements.

Aqueous Aerosol Microphysics
The treatment of aerosol microphysics follows closely the description in Kärcher (2017b), and we only briefly summarize the main features below. Aqueous particles contain inert soluble and optionally insoluble core components. As a two-moment scheme, n aer and q core together with predicted liquid water volume mixing ratio, q lw , it allows us to estimate the exact total size of aqueous particles from each size category as a function of time. These three variables are advected separately and further change due to ice nucleation; in addition, q lw increases (decreases) by condensation (evaporation). Hygroscopic and condensational growth are simulated, employing a parametric water activity model based on Köhler theory (Petters & Kreidenweis, 2007). The associated mass growth rates for particles in each size category, q lw ∕ t, are predicted from a diffusion equation with gas kinetic corrections. The amount of supercooled water is used to calculate the sizeand solubility-dependent water activity, If soluble cores are composed of multiple chemical species, then denotes the average core hygroscopicity; for insoluble cores, = 0. The liquid water saturation vapor pressure over solution droplets of a given size depends on the amount of water dissolved in them via A w . The mass balance between water vapor and supercooled liquid water is solved using the analytical predictor of dissolution scheme (Jacobson, 1999). Its solution returns at each vertical station a balanced set of gas/aerosol variables, q wv and q lw , from which the liquid water volume per aerosol droplet, V lw , and the wet particle radii are evaluated: Unless aerosols activate into nearly pure water (cloud) droplets, which occurs in detrainment zones of convective clouds but usually not during in situ cirrus formation, the amount of liquid aerosol water is much smaller than that in the gas phase.
At each time step at a given altitude, homogeneous aerosol freezing rates are estimated from temperature and water activity using the rate coefficient, J(A w , T) (Koop et al., 2000); J increases dramatically with decreasing T and is extremely sensitive to even small changes in A w (Murphy, 2014). The fraction of aqueous particles freezing homogeneously in a time step is given by 1 − exp(−V lw JΔt), determining the nucleated ice crystal number concentration, n, and the loss fraction of q aer , q core , and q lw .

SIPs
The total number of SIPs, NSIP, increases due to nucleation or when SIPs are injected into the model column with predefined properties mimicking ice crystal sources that are not represented in the model. We monitor the number of SIPs injected, nucleated, sublimated, and pushed or fallen out. The upper domain boundary, z top , will have to be increased should SIPs be pushed out at the cloud top during a simulation. NSIP decreases when SIPs settle through the lower domain boundary, z bot . We implement a randomized procedure to generate SIPs by homogeneous freezing (Unterstrasser & Sölch, 2014), resolving a minimum nucleated number concentration, n min . Lowering n min leads to a better resolution of homogeneous freezing events driven by low cooling rates and improves the representation of ice crystal size distributions but increases NSIP. On the one hand, average ice crystal properties in each grid cell containing NSIP particles converge ∝ 1∕ √ NSIP. Large NSIP increases the simulation time, on the other. The ideal value of NSIP is a trade-off between statistical accuracy and computational expense.
Essential SIP properties tracked across their life cycles include c, a, and n; the SIP location, z SIP ; and the radius, r core , of the aerosol particles on which ice crystals form. In addition, we may add any number of attributes for diagnostic purposes, such as the time of SIP generation or the local ice supersaturation at z SIP . After full ice crystal sublimation, SIPs are removed, and the associated values n and r core are added back to the aerosol grid modifying q aer and q core at the corresponding SIP location, respectively.
SIPs contain a certain number of real ice crystals. For instance, in a grid cell of thickness Δz = 1 m, a total number concentration of n = 1 cm −3 by volume corresponds to nΔz = 100 ice crystals per square centimeter of air. We work with the assumption that all crystals in a given SIP have the same physical properties and reside at the same predicted location. This assumption is reasonable, as SIPs are generated within short time steps (to accurately resolve ice nucleation events) during which the conditions controlling their initial properties do not vary much. Nonetheless, these properties tend to spread over time for any given SIP. To account for the effect of decorrelation, we might assign probability distributions to each SIP property once robust, predictive equations governing such distributions are available.
Advection and sedimentation cause SIP locations to change across the z-grid: While advection leaves the ice crystal number mixing ratio, = n∕ air , constant, sedimentation leaves n constant. We apply advection and sedimentation sequentially and adjust n after the advection step, followed by sedimentation. After locating the grid cells in which SIPs reside, we linearly interpolate values of gridded variables to the individual SIP locations to calculate process rates affecting SIPs. By limiting the time step (cf. section 2.6.4), we ensure that SIPs reside within a grid cell long enough to be able to release or deposit water vapor there.
We recall from section 2.1 that individual crystals are characterized by the semiaxes, a and c, of oblate and prolate spheroids, in which case analytical expressions for the ice crystal capacitance factors, C(a, c), are available (Lamb & Verlinde, 2011). Capacitances revert to volume-equivalent radii in the case of spherical particle geometry (c = a = D∕2). Terminal fall speeds are calculated via V = mg∕(6 C) (Westbrook, 2008), with the dynamic viscosity of air, (T). The local ice water content, M(z), is obtained from summing the ice mass concentrations per unit volume of air over all SIPs present in a given grid cell.

Water vapor fluxes
We use the spheroidal fluxes (per unit ice crystal surface area and unit time) towards the prism (x = c) and basal (x = a) faces of ice crystals from Zhang and Harrington (2014): with the kinetically limited diffusion coefficients In the above, C is the capacitance taken at distances a + and c + above the ice crystal surface with the vapor jump distance, , approximately equal to the mean free path of H 2 O in air; D wv is the diffusivity of H 2 O molecules in air; n sat = p sat ∕(kT) is the ice saturation vapor number concentration of H 2 O taken at the environmental temperature; T, s = S−1 is the associated ambient ice supersaturation; x are axis-dependent deposition coefficients; and = 4D wv ∕v is the diffusion length scale with the mean thermal speed of H 2 O molecules,v. The factors x are related to the surface resistance to H 2 O uptake, matching the transition from free molecular towards diffusion-limited growth where the impact of this resistance diminishes. As the fluxes do not contain angular dependencies, they basically model faceted growth.
Latent heat released during vapor growth increases the ice crystal surface temperature; hence, the local supersaturation directly above the crystal surface. This reduces further uptake of H 2 O. The factor in equation (7) accounts for diffusive transport of latent heat towards or away from the surface. In our notation: 10.1029/2019JD031847 with the kinetically limited H 2 O diffusivity (Zhang & Harrington, 2014) In equation (10), m wv is the mass of a water molecule,  is the enthalpy of sublimation, R wv is the gas constant for H 2 O, and K ⋆ is the kinetically limited thermal conductivity of air defined analogously to equation (11); in K ⋆ , we assume perfect thermal accommodation for each axis, as heat is transported efficiently towards or away from ice crystals by abundant air molecules. The approximate treatment of thermal effects underlying allows us replace the unknown ice crystal surface temperature with the ambient temperature in equations (7)-(11).

Deposition Coefficients
The above set of equations must be closed by an H 2 O uptake model providing x . Deposition coefficients contain information about the local ice supersaturation, s x , directly above the crystal surface that basically controls the kinetics of H 2 O attachment: wherêx follows from equation (9) but with C replaced by C (Zhang & Harrington, 2014). Two axis-specific factors in equation (12) modify the ambient supersaturation, s, driving deposition growth. The first accounts for the difference between ambient and ice crystal temperatures due to latent heating/cooling. The second accounts for the diffusional resistance of H 2 O moving in air. Both factors depend on x . Only one other published cirrus model couples surface supersaturation to an approximate model for surface layer nucleation, for spherical particles, and without latent heat correction (Cirisan et al., 2014).
Available attachment kinetic models describe faceted growth, either by dislocations or ledge nucleation (Nelson & Baker, 1996;Zhang & Harrington, 2014), with x values estimated from the transcendental equation and s x taken from equation (12); s cr,x (T) > 0 is a critical surface ice supersaturation, and is an adjustable parameter inferred from laboratory measurements (Harrington et al., 2019). Equation (13) is iterated twice for each SIP, yielding c and a depending on ambient conditions (p, T, s), particle dimensions (a, c), and attachment kinetics of H 2 O ( , s cr,x ).
Deposition coefficients obtained from equation (13) decrease with decreasing ambient supersaturation (s). Diffusional resistance diminishes the surface supersaturation relative to s as the pressure (D wv ) decreases and as ice crystals grow to larger sizes (C). High s cr,x values hinder growth, and increasing values interpolate between more efficient ( = 1 − 10, crystallographic defects) and less efficient ( = 10 − 30, ledge nucleation) H 2 O uptake. For large , uptake is basically shut off below s cr,x , indicating the sharp threshold behavior of the surface nucleation process. Treating as a constant assumes that the number of nucleated ledges on the crystal surface does not change with s.
We recall that while the model is formulated to solve the growth evolution of ice crystal habits, this cannot be adequately addressed since s cr,x values have not yet been reliably determined from laboratory measurements for cirrus conditions. We use the DiSKICE model to estimate a single critical supersaturation as a credible baseline employing crystal-averaged values, s cr = s cr,c = s cr,a , that is, identical yet variable deposition coefficients for both axes, = c = a . While values calculated in this way respond to changing environmental conditions, this assumption implies that is taken to be uniform over the ice crystal surface: the deposition coefficients for the basal and prism faces are assumed to be equal, preventing habit evolution. Consequently, spherical ice crystals remain spherical in the growth phase after nucleation. This is a sensible approach as past studies have shown that the mass uptake of nonspherical ice crystals is quite accurately reproduced in this way (Zhang & Harrington, 2014).
The incorporation of adsorbed H 2 O molecules into the ice crystal lattice derives from the molecular dynamics and structure of the crystal surface, which may contain chemically inert or reactive impurities 10.1029/2019JD031847 (Huthwelker et al., 2006). Little is known about how foreign (non-H 2 O) molecules, possibly dissociated into ions, interact with active ice crystal growth sites in atmospheric conditions. While the presence of surface impurities needs to be considered in laboratory measurements of ice crystal growth, this might be less problematic in the atmosphere due to the rather low abundance of potential coadsorbed species (such as HNO 3 and HCl) relative to H 2 O. However, as cirrus ice crystals may form from aqueous solution droplets, solute molecules may be expelled to, and retained at, the ice crystal surface during freezing, eventually blocking growth sites at the surfaces of small ice crystals.
During diffusion-limited ice crystal sublimation, we follow Harrington et al. (2019) and treat the aspect ratio as a constant by setting deposition coefficients to unity for both axes due to the effect of surface roughening. This preserves the habits of ice crystals they possess at the point where sublimation commences.

Ice crystal Growth
We first note the growth/sublimation equations that hold for the two axes in each SIP: with the volume of an H 2 O molecule in ice, = m wv ∕ , and the constant bulk mass density of solid ice, .
When summed over all SIPs present in a given altitude bin, the overall H 2 O mass uptake rate follows from representing the classical diffusion equation for nonspherical particles (Lamb & Verlinde, 2011). The local time scale for quenching of ice supersaturation is given by If fall speeds exceed updraft speeds, SIPs may travel across several grid cells within a time step since, unlike for advection, their motion is not constrained by a numerical stability criterion. We compare the time scale for supersaturation quenching from equation (16) with the average SIP residence time, t res = (Δz∕2)∕|w−V|. We limit the time step by t wv but only if the local interaction of H 2 O with the ice phase is faster than the residence time.
The net gain or loss of water vapor in each grid cell for use in equation (3) follows from the net microphysical source/sink term for used in equation (1) is given by The molecular mass ratio = m wv ∕m air appears in equations (17) and (18), as we have defined q values as number-based (instead of mass-based) mixing ratios.
The nonequilibrium partitioning of H 2 O between the vapor and ice phases is obtained by solving equations (14)-(17) at every vertical station. The solution yields water vapor and ice water volume mixing ratios, q wv and q iw , respectively. The latter is defined for each SIP by High vertical resolution of water vapor is needed to realistically simulate homogeneous freezing at the tops of continually forced cirrus (Lin et al., 2005). Moreover, we note that accurate simulations of homogeneously nucleated ice crystal concentrations based on particle tracking are sensitive to the resolution of small concentration values. This typically results in second to subsecond numerical time steps.

Model Applications
We illustrate the performance of our model with high resolution (Δz ≈ 1.5 m) simulations, addressing the interplay between deposition growth and sedimentation in a contrail (section 3.1) and additionally homogeneous freezing in cirrus (section 3.2). Sensitivity studies are performed to explore effects of diffusion, entrainment, and variable (low) deposition coefficients. As motivated in section 1, we represent ice crystals as volume-equivalent spheres to estimate their effect on H 2 O mass uptake.

Growth-Sedimentation (Contrail) Simulations
We simulate the evolution of a contrail 10 min past formation for 1.5 hr, covering the initial phase of the spreading stage.

Model Setup
We use the idealized meteorological setup and vertical grid underlying Figure 2 to initialize growth -sedimentation simulations, with H 2 O limited to ice saturation (S = 1) in the contrail layer. Guided by observations of young contrail behind a medium-sized aircraft , we initialize a contrail area across a 220-m thick layer (149 grid cells) between altitudes of about 9.9 and 10.1 km. We sample spherical SIPs in the diameter range 1 − 10 m from a log-normal size distribution with modal diameter of 2 m and geometric standard deviation of 1.5. The choice = 1 is supported by the fact that contrail ice crystals form from freezing of supercooled water droplets. We spread 25 SIPs randomly in each vertical bin where contrail ice is present, yielding NSIP = 3, 725. The initial total contrail ice crystal number concentration is largest (100 cm −3 at t 0 = 10 min) in the upper layer that is 50 m thick, decreasing exponentially downwards with a scale height of 100 m. The reduced number concentration in the lower contrail area accounts for sublimation losses that occur in aircraft wakes within the first few minutes past formation . The n profile represents a spatial average across the aircraft wake perpendicular to the flight direction. The parameter governing the entrainment rate is b = 0.9.
We prescribe perfect molecular attachment, = 1, to maximize the H 2 O uptake rate (case HI-). This option fully accounts for diffusive limitations to ice crystal growth but neglects surface kinetic limitations. We perform two additional simulations relative to this case to judge effects of ignoring entrainment on moisture and temperature distributions inside the contrail (but keep the effect of dilution on ice crystal number concentrations) and adding vertical diffusion with a turbulent diffusivity,  = 0.3 m 2 s −1 . Laboratory measurements of single ice crystal growth carried out at T and S similar to our simulation suggest low deposition coefficients, = 0.006 ± 0.002 (Magee et al., 2006). To address the issue of low growth efficiency, we discuss a simulation with = 0.006 (case LO-). We comment a posteriori on the limitations of using constant deposition coefficients.

Results and Discussion
We first discuss the HI-case in Figure 3 (black curves). The contrail is lifted by almost 300 m over the course of the simulation (1.5 hr). Dilution causes initial ice crystal number concentrations to decrease about tenfold over the lifting period. At the same time, the ice crystals grow in size due to the ice supersaturation generated in the updraft. Throughout the contrail core, the high ice crystal number concentrations in combination with rapid uptake of H 2 O result in very low supersaturation (displayed here at SIP locations) and rather uniform ice crystal growth. A steady-state supersaturation profile establishes in the contrail core, resulting from local balance between diffusive vapor loss, adiabatic cooling, entrainment, and latent heating. The latter causes temperatures to increase by few 0.1 K in the contrail volume. Supersaturation increases towards the contrail base due to substantially reduced ice crystal number concentrations. The vertical extent of the contrail increases by about 50 m at the base owing to ice crystal settling.
The number of SIPs per grid cell stays close to the initialized value (25), except in the lowest portion (10.1-10.15 km) into which only the largest (diameter D > 10 m) crystals settle. The newly formed contrail volume contains only few ice crystals (number concentration n < 0.1 cm −3 ) that are on average about three times larger than those in the contrail core. The decrease in ice water content, M, in the lower contrail region indicates that vapor uptake is diffusively hindered due to the low ice crystal number concentrations. The increase in M in the upper contrail region occurs because the uppermost ice crystals are exposed to strongly ice-supersaturated, entrained air.
The magenta curves in Figure 3 depict the LO-case. Reducing the deposition coefficient to a more realistic value on the one hand diminishes the rate at which H 2 O molecules are incorporated into the ice phase and on the other hand causes a higher steady-state ice supersaturation that, in turn, increases the vapor flux towards the ice crystals in the contrail core. As a result of these competing effects, mean ice crystal diameter and ice water content slightly increase there. In the contrail top and base regions, where nonequilibrium effects due to lower n prevail, ice crystals remain smaller and have lower fall speeds. Therefore, ice crystals are lifted to higher altitudes at the cloud top, and the fall streak depth is now only 20 m. At the cloud top, M is significantly smaller; together with smaller D (and possibly changes in ice crystal habits not simulated here), this alters shortwave optical depth and radiative flux changes due to the contrail. Owing to their potential importance for contrail radiative forcing, we suggest to study the magnitude of these changes in future work. Figure 4 shows the ambient ice saturation ratio profiles after 1.5 hr. While the contrail region (initially ice saturated, see Figure 2) stays close to saturation due to high n in case HI-, the contrail boundaries develop distinct ice-supersaturated states. High S develops in a thin layer at the contrail top above the altitude that is just void of contrail ice crystals. The contrail base also develops supersaturation, since the number of ice crystals settling into this is too low to deplete vapor in excess of ice saturation faster than it is generated. Entrainment has only little impact on the S-profile in this scenario. Slower H 2 O uptake in case LO-leads to a larger quasi steady-state supersaturation in the contrail correlating well with the decrease in n towards the contrail base (see Figure 3). It is interesting that turbulent diffusion is capable of creating a contrail layer with a base more strongly ice supersaturated than the top. Figure 5 shows contrail-averaged ice crystal size distributions after 1.5 hr. Turbulent diffusion has a large effect; the dispersion of supersaturation driven by diffusion at the contrail boundaries hinders ice crystal growth effectively. Ignoring entrainment of environmental air shifts the entire size distribution to smaller sizes because less water vapor is available for ice growth in the contrail. We note that the size distribution is more sensitive to the choice of the entrainment-dilution parameter, b, than associated changes in the S(z) profile indicate. Ice crystals in the contrail core grow larger in case LO-. As discussed above, the reduced H 2 O Figure 4. Vertical profiles of the ambient ice saturation ratio (black curve) in the contrail after 1.5 hr with a deposition coefficient set to unity. The effects of including turbulent vertical diffusion (blue), using a diffusivity of  = 0.3 m 2 s −1 suitable for young contrails, and ignoring entrainment (turquoise) are also shown. The solid magenta curve was obtained by reducing the deposition coefficient 0.006. uptake rate due to the low -value is offset by a larger steady-state supersaturation. At the same time, kinetic effects reduce the amount of large ice crystals in the fall streak. The size distributions are in qualitative agreement with in situ measurements (Schröder et al., 2000), as far as a direct comparison is possible.
We comment on the limits of modeling the contrail with a constant deposition coefficient. When using a more realistic model for that takes into account its supersaturation dependence, the higher supersaturation simulated in case LO-would increase and thereby enhance growth rates. We therefore expect that mean crystal diameters and levels of supersaturation differ from those discussed here when using a variable approach.
It is well known that persistent contrails are frequently present in situations where natural cirrus can also form or exist, in particular in regions of slow upper air ascent, as simulated here, ahead of a surface warm front (Kästner et al., 1999). With ongoing cooling (in our case beyond 1.5 hr), we envisage a situation where (i) new ice crystal formation commences by homogeneous freezing at the cloud top, whereby the magnitude of ice supersaturation, hence the timing of those freezing events and the nucleated ice numbers, is affected by turbulent diffusion; (ii) the ice crystal number in the contrail core decreases further by dilution allowing supersaturation to build up so that homogeneous ice formation may also take place there; (iii) smaller crystals with lower fall speeds stay in the contrail core and grow only slowly as long as supersaturation is maintained and ice crystal number concentrations stay high; and (iv) large ice crystals continue to settle far enough to sublimate, but smaller ice crystals remain in the contrail core as long as their associated sedimentation time scales into a dry layer below are longer than the contrail lifetime (as determined by ambient ice supersaturation). Exactly to which extent either process occurs depends on details of the initial and environmental relative humidity profile, the temporal development of the vertical wind forcing, and the presence of INP. Taken together, our findings emphasize the continuum nature of persistent contrails regarding their development into contrail cirrus and challenges our ability to distinguish them from co-existing natural cirrus.
Large-eddy simulations of the transition of contrails into contrail cirrus have been performed, either with hypothesized nucleation mechanisms under synoptic forcing conditions (Unterstrasser & Gierens, 2010) or without consideration of natural cirrus formation (Paoli et al., 2017). These earlier studies used conventional uptake models with fixed deposition coefficients. The present results suggest that contrails soon after formation in clear air are very likely found to be embedded in a cirrus cloud that formed in air surrounding 10.1029/2019JD031847 the contrail. In such circumstances, as well as in situations where contrails form within already existing cirrus (Tesche et al., 2016), it will be difficult to measure microphysical and optical properties of the embedded contrails in isolation. This is an issue for (i) measurements with fast-flying aircraft once contrails spread significantly or lose their linear shape due to the limited temporal resolution of cloud particle probes and (ii) remote-sensing measurements regardless of contrail geometry due to the sparse vertical resolution of passive sensors.

Homogeneous Freezing (Cirrus) Simulations
Effects of (i) ice crystal surface kinetics on homogeneous freezing and depositional growth (Zhang & Harrington, 2015) and (ii) temperature fluctuations induced by mesoscale gravity waves on homogeneous freezing  have been studied with parcel models. These studies indicate that (i) surface kinetic limitations do not significantly impact the number of homogeneously nucleated ice crystals owing to the high ambient supersaturation in which homogeneous freezing takes place and (ii) gravity waves significantly enhance the mean, homogeneously nucleated ice numbers by exposing air to a wide range of cooling rates in excess of synoptic values. Furthermore, parcel and 1-D simulations have been performed and compared to investigate the effect of sedimentation on homogeneously nucleated ice crystal number concentrations (Murphy, 2014). Below, we explore effects of sedimentation, vertical turbulent diffusion, and adaptive deposition coefficients on homogeneous freezing events.

Model Setup
To simulate homogeneous ice formation in cirrus, we use the idealized meteorological setup (w = 5 cm s −1 ) and vertical grid underlying Figure 2 to define a reference scenario with ice crystal sedimentation but without H 2 O diffusion. We set the minimum ice crystal number concentration resolved in this homogeneous freezing event to n min = 1 L −1 . In all scenarios, the initial aerosol size distribution is characterized by a total aqueous particle concentration of 500 cm −3 ; the underlying log-normal core size spectrum has a modal mean radius of 20 nm and a geometric standard deviation of 1.5. The solubility parameter is set to 0.5.
The base case (FIX-) assumes a fixed deposition coefficient. We prescribe = 0.7 according to aerosol/ cloud chamber studies of homogeneous freezing and subsequent growth (Skrotzki et al., 2013). Using constant values ≪ 1 during homogeneous freezing events is inconsistent with cloud physical theory (Lamb & Verlinde, 2011), field observations (Kärcher & Ström, 2003), and parcel model calculations (Gierens et al., 2003). We present one simulation that additionally includes the effect of a high background turbulent vertical diffusivity, setting  = 0.1 m 2 s −1 and another that ignores ice crystal sedimentation by setting V = 0 for all SIPs. To compare the outcome of these simulations in terms of nucleated ice number concentrations, we evaluate column-averaged ice crystal number mixing ratios, .
A low growth efficiency impacts depositional growth during and especially after nucleation (Zhang & Harrington, 2015), which is not captured with constant in the above cases. Therefore, we include a case VAR-using a variable deposition coefficient based on ledge nucleation ( = 10) as the primary growth mode. While we recall from section 2.1 that the current understanding of cirrus ice crystal growth from the vapor is limited, in this simulation, more realistically responds to changing ice crystal size and environmental conditions. We use variable critical supersaturations, s cr (T), taken from DiSKICE (Harrington et al., 2019) to predict ; in homogeneous freezing conditions, values s cr ≈ 0.2 indicate relatively inefficient H 2 O uptake.

Results and Discussion
We discuss by means of Figure 6 vertical profiles taken about 20 min after the first SIPs were generated in simulation VAR-. The profiles of n and NSIP match well, indicating the presence of a primary and a secondary nucleation region with local maxima at 10.25 and 10.35 km, respectively. The first ice crystals in the primary region form at the point where homogeneous freezing conditions are met in the prescribed moist layer. The nucleated ice crystals grow rapidly, quenching the high supersaturation and thereby shutting off freezing quickly, while the entire region is advected upwards. Growing ice crystals settle out of the nucleation region as soon as their fall speeds offset the updraft. This freezing event lasts only a few minutes. The secondary nucleation region is located in the upper branch of S(z), marking the cloud top. Contrary to the first freezing event, ice crystals are created continuously at this dynamic upper cloud boundary while they sediment towards the first zone. Together, this explains the occurrence of a broad ice crystal mode with mean diameter D = 30-50 m connecting the nucleation zones. Towards the cirrus base, ice crystals almost double their size as they pick up vapor while falling through unperturbed supersaturated air. If is initially high as in the case of fresh, homogeneously nucleated ice crystals, rapid deposition growth reduces ambient s and therefore surface supersaturation. This reduces , and the diminishing growth rates allow a new vapor gradient to build up increasing ambient s, in turn increasing . The deposition coefficient is highly variable throughout the ice-supersaturated growth zones. A quasi steady-state supersaturation builds up due to the nonlinear feedback between surface kinetics of H 2 O and ice crystal growth. Ice crystals begin to sublimate below 10.1 km, where transitions from values basically 0 around ice saturation to 1 (not shown) in ice-subsaturated conditions. These results illustrate the complexity that occurs solely due to the nonuniformity in nucleation and growth conditions typically present in a column of air. Figure 7 shows column-averaged ice crystal number mixing ratios as a function of time. In the base case FIX-, we see the sharp onset of the primary homogeneous freezing event right after t = 75 min, maximizing shortly thereafter. While in a parcel model, the ice crystal mixing ratio would stay constant after reaching its peak value, decreases in the 1-D simulation, and its evolution is largely determined by the depth of the nucleation and growth zones. The initial decrease in is caused by the nonuniformity of freezing and growth conditions discussed above. In addition, sedimentation spreads out ice crystal concentrations and reduces the peak ice crystal concentration in the first freezing event. Ice crystals formed in the secondary freezing event are too few to prevent from decreasing in this case. Due to the smoothing effect an enhanced turbulent diffusivity has on S, vertical diffusion shifts the onset of ice formation by ≈ 4 min and reduces the peak mixing ratio. This reduction of is probably overemphasized, as the prescribed diffusivity is at the high end in background conditions. By contrast, a supersaturation-dependent deposition coefficient from case VAR-generates more crystals relative to case FIX-since declining values retard the quenching of supersaturation and allow more ice crystals to nucleate. This process is very sensitive to the numerical time step due to the strong dependence of the freezing rate coefficient on water activity.
Taken together, peak ice number concentrations either using = 0.7 or variable overestimate column-averaged 1-D model predictions, depending on the time after first freezing. The inclusion of surface kinetics increases simulated peak ice number concentrations relative to a simulation with constant . It seems that the effect is small and that parameterizations based on parcel models with fixed might be approximately correct in predicting homogeneously nucleated ice numbers. However, this conjecture is Figure 7. Column-averaged mixing ratio of homogeneously nucleated ice crystal numbers versus time assuming (magenta curve) a constant deposition coefficient. Effects of (orange) ignoring sedimentation and (blue) vertical diffusion are also shown for the fixed case. The black curve was obtained by using variable deposition coefficients including sedimentation and neglecting diffusion. based on only one case using a single, low updraft speed, ignoring wave effects and employing the underlying initial conditions shown in Figure 2. Systematic studies are needed to quantify homogeneously nucleated ice numbers along with differences between parcel and 1-D model simulations. Effects of sedimentation, diffusion, and nonuniform nucleation and growth conditions on ice crystal number concentrations are interdependent and will be difficult to disentangle in observations.

Summary and Outlook
The Lagrangian trajectory column model primeice1D enables highly detailed process studies of aerosol-cirrus interactions. Like other 1-D cirrus models, it treats vertical advection, ice crystal nucleation by homogeneous aerosol freezing, and ice crystal growth from the vapor, sublimation, and sedimentation. It shares the tracking of individual SIPs and size-and composition-dependent (nonequilibrium) aqueous aerosol particle growth and evaporation with other column models. As novel features, primeice1D includes turbulent diffusion of water vapor and heat as well as spheroidal vapor fluxes combined with deposition coefficients that are constrained by laboratory measurements and respond to ice crystal size and environmental conditions. Once measurements become available to constrain or advance theories of faceted growth, the model is well poised to predict primary ice crystal habits. Besides describing the model physics, we demonstrate its potential for further in-depth studies with the help of simulations resolving meter scales.
We find that persistent contrails forming in clear air regions of sustained, slow upper air ascent develop into contrail cirrus within a few hours, either by homogeneous freezing of aqueous aerosols or in a more complex manner if other ice-forming particles are also present. As the associated meteorological situation signifies preferred contrail formation regions and aqueous aerosols are abundant, this type of transition, which embeds contrails in natural cirrus clouds, is likely to happen rather frequently. Our results further suggest that vertical diffusion and entrainment need to be considered to accurately simulate contrail ice crystal size distributions. Realistically low (but fixed) deposition coefficients enhance the steady-state supersaturation inside the contrail and increase the mean ice crystal size. More studies of young contrails are needed to address the impact of variable deposition coefficients.
We corroborate earlier findings that the nonuniformity in vertical profiles of moisture and temperature conditions gives rise to differences in homogeneously nucleated cirrus ice crystal number concentrations when compared to spatially uniform parcel simulations. In conjunction with sedimentation, this 1-D effect lowers nucleated ice numbers over the first minutes past freezing. We additionally find that in conditions of sufficiently strong turbulent diffusion, ice formation onset is delayed and peak nucleated ice numbers are reduced. Variable deposition coefficients enhance the quasi steady-state, in-cloud supersaturation, decrease mean ice crystal sizes, and increase homogeneously nucleated number concentrations. More studies are needed to fully map out the conditions in which such offsetting effects make ice numbers comparable to results obtained with nucleation parameterizations based on conventional parcel simulations using constant (high) deposition coefficients.
Our results suggest that variable growth efficiencies broaden the distribution of relative humidity in cirrus clouds. We find that supersaturation-dependent deposition coefficients can have opposite effects on water vapor mass uptake in ice clouds, depending on the number concentration of ice crystals. For high ice numbers, such as in young contrails where a low quasi steady-state supersaturation is established, variable deposition coefficients increase that level of supersaturation and hence ice crystal sizes. For comparatively low ice numbers, such as in natural cirrus, variable deposition coefficients decrease the rate of water vapor uptake in the declining supersaturation right after homogeneous freezing and thereby reduce ice crystal sizes. Exactly which effect dominates is influenced by the prevailing cooling rates.
The nonlinear interaction between gas phase and ice crystal surface processes that arises from the use of realistic kinetic models for water vapor uptake will also control the relative importance of various INPs in cloud formation for a given dynamical forcing. Capturing this issue is particularly important in simulations with multiple ice-forming aerosol particle types, as each particle will experience different growth conditions and trigger different supersaturation relaxation times during a nucleation event. Arguably, this strongly influences the total number of nucleated ice crystals present at the beginning of the cloud's life cycle. It will be interesting to see how this evolves in cloud-or cloud system-resolving models. Taken together, we expect that supersaturation-dependent deposition coefficients will be an important component in future simulations of aerosol-cirrus interactions, besides the effect of mesoscale gravity waves on the vertical wind speeds driving ice cloud microphysics. The prospect of predicting primary ice crystal habits will improve estimations of cirrus-induced shortwave radiative flux changes.