Evaluating an Ice Crystal Trajectory Growth (ICTG) Model on a Quasi‐Idealized Simulation of a Squall Line

Abstract

A major challenge in numerical weather prediction models is the ability to accurately simulate the microphysical properties and growth of ice hydrometeors in clouds. Eulerian bulk microphysics schemes in these models tend to obscure the properties and evolution of individual ice crystals, often resulting in inaccurate simulations of storm structures. To address this issue, this study presents a novel ice crystal trajectory growth (ICTG) model that simultaneously grows and advects individual ice crystals while tracking their evolving properties along their trajectories. The model is evaluated on a 3D quasi‐idealized leading‐convective, trailing‐stratiform squall line simulation. The ICTG model successfully produced a spatial distribution of ice crystal trajectories consistent with the simulated reflectivity structure of the storm above the melting level. Smaller initialized crystals (d ≤ 0.1 mm) were largely transported to the anvil and the trailing stratiform region. One primary trajectory involved sustained growth in the stratiform mesoscale updraft for ∼1.5 hr, resulting in a density reduction down to 600 kg m⁻³, a final particle size greater than 0.9 mm, and potential branching. In contrast, larger initialized crystals (d ≥ 0.5 mm) collected more rime and fell out primarily in the leading convective line. The ICTG model's realistic production of varied crystal growth properties owing to differences in transport and initial size suggests that it can be a valuable tool for learning about ice microphysical processes in a variety of cold cloud systems.

Introduction

Ice crystals and their growth processes (e.g., vapor deposition, riming) in mixed-phase clouds have been a longstanding source of uncertainty in the scientific understanding of organized cloud systems. Ice crystal growth processes can have a large impact on cloud structures, from modifying the cloud top radiative cooling to modifying the in-cloud vertical motions during phase change. The growth history of ice crystals can greatly determine which regions of a cloud the crystals are transported to, which can further influence the cloud's spatial extent and lifetime. Understanding these interactions is especially important for improving simulations and forecasts of organized precipitating clouds.

The leading-convective trailing-stratiform squall line (Parker & Johnson, 2000) is a classic example of a storm system whose cloud and precipitation structure is largely modulated by ice crystals. Its kinematic structure can be modified through various mechanisms, which include freezing droplets in the leading line strengthening convective updrafts (Fovell & Ogura, 1988; Tao & Simpson, 1989), sublimating crystals beneath the stratiform anvil strengthening the mesoscale downdraft (Rutledge et al., 1988; Yang & Houze, 1995), and melting crystals enhancing the negative buoyancy of the cold pool (Fovell & Ogura, 1988; Leary & Houze, 1979). Ice generated in the leading convective line is transported rearward to compose a large fraction of the stratiform precipitation, with both regions separated by a transition zone characterized by a relative minimum in radar reflectivity (Biggerstaff & Houze, 1991, 1993; Braun & Houze, 1994; Fovell & Ogura, 1988; Rutledge & Houze, 1987; Smull & Houze, 1985; Yang & Houze, 1995). There are several hypotheses on the mechanisms generating this reflectivity trough: ice crystal sublimation (Biggerstaff & Houze, 1991), lack of aggregation (Biggerstaff & Houze, 1993), suppressed crystal growth by downdrafts (Braun & Houze, 1994), and crystal size sorting in which the heavier particles fall out in the convective region and the lighter, smaller particles are transported rearward (A. A. Jensen et al., 2018; Braun & Houze, 1994; Fovell & Ogura, 1988; Rutledge & Houze, 1987). The fallout pattern of ice particles also impacts the intensity and across-line extent of the trailing stratiform precipitation (Rutledge & Houze, 1987; Smull & Houze, 1985), and without the rearward transport of ice from the convective cells, little to no precipitation would reach the surface in the stratiform region (Chin, 1994; Fovell & Ogura, 1988; Rutledge & Houze, 1987; Tao & Simpson, 1989). The varied results of these and other studies demonstrate that the structure of squall lines in numerical model simulations is sensitive to the microphysics parameterizations used (Morrison et al., 2009, 2012, 2015; Bryan & Morrison, 2012; Fovell & Ogura, 1988; Tao & Simpson, 1989; Yang & Houze, 1995). This sensitivity reveals the need to better determine fundamental controls on squall line precipitation features for their accurate representation by numerical models.

In numerical model simulations, most microphysics schemes represent small-scale ice processes using an Eulerian method to track the collective properties of a group of particles advecting through each grid box of the domain. Specifically, bulk microphysics schemes make implicit and explicit assumptions about how ice crystal properties (e.g., shape, density, fall speed) are distributed with respect to some measure of size. This Eulerian framework effectively averages these properties so that individual particle characteristics are lost, often resulting in erroneous, varied simulations of a given storm (e.g., Morrison et al., 2015). Numerical models are only now beginning to include microphysical methods that track the collective changes to the physical properties of ice crystals. These particle property schemes have been developed for both bulk (i.e., Jensen et al., 2017; Morrison & Milbrandt, 2015) and bin (Harrington et al., 2021) applications. However, the lifecycle of an ice crystal's properties is known to depend on the specific paths traced by the crystal through a cloudy atmosphere. Nelson (2008) has estimated that the great diversity in observed snow crystals is due to the heterogeneity of the trajectories traced by the crystals in clouds, suggesting that accurate modeling of ice populations requires Lagrangian modeling.

Lagrangian trajectory modeling provides a way to examine the diversity and evolution of ice particle properties by tracking the specific growth pathways taken by individual crystals, which is not possible in Eulerian frameworks. Lagrangian modeling has not yet been used on particle property schemes even though Lagrangian methods are perfectly suited to the task. Modeling the simultaneous transport and growth of individual ice crystals is challenging and computationally expensive, though recent “super-particle” Lagrangian methods are promising (e.g., Grabowski et al., 2019; Shima et al., 2020). The advantage of a Lagrangian approach is a more explicit representation of ice crystal growth processes, the accompanying crystal structural changes, and their interaction with the larger-scale cloud structure, which are all needed for studying cloud systems and improving bulk microphysics schemes.

Lagrangian particle trajectory studies have been used to examine the growth and transport of hydrometeors in a variety of cloud types. Numerous studies have conducted trajectory calculations on hailstones to understand storm controls on hailstone size and fallout location (e.g., Browning & Foote, 1976; Heymsfield, 1983; Kumjian & Lombardo, 2020; Miller et al., 1990; Nelson, 1983; Rasmussen & Heymsfield, 1987; Tessendorf et al., 2005). Similarly, liquid droplet trajectories have been analyzed in tropical cyclones (Xu et al., 2017), cloud mergers (Kogan & Shapiro, 1996), maritime clouds (Szumowski et al., 1999), arctic clouds (Harrington et al., 2000), and marine stratocumulus (Hartman & Harrington, 2005). Lackmann and Thompson (2019) computed snowflake trajectories in mesoscale snowbands to understand the mechanisms determining horizontal distributions of snowfall, and E. J. Jensen et al. (2018) computed trajectories of ice crystals detrained from deep tropical convection into the anvil. In comparison, there is a relative lack of ice crystal trajectory studies on squall lines. Past studies utilized Doppler radar observations to estimate ice particle trajectories in squall lines (Biggerstaff & Houze, 1991, 1993; Braun & Houze, 1994), but they lacked detailed information about the evolution of particle properties. Yang and Houze. (1995) briefly examined ice particle trajectories in a 2D model simulation of a squall line, but the particle growth trajectories were not the focus of their study.

The current study presents a novel ice crystal trajectory growth (ICTG) model that is used to examine ice processes and their impacts on squall line structure. Our ICTG model simultaneously grows and advects ice crystals in three-dimensional space, tracks physical changes of individual particles, and incorporates recent laboratory-based parameterizations of ice-vapor deposition (Harrington et al., 2019; Zhang & Harrington, 2014) and riming (Jensen & Harrington, 2015). The goals of this study are: (a) To introduce the new ICTG model, (b) To test the model's ability to accurately capture the spatial distribution and properties of ice particles in a simulated squall line in the context of past studies, and (c) To gain insight on how variability in an ice particle's initial size and location impacts the particle's trajectory, its evolution, and the squall line's precipitation structure. The remainder of this paper is organized as follows: Section 2 presents a description of the ICTG model, Section 3 describes how we employ the model to study squall lines, Section 4 and 5 show results from 2D and 3D ICTG simulations, and the conclusions are in Section 6.

Ice Crystal Trajectory Growth Model Description Model Overview

The ICTG model in the current study was developed to capture the natural evolution of ice crystal properties throughout their growth life cycle and learn how the crystal's evolution contributes to the microphysical structure of organized cloud systems. The model simultaneously grows and advects ice crystals by using the thermodynamic and kinematic output from 3D numerical model simulations. The model operates in a Lagrangian framework by tracking the steady-state thermodynamic conditions that each crystal experiences and calculating the resultant evolving particle properties (e.g., axis lengths, effective density, fall speed). The ICTG model runs completely offline and therefore does not produce any dynamic or thermodynamic feedbacks with the 3D simulation. All of the particles captured by the ICTG model are in the ice phase only, so from here on, the words “particle” and “crystal” will be used interchangeably.

The ICTG model is initialized with user-specified particle sizes, aspect ratios, densities, and locations in the domain. Then, over small, user-specified time steps (<20 s; see Harrington et al. (2013)), each particle's evolution from vapor deposition, sublimation, and riming of cloud droplets is calculated. The model calculates each ice particle's evolution individually, and it does not include melting, aggregation, or secondary ice production. Although aggregation is an important ice growth process in stratiform clouds, in this study we only use the ICTG model for proof-of-concept simulations to test the sensitivity of the results to changes in the crystal's initial size (see Section 3 for a description of how we employ the ICTG model in this study). After the particle's growth is calculated, the fall speed is updated, and the particle is advected some distance for each time step. The trajectory model uses first-order upwind advection, in which the velocity vector assigned to each particle comes from the nearest grid point at each time step. The particle's horizontal velocity is assumed to equal that of the horizontal wind, and the particle's vertical velocity is the difference between the vertical wind and the fall speed. Each particle is advected until it reaches the melting level or until its mass completely sublimates.

Vapor Deposition

Ice crystals can grow by vapor deposition to precipitation-sized particles with shapes determined by surface attachment kinetics that depend on the temperature and ice supersaturation. At temperatures above −20°C, crystal shapes generally fall into one of two primary habits: planar or columnar. Planar crystals usually have their major dimension parallel to the basal (hexagonal) facet, whereas columnar crystals have their major dimension parallel to the prism (rectangular) facet. The classical hexagonal and prism structure of single crystalline ice suggests describing crystal size in terms of at least two dimensions referenced to the basal and prism facets. We define these dimensions as c, which is half the maximum distance along the prism facet, and a, which is half the maximum distance along the basal facet. Crystal shape is represented by oblate (planar) and prolate (columnar) spheroids, which allows us to describe crystal shapes with two dimensions instead of one. The crystal's aspect ratio, defined as ϕ = c/a, is used to distinguish the habit, where 0 < ϕ < 1 indicates an planar crystal, ϕ > 1 indicates a columnar crystal, and ϕ = 1 indicates an isometric crystal. The rate at which each facet grows determines the crystal's habit, and the growth rate of a facet depends on the gas-phase vapor and thermal fluxes in combination with surface attachment processes. To refer to the crystal size, we use an equivalent volume spherical diameter, defined as $$d=2{\left({a}^{2}c\right)}^{1/3}$$ (Harrington et al., 2013).

Vapor Deposition Mass Growth Rate

In order for ice crystals to grow, vapor molecules must first adsorb onto the crystal surface and then find an attachment site. Steps in the crystal surface provide the attachment sites, and at least two mechanisms are known to govern the density of steps. Permanent screw dislocations are a continual source of surface steps, providing numerous attachment sites for adsorbed vapor molecules. In the absence of dislocations, two-dimensional nucleation produces surface steps and this mechanism has a strong supersaturation dependence. Since the surface mechanisms that control step formation and propagation are not well understood, surface attachment processes are usually modeled with growth efficiencies. The efficiency with which vapor molecules can incorporate into the ice lattice is given by a deposition coefficient α, where a value of unity means perfect incorporation. We use the parametric form of α proposed by Nelson and Baker (1996): [Image Omitted. See PDF]

Where s_surf is the surface supersaturation, s_char is the characteristic supersaturation, and m is an adjustable parameter that determines the crystal growth mechanism. Growth by permanent dislocations occurs for m = 1 (Burton et al., 1951), and growth by 2D nucleation occurs for 10 ≤ m ≤ 30 (Nelson & Baker, 1996). Ice crystal growth is not largely a function of m once m > 10, though notable differences in the crystal semi-axis length (∼800 μm) can occur at −15°C when growing crystals for 15 min via permanent dislocations versus 2D nucleation (Harrington et al., 2019; Zhang & Harrington, 2014). For the simulations in this study, we follow Zhang and Harrington (2015) and set m to 15, corresponding to growth by 2D nucleation, which is necessary to realistically produce thin ice crystals such as dendrites (Frank, 1982; Harrington et al., 2019; Nelson & Baker, 1996).

The capacitance model has traditionally been used to calculate the ice crystal mass growth rate from vapor deposition. Although the capacitance model assumes a rough crystal surface and perfectly efficient growth, to better approximate the growth of real crystals, we instead use the model of Zhang and Harrington (2014) that combines the capacitance model with a model for surface attachment through steps. Zhang and Harrington (2014)'s model approximates faceted crystals as spheroids and includes axis-dependent attachment kinetics. In their model, the mass growth rate of an ice crystal from vapor deposition is [Image Omitted. See PDF]

where C is the capacitance, s_i is the ambient ice supersaturation, ρ_eq is the equilibrium vapor density at the ambient temperature T, l_s is the enthalpy of sublimation, R_v is the gas constant for water vapor. The variables $${D}_{v}^{\prime }$$ and $${K}_{T}^{\prime }$$ include impacts from attachment kinetics on each facet, and they take the following forms: [Image Omitted. See PDF] [Image Omitted. See PDF]

where C_Δ is the capacitance a small distance Δ above the crystal's surface, D_v is the vapor diffusivity of air, K_T is the thermal conductivity of air, f_v and f_T are the ventilation coefficients for vapor diffusion and thermal conduction, c_p is the specific heat capacity of air at constant pressure, ρ_a is the air density, and $${\bar{v}}_{v}$$ and $${\bar{v}}_{a}$$ are the mean molecular speeds of water vapor and air molecules, respectively. Note that Equation 3 also contains deposition coefficients for each axis of the crystal, α_a and α_c, which control the efficiency with which each facet grows. In Equation 4, α_T,a and α_T,c are the thermal accommodation coefficients for each axis. We set both α_T,a and α_T,c to unity in the ICTG model, following Shaw and Lamb (1999).

Aspect Ratio and Density Evolution From Vapor Deposition

The aspect ratio of an ice crystal evolves through differences in α_a and α_c, which have the functional form of Equation 1. Differences between their magnitudes come from s_char and s_surf, which is calculated above each axis (see Zhang and Harrington (2014), pg. 378). Aspect ratio evolution follows the hypothesis from Nelson and Baker (1996), which states that the ratio of the axis growth rates is equal to the ratio of the deposition coefficients for each axis, Γ(T): [Image Omitted. See PDF]

We use this hypothesis because it is consistent with growth of faceted ice with steps nucleated at the crystal edges (Harrington et al., 2019).

At high ice supersaturations, secondary ice crystal habits such as dendrites, hollowed features, and rosettes can develop. Such secondary habits are not explicitly represented in microphysical models. These features are represented through a deposition density (ρ_dep) that accounts for the branches and hollows that occur in real crystals (Chen & Lamb, 1994): [Image Omitted. See PDF]

Here, ρ_dep depends only on the temperature and the vapor density excess above the equilibrium value (Δρ_v = ρ_v − ρ_eq). We add the deposition density to the particle effective density by using volume-weighting: [Image Omitted. See PDF]

Where ρ_i,1 (ρ_i,2) and V_i,1 (V_i,2) are the ice crystal effective density and volume before (after) one time step. During periods of sublimation, the aspect ratio is assumed to remain constant, which has been documented in laboratory measurements (Nelson, 1998).

Ventilation

Once an ice crystal becomes large enough to appreciably fall, it can amplify the thermal and vapor fluxes around it, which can enhance the particle's growth rate through ventilation. In the ICTG model, we assume that ventilation only occurs as a result of the crystal's vertical falling motion and neglect ventilation from the crystal's horizontal motion. From Hall and Pruppacher (1976), the ventilation coefficient for vapor diffusion is given by [Image Omitted. See PDF]

where a₁, a₂, and b are constants, N_Re is the Reynolds number of the ice crystal, and N_Sch is the Schmidt number of water vapor. Reynolds numbers are calculated using the Reynolds number-Best number power law relationship from Mitchell (1996). We use power law coefficients with dependencies on a particle's projected cross-sectional area and mass-diameter relation such that the coefficients are general for all ice types (Mitchell & Heymsfield, 2005). The ventilation coefficient for thermal diffusion follows a similar functional form: [Image Omitted. See PDF]

where a₁, a₂, and b are constants, and N_Pr is the Prandtl number of air.

Riming

Riming is the process of an ice crystal collecting supercooled liquid droplets that freeze onto the crystal surface, which can drastically increase the crystal's mass and fall speed while only mildly increasing its maximum dimension. Riming in the ICTG model is only from collection of cloud droplets and is based largely on Jensen and Harrington (2015), who developed a method to capture the evolving aspect ratio of ice crystals during riming. The following is a brief summary of their riming model, with slight adjustments made in the ICTG model.

Riming Mass Growth Rate

The mass growth rate of an ice crystal from riming is given by [Image Omitted. See PDF]

where E_il is the efficiency with which the ice crystal collects cloud droplets of radius r_l, A_g is the combined geometric cross-sectional area of the crystal and droplet, v_i is the fall speed of the ice crystal, and v_l, m_l, and n_l are the fall speed, mass, and number concentration of a cloud droplet of radius r_l, respectively. In Equation 10, crystals are assumed to fall with their major axis perpendicular to their fall direction. The cloud droplets are assumed to follow a log-normal distribution with 200 size bins, and the distribution is calculated from the cloud droplet mass and number mixing ratios output from the 3D numerical model simulation. Because the fall speeds of cloud droplets are generally much smaller than those of precipitating ice crystals, the cloud droplet fall speed in Equation 10 can be neglected. Depending on the crystal habit, different collision efficiencies are calculated. Spherical collision efficiencies are from Beard and Grover (1974), plate efficiencies are from Hall (1980), and columnar efficiencies are based on a modified spherical efficiency developed in Jensen and Harrington (2015).

Aspect Ratio Evolution From Riming

Heymsfield (1982) found that riming acts to primarily increase the minor axis of an ice crystal with little changes to the major axis length. So as ice crystals become heavily rimed, they lose their initial shape. Once their underlying crystal shape is no longer discernible, they are considered to be graupel, which has an aspect ratio of approximately 0.8 (Böhm, 1992; Heymsfield, 1978). Given these factors, Jensen and Harrington (2015) only add rime to the minor axis and allow crystals to rime with an evolving aspect ratio until it reaches a limiting value, with plates riming to an aspect ratio of 0.8 and columns to an aspect ratio of 1.25 (the inverse of 0.8). Once these limiting values are reached, the aspect ratio remains constant thereafter, provided that growth only occurs by riming. Note that in the ICTG model, the ice crystal aspect ratios can evolve freely after these limiting values are reached if the crystals are advected into regions with no liquid cloud water but continue to undergo vapor growth.

In the ICTG model, we compute the density of the rime mass following Macklin (1962): [Image Omitted. See PDF]

where r_i is the radius of the ice crystal. The rime density is added to the full crystal density using volume-weighting. Limits are imposed on the full crystal density (whether it is undergoing riming and vapor growth or vapor growth only in a given time step), with a lower limit of 50 kg m⁻³ and an upper limit of 917 kg m⁻³. The lower limit has been used to approximate the density of thin dendrites (Ovchinnikov et al., 2014).

Fall Speeds

Ice crystals in the ICTG model are initialized as stationary with 0 m s⁻¹ fall speed. After the first time step, ice crystals fall at terminal velocity. Terminal velocities are computed following Mitchell (1996): [Image Omitted. See PDF]

Here, η_a is the dynamic viscosity of air, ρ_a is the air density, and L is a characteristic length of the ice crystal. Böhm (1992) determined that $$L=2\sqrt{ac}$$ for columns, L = 2a for plates, and L = 2r for spheres. This characteristic length scale allows the fall speed to be habit-dependent and assumes the crystals fall with their major axis perpendicular to their fall direction.

Simulation Design 3D Quasi-Idealized Squall Line Simulation Setup

To test the performance of the ICTG model, this study uses output from the ARW configuration of the Weather Research and Forecasting (WRF) Model, version 4.1.2 (Skamarock et al., 2019), which is a non-hydrostatic, fully compressible model. We use WRF to simulate a three-dimensional quasi-idealized squall line based on a leading-convective, trailing-stratiform squall line that tracked through central Oklahoma on 19 June 2007. This particular case has also been used to examine electrification and test microphysics and dynamics parameterizations in squall lines in previous studies (A. A. Jensen et al., 2018; Lang et al., 2010; Lebo & Morrison, 2014; Morrison et al., 2012, 2015).

We use the same simulation setup as Morrison et al. (2015), except that we use the aerosol-aware microphysics scheme from Thompson and Eidhammer (2014). Radiation, boundary layer physics, and surface fluxes are turned off for simplicity. The domain is 612 × 122 km with a horizontal grid spacing of 1 km. The domain has 100 vertical levels with a top at 25 km, and Rayleigh damping is applied to the top 5 km. The bottom and top boundaries are rigid, the zonal (across-line) boundaries are periodic, and the meridional (along-line) boundaries are open.

The squall line is initialized with a combined sounding from Lamont, Oklahoma beneath 700 hPa and Norman, Oklahoma above 700 hPa at 0000 UTC that was smoothed using a 1-2-1 moving average (Morrison et al., 2015). The initial environmental wind shear is 12 m s⁻¹ from 0 to 5 km altitude in the across-line direction and 0 m s⁻¹ in the along-line direction. Convection initiates from forced horizontal convergence at the surface and domain center that decreases with height for the first hour. Additional details on the simulation setup are found in Morrison et al. (2015). The model is run with a 2.5-s time step and reaches a quasi-steady state by hour 6. At hour 6, the simulated squall line is a mature storm with a distinct low-reflectivity transition zone, an upshear-tilted convective updraft, and a broad, mesoscale updraft in the stratiform region (Figure 1). The convective region is defined as the region with a maximum in line-average reflectivity from x = 260–300 km, the transition zone as the region with a relative minimum from x = 240 to 260 km, and the stratiform region as the region with a relative maximum from x = 190 to 240 km.

Figure 1. (a) Plan view of simulated reflectivity at 6 hr and 0.5 km altitude using the Thompson-aerosol microphysics scheme (Thompson & Eidhammer, 2014). Line-average (b) reflectivity, (c) vertical velocity, with positive values indicating upward motion, and (d) across-line velocity, with positive values indicating rear-to-front flow.

We interpolate the vertical dimension of the 3D WRF simulation to Cartesian coordinates with 0.25-km vertical grid spacing, and use output from the simulation at hour 6 for the ICTG model. Therefore, the thermodynamic and kinematic fields from the WRF simulation that are ingested into the ICTG model are static in time. Because the ICTG model does not provide feedbacks to the WRF simulation, the thermodynamics, microphysics, and dynamics from the WRF simulation do not evolve in response to particle growth in the ICTG model.

Initializing the ICTG Model

The ICTG simulations assume crystals are growing from frozen cloud or drizzle drops initiated in the leading convective edge of the squall line. This reasoning is based on the squall line's thermodynamic and precipitation structure, in which considerable rain and cloud water are present above the melting level in the convective line under highly ice-supersaturated conditions (Figures 2a–2d). Additionally, above 6 km, the convective region contains a sharp vertical increase in snow and cloud ice mixing ratios and a decrease in cloud droplet and rain mixing ratios, indicating ongoing freezing of liquid droplets (Figures 2e and 2f). Given the presence of copious ice particles in this region, it is reasonable to assume that supercooled droplets are freezing by immersion freezing and contact nucleation and are subsequently growing into larger ice particles. It is worth noting that this ice nucleation method is used merely as a way to select reasonable initial crystal sizes. Ice crystals will certainly also form through deposition nucleation, but these crystals will be initially small and of similar size to crystals formed from small cloud drops.

Figure 2. The shading indicates line-average (a) supersaturation with respect to ice, (b) temperature, and mass mixing ratio for (c) cloud droplets, (d) rain, (e) cloud ice, (f) snow, and (g) graupel. In all panels, the line-average reflectivity is contoured in black at 10, 20, 30, and 50 dBZ. The cyan contour is the line-average 0°C isotherm, and the magenta square outlines the crystal initialization domain.

We simulate 3D ice crystal trajectories using the full 3D WRF output, and we also simulate 2D ice crystal trajectories using an along-line average of the WRF output across the y-dimension of the domain. To decide where to initialize the ice crystals, we define a set of criteria based on the line-average temperature (Figure 2b) and vertical velocity (Figure 1c). Crystals are initialized above the maximum updraft within the temperature range of −10 to −24°C, which is conducive to ice crystal nucleation and growth from frozen cloud droplets (Auer, 1972). In the across-line direction, crystals are initialized 10 km ahead to 20 km rearward of the line-average maximum updraft. These criteria yield an ice crystal initialization domain with vertical limits of 5.75–7.75 km and across-line limits of 267–297 km. Ice crystals are initialized every 1 km horizontally and every 0.25 km vertically. This choice of initial particle placement results in 30,969 total ice crystal trajectories calculated in a single 3D simulation. For a given 2D simulation, 279 total trajectories are calculated.

Ice crystals are initialized as spheres with a density of 917 kg m⁻³. We focus on four different initial crystal sizes for both the 2D and 3D simulations. These initial sizes were chosen based on the cloud droplet and raindrop size distributions within the crystal initialization domain (Figure 3). These distributions do not vary considerably with height, so the initialization domain-average distribution is used to determine initial crystal sizes. Two initial sizes were chosen to be the 50th percentiles of the domain-average cloud droplet and rain drop size distributions (based on number of drops). These diameters are d = 0.04 and 0.1 mm, respectively, which represent cloud- and drizzle-sized particles. The cloud droplet size distribution was calculated using the non-zero average cloud mass and number mixing ratios in the 3D crystal initialization domain. An inverse exponential size distribution was calculated for rain, which is the distribution type used in the Thompson-aerosol microphysics scheme (Thompson & Eidhammer, 2014). We removed part of the raindrop distribution at diameters less than 0.05 mm as this size range is typically within the range of cloud droplet sizes (Figure 3). The 50th percentile of the raindrop size distribution was then calculated. The diameter that we select to separate the cloud droplets from the precipitation-sized particles is somewhat arbitrary, as there is no specific size that delineates between these particles.

Figure 3. (a) The horizontally averaged log-normal cloud droplet size distribution at every 0.5 km in altitude in the crystal initialization domain. The vertically and horizontally averaged distribution over the full crystal initialization domain is indicated by the dotted black line. (b) As in (a), but for the raindrop size distribution. The vertical black lines in both panels indicate the initial particle sizes used for the trajectory simulations.

To examine additional variability in realistic growth trajectories, larger initial particle sizes were also chosen: d = 0.5 mm (drizzle/raindrop size) and 1 mm (raindrop size). Trajectories from other interim sizes (0.01, 0.07, 0.3, 0.7 mm) were also examined and are briefly discussed, but we determined that the aforementioned four initial sizes best capture the largest variations in particle properties and trajectories. As shown in Table 1, the four main trajectory simulations are named according to their dimensional size (2D or 3D) and their initial particle size (cloud, drizzle, drizzle/rain, rain). Both the 2D and 3D trajectory simulations are integrated with a 15-s time step for a maximum of 6 hr. The simulations are stopped sooner if the particle sublimates to extinction or reaches 4.5 km altitude (just above the melting level) before 6 hr.

Table 1 Simulation Names for the Four 2D and 3D Simulations of Focus in This Study

	d = 0.04 mm	d = 0.1 mm	d = 0.5 mm	d = 1 mm
2D	SIM2D-CL	SIM2D-DRZ	SIM2D-DRZRA	SIM2D-RA
3D	SIM3D-CL	SIM3D-DRZ	SIM3D-DRZRA	SIM3D-RA

2D Ice Crystal Trajectory Simulations Ensemble of Trajectories

To test the ICTG model, we first examine the trajectory patterns from the four main 2D ice crystal trajectory simulations (Table 1), with interim particle sizes added for completeness. Shown in Figure 4, the 2D simulations have distinctly different particle trajectory patterns, which vary according to the initial crystal diameters. The majority of the smaller crystals (d ≤ 0.1 mm) are transported rearward into the trailing stratiform region. These crystals also end up in the leading anvil, but only crystals initially smaller than 0.07 mm are light enough to be transported to the trailing anvil. The particles larger than 0.3 mm are too heavy to be lofted rearward to the stratiform region and instead reach the melting level either in the transition zone or the leading convective line. Only the largest particles (d ≥ 0.7 mm) reach the melting level in the leading convective line, consistent with the greater reflectivity values in this region. In all of the simulations, particles initiated in the highly ice-subsaturated front edge of the convective line quickly sublimate to extinction, especially for the particles initially smaller than 0.1 mm. These collective trajectory patterns provide an initial indication that the ICTG model is performing as intended.

Figure 4. Ice crystal trajectories from the four 2D simulations of focus (Table 1) in the right column, with trajectories from additional interim sizes in the left column. The crystal trajectory colors vary by initial across-line position; initial positions and color coding of crystals for all simulations is shown in the inset of (h). Black (white) circles at the end of trajectories mark the final position of particles that did not sublimate (did sublimate) to extinction. The gray square outlines the crystal initialization domain. Shading represents the line-average supersaturation with respect to ice. Black contours are line-average reflectivity at 10, 20, 30 and 50 dBZ. Magenta contours are line-average vertical velocity at 0.5, 1, 2, and 4 m s−1. The lower height cutoff in each panel is 3 km.

From here on, our discussion will focus on results from the four main particle sizes as these sizes best capture the largest distinctions in the trajectory patterns and resultant particle properties. Figure 5 shows the trajectory density (number of times a particle passes over a grid box non-consecutively) for the four main 2D simulations, indicating the same trajectory patterns but more clearly showing the highly trafficked regions and the spatial spread in trajectories. In SIM2D-CL (Figure 4b), the particle trajectories are organized according to their initial position relative to the convective updraft. The particles initiated rearward of the updraft (dark blue trajectories in Figure 4b) between 7.25 and 7.75 km altitude are most likely to be transported into the forward and rear anvils. Many of these particles have not yet reached the melting level after 6 hr of ICTG model integration, consistent with the trajectory retrieval in (Braun & Houze, 1994), which found that particles in the trailing anvil had the longest trajectory times. Their slower fall speeds suggest that they have fall speeds nearly balanced by the weak mesoscale updraft in the stratiform region (cf. Figure 1c). The updraft-rearward particles initiated between 5.75 and 7.25 km altitude are transported directly rearward and approach the melting level near the rear end of the stratiform region between x = 185 and 215 km. These trajectories have some spread that may be due to their different initial altitudes. Higher altitude particles have longer interaction with the mesoscale updraft, thus prolonging the time that they remain lofted while being transported rearward. These results are consistent with previous studies that found that ice crystals falling into the stratiform region originated from immediately rearward of the convective updraft at mid-to upper levels (Biggerstaff & Houze, 1991, 1993; Braun & Houze, 1994; Yang & Houze, 1995).

Figure 5. Shading represents the trajectory density (number of times a particle passes over each grid box non-consecutively) for the 2D simulations (a) SIM2D-CL, (b) SIM2D-DRZ, (c) SIM2D-DRZRA, and (d) SIM2D-RA. Black contours are line-average reflectivity factor at 10, 20, 30 and 50 dBZ, and gray contours are line-average vertical velocity at 0.5, 1, 2, and 4 m s−1. The colored lines are representative trajectories in a given simulation, and the black dot at the end of the trajectory is the particle's final position above the melting level. The gray square outlines the crystal initialization domain.

For the SIM2D-CL particles initiated within the updraft (green and yellow trajectories in Figure 4b), there are two distinct trajectory regimes. Particles initiated between 7.25 and 7.75 km altitude are transported to the forward anvil where they sublimate to extinction, and particles initiated between 5.75 and 7.25 km are lofted upward and then transported rearward, where they reach the melting level in the trailing stratiform region. This second regime of particles reaches the melting level at the same location as many updraft-rearward particles, but these particles do not have similar trajectories; rather, the updraft-centered particles follow a more straight-line, diagonal path to the stratiform region, suggesting that they have larger fall speeds earlier in their trajectories or that they encounter downdrafts or weaker updrafts in the transition zone. The maximum in trajectory density near z = 6 km and x = 240 km marks the overlap of updraft-rearward and updraft-centered trajectories, as shown by the two representative trajectories in Figure 5a. Finally, the SIM2D-CL particles initiated forward of the updraft (orange and red trajectories in Figure 4b) sublimate to extinction either immediately or after following an elliptical path to the forward anvil.

The trajectories in SIM2D-DRZ (Figure 4d) are nearly similar to those in SIM2D-CL with a few exceptions. One difference is that the updraft-rearward particles (dark blue trajectories) are now too heavy to be lifted by the mesoscale updraft into the anvils; instead, all of the particles reach the melting level in the stratiform region, slightly forward of the particles in SIM2D-CL. The relatively larger fall speeds in SIM2D-DRZ cause some of the particles initiated at lower altitudes to fall immediately into subsaturated air and sublimate to extinction just above the melting level at x = 225 km. These sublimated particles yield a gap in the trajectories between the updraft-rearward particles (dark blue trajectories) and the updraft-centered particles (yellow trajectories). A second gap in trajectories reaching the melting level appears at the transition zone (x = 240–260 km), though this gap results from a smaller group of quickly falling particles in the leading line due to their accelerated growth aloft. An example trajectory from this group of particles is shown in Figure 5b.

The trajectories of particles from SIM2D-DRZRA (Figure 4f) and SIM2D-RA (Figure 4h) have little dependence on the initial particle location. The majority of the particles approach the melting level toward the rear end of the convective reflectivity tower, with the SIM2D-DRZRA particles traveling farther rearward than the heavier SIM2D-RA particles. This finding is consistent with past studies asserting the importance of particle size sorting in generating the reflectivity structure of the squall line. In both SIM2D-DRZRA and SIM2D-RA, the updraft-forward particles (red trajectories in Figures 4f and 4h) initiated in subsaturated air are lofted up to 2 km by the convective updraft and travel farther than the updraft-centered and updraft-rearward particles, suggesting that their mass was reduced by partial sublimation. The overall result is a downward and rearward-extending maximum in trajectory density for both simulations (Figures 5c and 5d). The reduced trajectory density in the transition zone is consistent with the reduced snow and graupel mixing ratios in this region (cf. Figures 2f and 2g). Notably, some SIM2D-RA trajectories ending in the transition zone sublimate to extinction after prolonged periods in ice-subsaturated air (Figure 4h). The observed locations of widespread particle sublimation (underneath the leading anvil, above the stratiform melting level, in the transition zone) in all of the simulations are consistent with the locations identified in Braun and Houze (1994) and A. A. Jensen et al. (2018).

Representative Trajectories

The Lagrangian framework of the ICTG model allows us to learn how the particle properties evolve in response to their varied trajectories through different regions of the storm with different growth conditions. The representative trajectories shown in Figure 5 were selected to best represent common trajectory patterns that result in notable differences in the particle evolution and final particle properties. The five representative particles are named according to their simulation and initial size: P2D-CL1, P2D-CL2, P2D-DRZ, P2D-DRZRA, P2D-RA. In Figure 6, we examine time series of the environmental conditions and microphysical processes that give rise to each representative particle's properties. The goal is to understand the role of various microphysical processes in producing the modeled evolution of the representative particles.

Figure 6. Time series of particle (a,c,e,g,i,j) equivalent diameter, mass, accumulated rime mass, effective density, aspect ratio, and height for the representative particles shown in Figure 5 (b,d,f,h) ambient temperature, ice supersaturation, cloud water mixing ratio, and vertical velocity. The black dots at the end of each line is the value at the crystal's final position.

The P2D-CL1 and P2D-CL2 trajectories exemplify the two types of paths that the smaller particles (d ≤ 0.1 mm) follow toward the stratiform melting level. Particle P2D-CL1 is lofted by the updraft during the first 0.5 hr (Figures 6h and 6j). During this time, it is situated in strongly ice-supersaturated air (s_i up to 0.18, Figure 6d) and abundant cloud water (Figure 6f) in the convective updraft, but it only experiences minimal riming due to its low collection efficiency (Figure 6e). Around 0.25 hr, the particle's density begins to stabilize around 760 kg m⁻³ as it travels away from the abundant liquid cloud water (Figure 6f) in the convective zone. From 0.5 to 1.5 hr, the particle passes through the transition zone, where the ice-supersaturation remains near zero (Figure 6d), and the particle's diameter and mass (Figures 6a and 6c) remain relatively small (<0.25 mm and <0.02 mg, respectively). Just before reaching the stratiform melting level at 1.5 hr, the majority of the particle sublimates (due to ice-subsaturation), and its diameter shrinks by 0.16 mm.

In contrast, the updraft-rearward Particle P2D-CL2 has a different evolution. This particle does not experience an initial increase in size because it is initiated in nearly unsaturated air and low cloud water content (thus limited riming). This particle experiences moderate supersaturations (s_i < 0.08) throughout its entire trajectory, yet becomes 0.85-mm larger and 0.35-mg heavier than Particle P2D-CL1. As shown in Figure 6b, p2D-CL2 remains mostly in the planar crystal temperature range (−8 to −22°C), and this is consistent with the particle's oblate aspect ratio evolution (Figure 6i). At 1.25 hr, this particle's density decreases to 675 kg m⁻³ once the particle reaches the mesoscale updraft. Because the particle is unrimed, this density reduction is associated with dendritic growth (i.e., branching). Branching is further confirmed by the increasing ice supersaturation and steady air temperatures near −16°C, which falls in the dendritic growth temperature range (−12 to −18°C). On the other hand, Particle P2D-CL1 did not experience significant dendritic growth because it quickly passed through the appropriate temperature range (near −15°C) which had low supersaturations. Although the current iteration of the ICTG model does not include aggregation, the diversity of particles being transported into the stratiform region supports the likelihood of aggregating particles through mechanical interlocking of ice crystal branches. This finding implicitly supports past work associating the enhanced vapor deposition in the mesoscale updraft with a large prevalence of aggregates composed primarily of dendrites (Braun & Houze, 1994; Leary & Houze, 1979; Willis & Heymsfield, 1989). While aggregation is beyond the scope of the present study, future trajectory modeling studies will examine the effects of aggregation.

Because SIM2D-CL and SIM2D-DRZ have similar trajectory patterns, the next representative particle, P2D-DRZ, was chosen to be distinct from the previous two representative trajectories; thus, it does not follow the maximum trajectory density. Particle P2D-DRZ is initialized toward the front edge of the convective updraft and experiences an initial decrease in density to 375 kg m⁻³ during the first 0.1 hr, owing to light riming in the updraft. From 0.1 to 0.5 hr, ice-supersaturated conditions lead to light vapor deposition. From 0.5 to 3 hr, the particle's fluctuating density evolution suggests there are competing effects from riming and vapor deposition. From 2 to 3 hr, the particle actually passes through the convective updraft a second time, as seen by the trajectory loop (Figure 5b), vertical velocity (Figure 6h), and particle height (Figure 6j). Between 1.5 and 2.5 hr, riming dominates the relatively small growth in the particle's mass. After 2.5 hr, the particle continues to rime, but its growth is dominated by vapor deposition (Figure 6d), as indicated by the increasing particle density and the final particle mass (∼0.125 mg) being much larger than the accumulated rime mass (∼0.015 mg). This accelerated growth aloft results in the particle becoming too heavy for the updraft to support at 2.8 hr, resulting in a fast, 0.5-hr descent to the melting level at the rear end of the convective line.

Particle P2D-DRZRA and P2D-RA are both initiated at the same location as Particle P2D-CL1, yet they have considerably different trajectories due to their larger size. Instead of being lofted by the updraft, the particles immediately fall through the updraft, resulting in their shorter trajectory times (0.55 hr for P2D-DRZRA, 0.25 hr for P2D-RA). Their larger size also enables significant accumulations of rime mass. Particle P2D-DRZRA has a larger reduction in density (by 200 kg m⁻³) than Particle P2D-RA because the smaller initial size of P2D-DRZRA makes it more susceptible to density changes from riming. During the first half of its trajectory, Particle P2D-DRZRA experiences ice-supersations up to 0.09 in the planar crystal temperature range and thus simultaneously grows by vapor deposition into an oblate shape (Figure 6i). During the second half of its trajectory, conditions are ice-subsaturated down to −0.1, and Particle P2D-DRZRA sublimates during its descent to the melting level to nearly its initial size and mass at 0.5 mm and 0.06 mg, respectively. Despite riming, the strongly ice-subsaturated conditions down to −0.15 result in Particle P2D-RA sublimating to a smaller diameter (0.95 mm) and mass (0.4 mg) than it initially had.

We also tested the sensitivity of the 2D trajectories to a decrease in the initial density to 500 kg m⁻³ for larger (≥0.5 mm) crystals (not shown), and found that the trajectories shifted rearward by 6–12 km and residence times increased by ∼7 min. Overall, the particle properties remained similar to those with initial densities of 917 kg m⁻³.

3D Ice Crystal Trajectory Simulations

In this section, we expand our trajectory analysis to use the full three-dimensional output from the WRF simulation. This approach will introduce along-line variability in the trajectories and provide insight on how the storm's along-line variability impacts the spatial distribution of ice crystals and their properties. Moreover, the 3D trajectory simulations will show that the ICTG model can be utilized on 3D output.

Representative Trajectories

Figure 7 shows the trajectory density from the 3D simulations and select representative trajectories with the same starting positions as those in Figure 5. The representative particles are named with the same convention as those in the 2D simulations. Overall, the trajectory density patterns are similar to those from the 2D trajectories which were computed using a line-average of the squall line's structure. One notable difference is that the trajectory densities in SIM3D-DRZRA and SIM3D-RA extend up to 12 km altitude in the convective region, as opposed to 8 km altitude in SIM2D-DRZRA and SIM2D-RA, owing to the stronger updrafts in the 3D simulations, which become weaker when averaged for the 2D representation of the flow. The representative trajectories in all four simulations are also similar to their 2D counterparts, with the exception of P3D-DRZ, which travels a maximum of 12 km in the along-line direction (Figure 8a). However, the majority of the particles in all 3D simulations are not advected an along-line distance greater than 5 km (Figure 8b). In fact, the majority of the larger particles (≥0.5 mm) are advected no greater than 0.5 km in the along-line direction. These results show that the 2D simulations using line-averaged thermodynamic and kinematic fields maintain the salient trajectory patterns of fully 3D particle transport. We now examine whether this averaging also maintains comparable particle growth signals.

Figure 7. As in Figure 5, but for the 3D simulations.

Figure 8. (a) Representative trajectories from Figure 7, but in three dimensions. Their horizontal projections are plotted onto the along-line/across-line plane, as well as the reflectivity at 0.5 km altitude. The gray rectangle outlines the across- and along-line limits of the ice crystal initialization domain. (b) Normalized frequency of the maximum along-line distance traveled by ice crystals during their trajectories. Normalized frequency is defined as count of particles with a given maximum along-line distance divided by the maximum count for that simulation.

Figure 9 shows the time series of the representative particle properties from the 3D simulations. Particles P3D-CL1 and P3D-CL2 both show similar properties as the corresponding particles in the 2D simulations: P3D-CL1 ends up 0.9-mm smaller and 0.6-mg lighter than P3D-CL2 upon reaching the melting level. Although the first 0.75 hr of P3D-CL1's trajectory are in ice-supersaturated conditions, the temperatures are also extremely low (<−25°C), which slows vapor deposition growth rates. Even after the particle descends beneath 8 km into higher subfreezing temperatures (t = 1 hr), particle growth remains inhibited, now by supersaturation values that are near zero. The particle ends up with an aspect ratio >1, indicating a columnar crystal. A columnar crystal occurs here because Particle P3D-CL1 is immediately lofted ∼1 km higher than Particle P2D-CL1 by the stronger convective updraft in the 3D simulation. These higher altitudes are in the columnar temperature regime (<−22°C). Afterward, Particle P3D-CL1 enters the transition zone where potential growth into other habits is inhibited by the fluctuating levels above and below ice saturation.

Figure 9. As in Figure 6, but for the 3D simulations.

In contrast, Particle P3D-CL2 grows to be larger and heavier than P3D-CL1. P3D-CL2 grows into a planar crystal with an extremely low aspect ratio (down to 0.05), and it begins branching once it reaches the stratiform mesoscale updraft at t = 1.25 hr, as shown by the corresponding diameter increase, density decrease, and temperature of −16°C (within the dendritic growth temperature range). The particle continues to quickly grow thereafter as it encounters persistent ice-supersaturated conditions provided by the mesoscale updraft, reaching the melting level with a final diameter of 1.15 mm and a final mass of 0.6 mg. This result agrees well with past modeling and observational studies remarking on the faster deposition growth of ice particles within the mesoscale updraft of the trailing stratiform region (Braun & Houze, 1994; Rutledge & Houze, 1987).

Particle P3D-DRZ grows to be the third largest 3D representative particle, with a final diameter of 1 mm upon reaching the melting level. It becomes 0.13-mg heavier than its 2D counterpart, with most of its mass acquired by a combination of riming and vapor deposition when the particle is in warmer temperatures below 7 km (at t = 0.5 and 1.2 hr). Particle P3D-DRZ ends up oblate (in contrast to the slightly prolate P2D-DRZ), owing to reduced growth of P3D-DRZ while it is at higher altitudes in the columnar temperature regime (<−22°C).

Particles P3D-DRZRA and P3D-RA are both initially lofted ∼0.5 km by the stronger convective updrafts in the 3D simulations instead of immediately falling downward as their 2D counterparts did. They quickly rime while in the convective updrafts and then approach the melting level in the leading line instead of being transported to the stratiform region. Similar to their 2D counterparts, both particles have trajectory times between 0.25 and 0.5 hr and partially sublimate upon quickly falling into ice-subsaturated air above the melting level. However, the stronger subsaturations in the 3D simulation causes more of Particle P3D-RA's mass to sublimate than that of Particle P2D-RA, where P3D-RA has a 0.35-mg mass reduction at the end of its trajectory, in contrast to the 0.13-mg reduction for P2D-RA.

Thus far, the ICTG model has demonstrated the ability to realistically track the locations and changes of evolving ice crystals in the simulated storm. Differences in the particle properties were shown to depend on their initial size and location, which is consistent with past observational and modeling studies. The 2D and 3D simulations produce similar bulk trajectory patterns, but the properties of the cloud- and drizzle-sized particles in the 3D simulations are generally more extreme in magnitude than those from the 2D simulations. At appropriate temperatures, the more extreme supersaturations in the 3D simulations help generate considerably larger particles than in the 2D simulations, but there are also comparable regions of extreme subsaturation in the 3D simulations that can produce smaller particles through greater sublimation. The similarity between the 2D and 3D trajectory patterns likely arose due to a balanced interaction between the particle's fall speed and the simulated vertical and across-line wind speeds; the faster winds and more extreme particle growth in the 3D simulations balanced each other in the same way as the slower winds (due to line averaging) and less extreme particle growth in the 2D simulations. The ICTG model has thus demonstrated its versatility across different source model configurations.

Particle Recycling in the Convective Updraft

A small group of particles with a distinct trajectory type in both 2D and 3D simulations grew to significant sizes, some even larger than the selected representative particles. The updraft-forward upper-level particles from Figures 4a–4d, as well as Particles P2D-DRZ (Figure 5b) and P3D-DRZ (Figure 7b), all had elliptical trajectories with repeated passes through the convective updraft. We now explore the impact of these repeated passes (that we call particle recycling) on a particle's final properties upon reaching the melting level.

Figure 10 shows select particle trajectories from SIM3D-CL and SIM3D-DRZRA that have particle recycling patterns. Both particles make multiple passes through an updraft-downdraft couplet between 8 and 10 km altitude before reaching the melting level within the convective line. During these repeated passes, the particles grow by both vapor deposition and riming, as supported by the multiple peaks in ice supersaturation and cloud liquid water (Figures 10d and 10f) occurring at the same time as small increases in the particle diameter, mass, and accumulated rime mass (Figures 10c–10g). In these examples, the collected rime mass is an order of magnitude greater than that of the representative particles in Figure 9. Due to its larger size, the particle from SIM3D-DRZRA has a larger collision efficiency and accretes more cloud droplets than the smaller particle from SIM3D-CL. The fastest increase in both particles' mass occurs during their final descent toward the melting level, when the particles quickly collect rime (at t = 2.4 hr for SIM3D-CL, t = 1.2 hr for SIM3D-DRZRA), yielding a final mass of 2.8 and 1.5 mg, respectively. With repeated passes through the vertical velocity couplet, the particles have more time to grow (via vapor deposition) large enough to considerably collect cloud droplets once they fall into the higher temperatures below. It is worth noting, however, that riming only contributes to ∼20% of their total mass prior to partially sublimating at the end of their trajectory; the other ∼80% is from vapor deposition.

Figure 10. Representative trajectories showing particle recycling in the (a) SIM3D-CL and (b) SIM3D-DRZRA simulations. Cross sections of vertical velocity (shading) and reflectivity (black contours at 30 and 40 dBZ) are shown at y = 15 km in (a) and y = 19 km in (b) to match the initial along-line position of each particle and show the kinematic conditions contributing to their recycling trajectories. The gray box outlines the crystal initialization domain. Panels (c)–(h) show time series of the ice crystal properties and environmental conditions over the particle's trajectory, with the line colors corresponding to the trajectories in (a) and (b).

Bulk Particle Properties

We now examine the final properties of all particles from the SIM3D-CL and SIM3D-DRZRA simulations. These two simulations were chosen because they yield distinctly different final particle properties and have initial crystal sizes that are represented in the cloud and raindrop size distributions (cf. Figure 3). Figures 11 and 12 show the line-average final particle properties based on their initial locations in the crystal initialization domain, excluding particles that sublimate to extinction. This examination provides a more holistic analysis of the bulk particle properties than is possible by examining only select representative particle trajectories.

Figure 11. (a) The average final equivalent diameter of the crystals that did not completely sublimate, plotted on the crystals' initial position, using data from SIM3D-CL. Black contours are line-average supersaturation with respect to ice, and white contours are line-average vertical motion at the 0.5, 1, 2, and 4 m s−1 levels (b)–(f) As in (a), but for final mass, accumulated rime mass, final particle density, final aspect ratio, and total trajectory time. Total trajectory time is the time it takes the crystal to reach the melting level or 6 hr, whichever comes first.

Figure 12. As in Figure 11, but for SIM3D-DRZRA.

In SIM3D-CL (Figure 11), there are two initialization regions with particles that grow considerably larger and heavier than the rest: ahead of the convective updraft near 7.5 km altitude and rearward of the updraft near 6.75 km. Nearly all of the updraft-rearward particles (x < 277 km) grow exclusively by vapor deposition and have a larger density (700 kg m⁻³, Figure 11d) and lower aspect ratio (<0.6, Figure 11e), likely related to their lack of riming (Figure 11c). The largest (≥0.9 mm) and heaviest (≥0.6 mg) updraft-rearward particles are most similar to Particle P3D-CL2 (cf. Figure 7a), which travels directly rearward to the stratiform region and undergoes dendritic growth. In contrast, the relatively smaller updraft-rearward particles initialized near 7.5 km altitude are transported to higher altitudes, into the leading and trailing anvils, where the lower temperatures and higher characteristic supersaturations produce slower vapor deposition growth rates. This group of particles also has the longest trajectory times (>2 hr, Figure 11f), in agreement with the 2D simulations, which show that the particles entering the anvil regions do not complete their trajectories after 6 hr. The updraft-rearward particles initialized near 6 km do not grow as large as those initiated near 6.75 km due to prolonged exposure to ice-subsaturated air.

The updraft-centered particles (277 < x < 290 km) are the smallest (∼0.4 mm) and lightest (<0.1 mg) of all of the particles. The updraft-centered particles initiated above 7 km mostly travel into the leading anvil region and experience slowed growth. The particles initiated below 7 km have trajectories similar to Particle P3D-CL1 (cf. Figure 7a), whose growth is limited by low supersaturation and sub-saturation in the transition zone.

Most of the updraft-forward particles (x > 290 km) sublimate to extinction (not shown), but of those that remain, they grow relatively large. Their relatively large size (≥0.9 mm), mass (≥0.4 mg), and accumulated rime mass (≥0.1 mg) suggest that these particles undergo recycling, consistent with the initial locations of the individual particles in Figure 10 and Particle P3D-DRZ. The aspect ratio of these particles is near 0.8, and their density is between 400 and 500 kg m⁻³, suggesting that these particles have become small, graupel-like particles that reach the melting level in the convective line, in agreement with the spatial distribution of the graupel mixing ratio (cf. Figure 2g).

For the bulk particle characteristics of SIM3D-DRZRA (Figure 12), a gradient in the final particle properties occurs across the convective updraft. On average, the updraft-rearward particles (x < 277 km) generally end up smaller (<0.6 mm), lighter (<0.1 mg), and denser (>700 kg m⁻³) than the remaining particles due to growth at low supersaturation and limited cloud water. These particles are initially too heavy for the convective updrafts to considerably loft them into regions favorable for prolonged deposition or riming growth and thus do not grow as large or heavy as the initially smaller, updraft-rearward particles from SIM3D-CL. This finding is consistent with (Braun & Houze, 1994), who found greater growth of smaller ice particles (relative to larger particles) initialized at upper levels rearward of the convective updraft in their 2D model. As a result, these larger updraft-rearward particles follow quasi-straight paths down to the melting level, similar to Particle P3D-DRZRA. The particles initiated below 6.5 km generally have aspect ratios >1 as a result of them falling to lower altitudes and into a temperature range favorable for columnar growth (−8 to −3°C). This finding also agrees with (Braun & Houze, 1994), who found that needles were produced near 5.5 km altitude in their squall line simulation.

The updraft-forward particles (x > 290 km) grow relatively larger. Given that these particles are initiated in subsaturated air, they partially sublimate before reaching the convective updrafts. This early sublimation allows the updrafts to increase the particles' residence time in ice-supersaturated air and the probability of them coming into contact with supercooled liquid cloud droplets. The majority of the updraft-forward particles have a reduced density between 450 and 550 kg m⁻³, suggesting that they have some degree of riming or branching. Similar to SIM3D-CL, the particles initiated above 6.75 km have an increased size (≥0.9 mm) and mass (≥0.5 mg), though the particles from SIM3D-DRZRA have larger masses and collect more rime than the corresponding particles in SIM3D-CL. These particles are also likely graupel-like particles that resulted from recycling multiple times through the convective updraft. Their initial location, relatively large accumulated rime mass, and longer trajectory time suggest that they are similar to the recycling trajectory from SIM3D-DRZRA in Figure 10b.

The distinct differences in final particle properties between SIM3D-CL and SIM3D-DRZRA collectively show that capturing size-dependent evolution of crystal properties is important for simulating the microphysics in squall lines. Smaller particles ejected from the convective line have longer trajectories and are more likely to branch, have reduced densities, and end up in the upper parts of the cloud. In contrast, initially larger particles are more likely to heavily rime, have less spread in their trajectories, and reach the melting level sooner and closer to the convective line. Particle densities and fall speeds can have impacts on the kinematics through precipitation loading, and differences in the spread of trajectories hold implications for the latent heating structure of the storm. These notable differences in particle properties and transport highlight the importance of capturing accurate ice initialization in simulated squall lines.

Conclusions

This study introduces a novel ICTG model and evaluates its performance on a quasi-idealized WRF simulation of a leading-convective, trailing-stratiform squall line. The ICTG model includes recent laboratory-based parameterizations of ice crystal growth via vapor deposition and riming, as well as ice crystal decay via sublimation. The ICTG model tracks the trajectories of individual ice particles in space and allows the properties of those particles to evolve naturally based on their environment. By following the specific growth pathways taken by the ice crystals, the ICTG model simulates ice crystal characteristics within a squall line in a way not possible with Eulerian bulk microphysical models. The Eulerian modeling framework effectively averages the particle properties of a given ice crystal population, while the Lagrangian modeling framework maintains and tracks the evolving characteristics of individual particles throughout each particle's lifetime. This Lagrangian approach is important because ice crystals can develop a large diversity in crystal forms that affect their subsequent growth, fall speeds, and collection rates. Indeed, Nelson (2008) has estimated that the heterogeneity in ice crystal trajectories accounts for much of the estimated diversity in snowfall. Our modeling framework is designed to explore this diversity in ice crystal evolution and examine how such variety determines the microphysical properties of cold cloud systems.

Our analysis focused on 2D and 3D simulations of ice crystals initialized in the leading convective updraft of a squall line with diameters of 0.04, 0.1, 0.5, and 1 mm. In both the 2D and 3D simulations, the ICTG model was capable of producing an ice crystal trajectory pattern with particle properties that were consistent with past observational studies and the overall reflectivity structure of the simulated squall line. The particles with smaller initial diameters (d ≤ 0.1 mm) had trajectories that were generally organized by their initial across-line position relative to the convective updraft. Those smaller particles initiated immediately rearward of the updraft were transported to the stratiform region where they grew rapidly via vapor deposition in the supersaturated mesoscale updraft. These particles became the largest in diameter (≥0.9 mm) and mass (≥0.6 mg), while simultaneously evolving into secondary habits, consistent with past studies identifying increased dendritic growth in the stratiform region. The smaller particles initiated within the convective updraft ended up in the leading anvil or in the trailing stratiform region; these particles did not grow as large as the updraft-rearward particles, owing to their limited interaction with the mesoscale updraft. In contrast, the initially larger particles (d ≥ 0.5 mm) fell through the convective updraft and followed quasi-straight paths down toward the melting level, regardless of their initial position relative to the convective updraft. Their shorter trajectories in mostly subsaturated air resulted in these particles having a final mass 0.4-mg less than that of the initially smaller particles. In all of the simulations, the particles initiated along the front edge of the convective updraft were most likely to rime and reach the melting level closer to the convective line. Of these rimed particles, those initiated above 7.5 km became large (d ≥ 0.9 mm) graupel-like particles, owing to repeated passes through the convective updraft.

These results show that initially small particles (d ≤ 0.1 mm) can grow large either by vapor deposition when transported to the stratiform region or via recycling in the convective updraft, whereas initially larger particles (d ≥ 0.5 mm) with faster fall speeds have limited time and growth potential and thus remain primarily in the convective region. The ICTG model's production of varied particle trajectories based on size is consistent with past modeling and observational studies documenting size sorting of ice particles ejected from the leading convective line. One potential result of this size sorting is the commonly observed minimum reflectivity zone between the leading line and stratiform precipitation regions.

The above findings collectively reveal the importance of riming and vapor deposition in organizing ice particles in leading-line, trailing-stratiform squall lines. The ICTG model neglected aggregation and secondary ice production, yet produced a spatial distribution of ice particles in agreement with the simulated precipitation pattern. The ICTG model's production of realistic particle properties suggests that the model can be used to infer particle growth mechanisms in squall lines and other storm systems. This is a much-needed analysis because understanding the ice transport and growth in these systems at a basic level can improve forecasts of these storms as well as quantitative precipitation estimation.

Although the ICTG model is still under development, future use of the model can guide bulk microphysics scheme development. Differences between particle properties produced by the ICTG model versus Eulerian numerical models can reveal which assumptions and microphysical processes may be contributing to inaccuracies in precipitation features produced by bulk microphysics schemes. The Lagrangian nature of the ICTG model can provide further information on which regions of a storm bulk microphysics schemes struggle to reproduce and can identify the microphysical processes most responsible for potential inaccuracies. These comparisons can improve the microphysics of operational forecast models, leading to better forecasts of storm evolution as well as quantitative precipitation estimation. Future work will involve the addition of aggregation to the model, accounting for variations in particle number concentrations, and testing the robustness of the results using squall line simulations with different microphysics schemes.

Acknowledgments

The authors thank Matthew Kumjian and the three anonymous reviewers for their insightful comments that greatly improved this manuscript. C.N. Laurencin was supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1255832. A.C. Didlake was supported by the National Science Foundation under grant AGS-1810869. J.Y. Harrington was supported by the National Science Foundation under Grant AGS-1824243. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Data Availability Statement

Model output data are available on the Data Commons repository at Penn State University: https://www.datacommons.psu.edu/commonswizard/MetadataDisplay.aspx?Dataset=6299.