Mixed-phase Direct Numerical Simulation: Ice Growth in Cloud-Top Generating Cells

.


Introduction
It is well-recognized that the model representation of natural ice is not accurate due to a poor understanding of the physics related to, and therefore the poor model representation of, the ice and mixed-phase microphysical processes such as the ice 25 nucleation mechanisms, representation of supercooled liquid water, and the interactions between the ice-liquid phase hydrometeors (Korolev and Milbrandt 2022;Korolev et al. 2017;Koop and Mahowald 2013;Ovchinnikov et al. 2014). This poor model representation is exacerbated by the insufficient measurements used for constraining the microphysical parameterization at sub-grid scales (Fridlind et al. 2012;Baumgardner et al. 2017;Heymsfield et al. 2017;Korolev et al. 2017;Field et al. 2017). This uncertainty in the growth of ice can further impact the forecast of the supercooled liquid water 30 amount in mixed-phase clouds (e.g., Morrison and Pinto 2006;Barrett et al. 2017) via the Weber-Bergeron-Findeisen (WBF) mechanism (i.e., ice grows at the expense of liquid evaporation as a result of the difference in saturation vapor pressures with respect to liquid and ice, e.g., Storelvmo and Tan 2015) and thus the subsequent precipitation process.
The interactions between the ice particles and the supercooled liquid water in a turbulent environment is a multi-scale issue, and detailed information, such as the Lagrangian history of the particle growth, the spatial distribution of the particles of 35 distinct phases at fine scales are key information to the decay and maintenance of mixed-phase clouds (Korolev and Milbrandt 2022). Theoretical, modeling, and observational studies have been conducted to examine the influence of various factors on mixed-phase clouds such as the entrainment-mixing, WBF process, large eddies, vertical oscillations, and ice number concentrations (Hoffmann 2020;Korolev 2007Korolev , 2008Korolev and Field 2008;Khain et al. 2022;Morrison et al. 2012; Lee et al. 2021). Nevertheless, the native scales at which interactions between individual ice particles, liquid droplets, 40 and turbulence were rarely explored due to the limitation of the spatial resolution of the models and in-situ measurements.
Lagrangian-particle-based direct numerical simulations (DNS) can fill this knowledge gap in fine-scale microphysics and dynamics by resolving the turbulence and cloud particle growth.
DNS has been applied in cloud microphysics modeling since the early 2000s (Vaillancourt et al. 2001(Vaillancourt et al. , 2002 to study the turbulent impact on particle growth and spectral broadening. The bulk of the research has centered on warm-phase clouds 45 (e.g., Ayala et al. 2008;Grabowski and Wang 2013;Li et al. 2020;Gotoh et al. 2016Gotoh et al. , 2021Saito et al. 2019;Chen et al. 2018aChen et al. , b, 2020Chen et al. , 2021 with a limited number of studies investigating ice-phase clouds (Vowinckel et al. 2019a, b). In contrast, mixed-phase processes have received very little attention to date, despite their importance and complexity.
In this study, a Lagrangian-particle-based mixed-phase microphysics was implemented in the DNS model to untangle the interactions between ice, droplets, and turbulence at true native scales relevant to microphysics and to improve the 50 fundamental understanding of ice growth in mixed-phase clouds. This paper aims to present a mixed-phase DNS model with the capability of simulating the microphysical processes key to the development, maintenance, and decay of mixed-phase environments and the evolution of the ice and droplet size distributions.
Specifically, this study investigates the influence of cloud top generating cells (GCs) on ice growth in mixed-phase clouds.
GCs are fine-scale structures commonly observed in different types of clouds such as the winter orographic clouds 55 (Tessendorf et al. 2019;Kumjian et al. 2014), midlatitude cyclones (Keeler et al. 2016;Rauber et al. 2015), and mixed-phase clouds in the Southern Ocean (Wang et al. 2020). The GCs typically contain more liquid and ice and stronger updrafts than the adjacent portion of the clouds. Therefore, it is a key region for the formation of ice and precipitation. To resolve processes in GCs, an accurate representation of microphysics and dynamics are needed. The processes involve interactions between ice particles and liquid particles at centimeter scales. The structure of the flow at this length-scales is difficult to 60 measure using in-situ data and is also unable to resolve in mesoscale models. Therefore, a mixed-phase DNS is beneficial to obtain a process-level understanding at such fine scales.
With the overarching research question of "how can we use this innovative tool to design experiments and improve ice parameterizations in models that do not resolve mixed-phase microphysics explicitly?", our model development efforts aim to: 65 1) Present a DNS model studying the mixed-phase cloud microphysics with resolved fine-scale turbulence and Lagrangian cloud particles (aerosols, droplets, and ice crystals) and demonstrate the capability of such a model.
2) Investigate the impact of turbulence and environmental conditions on the ice growth and cloud glaciation at the native scales of droplet-ice interactions. We explored the potential mechanism of efficient ice growth in GCs.
Controlled experiments are conducted by the mixed-phase DNS to quantitatively examine the impact of the GCs 70 environment on the ice growth. In-situ measurement data from the Seeded and Natural Orographic Wintertime clouds -the Idaho Experiment (SNOWIE) (Tessendorf et al. 2019) are used to set up the initial conditions for the microphysics and environmental conditions which resemble those inside and outside GCs. Finally, the atmospheric implications of the simulations will be discussed.

Mixed-phase DNS model 75
As a part of our continued efforts to develop a comprehensive cloud microphysics model (Chen et al. 2016(Chen et al. , 2018a(Chen et al. , 2018b(Chen et al. , 2020(Chen et al. , 2021, this study presents the first version of a mixed-phase DNS model. The aim of this model is to address the intricacies of the interactions between droplets, ice, dynamics, and thermodynamics at microscales. The modeling framework of the warm-phase DNS model developed by Chen et al. (2018b) was adopted for this purpose. In this modeling framework, droplets are modeled as Lagrangian particles and are treated as inertial particles, subject to both turbulent flow and gravity. 80 The simulation domain is defined by a triply-periodic boundary condition. The following two model improvements were made to consider the ice depositional growth in a turbulent-mixing environment: (1) added depositional growth of spherical ice particles, and (2) added the external forcing on the thermodynamic fluctuation fields (i.e., temperature, ′, and water vapor mixing ratio, # ′) to account for the effects of turbulent mixing cascading down from larger scales. 85

Depositional growth of ice particles
The purpose of this model is to simulate the initial ice growth ( < 100 ) in mixed-phase clouds to examines the impacts of microphysical/environmental conditions on the production of larger ice particles (i.e., ≥ 100 ) and to understand the WBF process, which refers to the rapid growth of ice at the expense of the evaporation of droplets when supercooled droplets and ice crystals co-exist, but also impacted by fine-scale turbulence in clouds. In this study, the deposition growth of newly 90 nucleated ice crystals was considered in the model as the only mechanism for ice growth. Other nucleation and growth mechanisms such as the freezing of supercooled droplets, secondary ice productions, ice aggregation, and riming were not considered in this study. The ice depositional growth rate in mass of each individual spherical ice crystal was calculated in the model (Rogers and Yau 1989, equation (9.4)): / is the mass of a spherical ice crystal with a radius of r, = is the capacitance for spherical ice crystal, / is the supersaturation over the ice surface, and 0 = 2.83 × 10 1 !, is the latent heat required to deposit per unit mass of water vapor to ice. # = 467 !, !, is individual gas constant for water vapor, ′ and ′ are, respectively, the thermal conductivity and water vapor diffusivity that include kinetic effects (see equation 11a-b in Grabowski et al. 2011 for the detailed expressions). Note that the ice crystals are assumed to be spherical in this study and therefore the kinetic effects of the 100 diffusivity and conductivity on ice crystals are treated in the same way as water. Zhang and Herrington (2014) demonstrated that when ice crystals take non-spherical shapes, the kinetic effects can be very different, and their shape parameter needs to be considered. 02&,/ is the saturation water vapor pressure. The assumption of a circular ice particle may not accurately represent the true geometry of observed ice particles whose true habit depends on temperature and humidity (Korolev and Isaac 2003). Model development is undergoing to consider ice habits in future simulations. 105

Thermodynamic forcings on fluctuation fields
To account for the effects of turbulent mixing between different parts of the cloud cascaded down from scales larger than DNS domain size (i.e., ≫ 1 ), a forcing scheme was implemented to modulate the fluctuations of the temperature field ( ′) and the water vapor mixing ratio field ( # ′). In this way, statistical stationarity of thermodynamic fluctuations imposed by largescale turbulence is achieved: 110 In these equations is the flow velocity, ′ is the vertical perturbation velocity, $ ′ and / ′ are the differential droplet condensational rate and differential ice depositional rate, respectively, between the grid cell and the entire domain. & and # are thermal diffusivity and water vapor diffusivity. 5 and 9 % are forcing terms. The two fluctuation fields were forced in a 115 low-wavenumber band in the Fourier space, similar to the way the velocity field is forced (Chen et al. 2016). The energy cascades down to smaller and smaller scales via the non-linear interactions between eddies is eventually dissipated by viscosity and thermal diffusivity at millimeter scales. The forcing intensity, which measures the strength or efficiency of the turbulent mixing of the scalar fields ( and # ), is determined by the prescribed standard deviation ( , # ) which can be obtained from flight measurements: 120 Here 5 ) and 9 % ) are constants determined by the domain size of the simulation.

Numerical experiments
To set up the model initial conditions, in-situ measurements during the intensive observation period (IOP) 22 of the Seeded 125 and Natural Orographic Wintertime clouds -the Idaho Experiment (SNOWIE, see Tessendorf et al. 2019) was used. IOP22 was observed from 2017-03-09 14:00-18:00 UTC. We utilized the data from flight segment "A" to extract the environmental conditions of the generating cells. In segment A (14:08-14:10 UTC, with flight distance of ~ 12 km given a mean aircraft speed of 98.9 m/s), 13 GCs in total were identified. Six of the GCs had at least 5 seconds of flight measurements, which translates to a horizontal length scale of roughly 500 m or longer. We selected the time series of GC that lasted the longest in the 130 measurement period (i.e., the largest generating cell of the five) as the statistics were believed to be more representative of GC conditions than the rest of the timeseries. Table 1 summarizes the environmental conditions measured in this segment. As demonstrated, the observed number concentrations of drizzle and large ice particles (d > 100 ) in the GC environment are much higher than that of non-GC environments.
Two sensors (Rosemount and DMT100) were available to measure the water vapor mixing ratio which has uncertainty in the mean values (1.66 !, from DMT100 and 2.17 !, from Rosemount), therefore, both values are tested to account for the measurement uncertainty. To study the WBF process, a diagnostic relation between the supersaturation and temperature 140 was calculated using a simple steady-state parcel model to identify the range of relative humidity (RH) favoring the WBF process (as shown in the triangular shaded area in Figure 1). To ensure that (1) the tested range covers both conditions that are favorable or unfavorable for the WBF process, and (2) the corresponding mixing ratio at this condition reflects the measured range and its uncertainty (Table 1), we chose a range in supersaturation between : = −20% ∼ 5%, corresponding to RH = 80∼105% (purple stripe in Figure 1) at the given temperature and pressure condition of the GCs was chosen as the initial 145 conditions in the DNS. The microphysical response to varying moisture levels was therefore investigated.

Figure 1.
Diagnostic relationship of the temperature (x-axis), supersaturation over the liquid water (y-axis), and supersaturation over ice (color contours). The shaded gray triangular area indicates the condition in favor of the WBF process, and the purple column marks the mean temperature of the generating cell measured in IOP22 (T~ -13.6 C). The 150 range of supersaturation (-20%-5%, i.e., RH=80-105%) in the y-axis will be examined in the current study.
Six groups of experiments, including the control group, were performed, which consist of 18 simulations in total (see summary of the experiments in Table 2). Within each group, only one condition is varied and the rest remains the same. The initial conditions for the controlled experiment are listed in Table 2. Droplets have a dry radius of 1 μm with a hygroscopicity 155 parameter = 0.3. The purpose of conducting these sensitivity tests is firstly to validate the model on whether it could produce reasonable results, and secondly, by tuning different controlling variables, to understand what conditions favor ice growth and explore why the generating cell environment tends to produce more ice. The giant dry radius of 1 μm is not a typical aerosol radius, however it is used in this experimental setup for the purpose of examining the upper bound impact of the aerosol's hygroscopic effects on the ice growth and the WBF process. 160 The GCs were reported to contain a higher liquid water content (LWC) than the surrounding cloud body. Therefore, two groups of experiments were conducted, which highlight a high LWC and a low LWC condition, respectively. The high LWC experiments were initialized with LWC = 0.55 !" ( :;& = 10 , $ = 100 !" ), a typical value observed in a GC during SNOWIE, and therefore are named as RH_GC. Likewise, the low LWC experiments were initiated with LWC = 0.002 !" ( :;& = 1.5 , $ = 100 !" ), a much lower value than inside the generating cell, and are named RH_noGC (see 165 Table 2 for setup detail). The primary objective of conducting experiments in low LWC environments is to demonstrate the crucial role of large drops in ice growth under the same RH conditions, as well as the susceptibility of ice growth to changes in RH at varying LWC environments. Another two groups of experiments were conducted with different initial ice number concentrations. The first group is named High_IN in which the small ice number concentrations vary from 1 to 100 !" with the same initial radius of 1 ( Table 2).
The second group is named Low_IN with / = 8, 80, 800 !, . This range is commonly observed in many natural mixed-phase clouds. 175 To examine how turbulent-mixing-induced supersaturation fluctuations affect the ice growth and the glaciation process, the last group, named Group Turb, compares simulated results with and without turbulence. The model configurations are the same as the control group except that the experiment without turbulence has no fluctuations in the velocity, temperature, and mixing ratio fields. In other words, the model is treated as a traditional cloud parcel with no external forcing term and turbulent transportation term in the temperature and mixing ratio equations (Equation (2-3)) to generate supersaturation fluctuations. 180 To accurately capture all turbulent scales (i.e., scales down to Kolmogorov length scales) in the computational domain, the spatial resolution was set to 0.15625 mm. The time step was set to be very small (Δ = 0.472 × 10 !' ) to ensure that the smallest eddies and the particle relaxation time scales could be resolved. The droplet relaxation time scale defines how fast a droplet's surface temperature adjusts to equilibrium (see also the discussion of relaxation timescales in Vaillancourt et al. 2001). The domain size was 20 cm x 20 cm x 20 cm, and with a droplet number concentration of 100 !" , 800,000 droplets 185 were simulated and tracked within the domain.

Macrophysical impact
We first investigate the impacts of initial macrophysical conditions (i.e., mean quantities) on the ice growth. We emphasize the relative humidity and liquid water content, as the two are directly related to the water availability and immediately affect 190 the ice growth. Figure 2 shows the response of the mean particle radius and mean supersaturation to various initial RH. In the RH_GC group, the mean supersaturation adjusts to quasi-stationarity within 20-30 seconds and remains relatively unchanging for the duration of the simulation. (Figure 2 (c-d)). All air parcels become nearly saturated with respect to water and supersaturated over ice ( / ≈ 10%) and remain statistically steady thereafter. In fact, this quasi-stationary saturation ratio is maintained by an 195 unstable microphysical condition under which the WBF process (the ice grows rapidly at the expense of droplet evaporation) is active (Figure 2 (a-b)). The loss of water vapor due to the ice growth is continuously replenished by the evaporation of droplets, leading to a quasi-stationary RH.
In weather models or LES, a time step of 1-20 seconds is typically used with the assumption that liquid water is in thermodynamic equilibrium, but not for the ice. This assumption is based on the approximation that the phase relaxation time 200 for droplets is faster than that of ice. This simplified assumption that water condensation/evaporation is faster than ice processes is based on the single-phase cloud (i.e., warm cloud or ice cloud) where droplets or ice respond to moisture fields independently of one another. The phase relaxation time ( ?@20; , the time it takes for the saturation ratio to relax exponentially to 100%), can be calculated as ?@20; = (4 D A d ) !, , where D A = 2.55 × 10 !B 8 / is the diffusivity of water vapor, and d is the mean droplet radius. For a liquid cloud with a number concentration of = 100 !" , and a mean radius of d = 10 , ?@20; = 205 3.12 ; for an ice cloud with a much lower number concentration ( = 1 !" and = 1 in our case), ?@20; = 3121 .
However, this assumption does not hold true for evolving microphysics (i.e., number concentration and mean radius changing over time) or in mixed-phase environments where water vapor exchange occurs between ice and water drops. In mixed-phase conditions, ice, droplets, and thermodynamic fields interact concurrently. Both particles react to the thermodynamic fields in the same phase-relaxation timescale that evolves with time in the system, resulting in a slower reaction timescale to reach 210 quasi-stationarity than the one defined by the droplet phase relaxation time scale.
The RH_GC experiments demonstrate that the saturation at its quasi-stationarity is insensitive to the initial RH when at high LWC conditions (Figure 2c-d). During the quasi-stationarity, all ice grows at a relatively constant rate and is less sensitive to the initial relative humidity compared to droplet growth (Figure 2a). All air parcels tend to be glaciated, and the rate of glaciation depends on the initial RH. As Figure 3 demonstrates, a low RH parcel glaciates faster than that of a high RH but 215 produces a lower ice water content (IWC) due to less available total water vapor. In the atmosphere, when patches of droplets or drizzle drops are falling through or transported to a sub-saturated environment from elsewhere through processes such as strong turbulent-mixing and fall streak of hydrometeor (drizzle or melted snow) from the GCs, they can quickly replenish the surrounding air and create favorable conditions for ice growth. In the meanwhile, the cloud region with high RH environment, such as in the GCs due to its more upward motion than the surrounding environment, can maintain mixed-phase longer, but 220 also generate a high IWC due to the high availability of water vapor and liquid water.
In RH_noGC environments, due to insufficient liquid water in the parcel, ice growth rate is highly dependent on RH (Figure 2e).
For RH <90%, the air parcel is too dry to maintain ice growth due to low availability in both liquid water and water vapor, reaching a negative quasi-stationary supersaturation (Figure 2g). For RH>=90%, the air parcel is still supersaturated with respect to ice ( Figure 2g) and thus positive ice growth is maintained. The ice grows at a much slower rate compared to its 225 counterpart in the RH_GC group. It is shown that the air parcel generally glaciates faster regardless of its lower IWC (Figure 3c  The impact of RH on the glaciation rate is also prominent. Compared to the low RH cases, the air parcel initialized with a high RH glaciates slower but produces ice water faster (Figure 3). This is because a high RH can delay and slow down the WBF process. It is shown that maintaining a high RH is important for efficient ice growth while maintaining mixed-phase for a 235 longer period. In RH_GC group, droplets with RH=105% maintain growth for the first few seconds ( Figure 2b); in RH_noGC, droplets maintain high growth rates in the first 30 seconds (Figure 2f). This difference in the two groups is caused by the difference in the initial droplet wet radius. RH_noGC was initiated with droplets with wet radius of 1.5 and RH_GC with 10 , while droplets in both groups have the same dry radius of 1 . Following the -Kohler theory (Petters and Kreidenweis 2007) , the effective saturation ratio that considers the solute and curvature effects at the droplet surface is 02& = 240 ) is the curvature term. LK=# ≈ 1.0 in our case, and 02& ≈ 0IJK . The droplet growth rate is proportional to − 02& , where S is the environmental saturation ratio. Therefore, either a large S or a small 02& can lead to an effective droplet growth.
In RH_noGC, the droplets have a large dry radius a thin wet layer to begin with ( 0IJK&; ≈ 0.89 << 1.0 with :;& = 1.5 245 and $=> = 1 ), the solute term has an immediate strong impact on the growth rate in cases with RH≥ 100%, i.e., − 02& >>0. One can observe the prominent droplet growth in Figure 2f at the beginning, which is absent in RH_GC case (Figure 2b). When droplets grow to a larger size, :;& , the growth rate slows down due to a higher effective saturation ratio and a lower environmental saturation ratio (i.e., a smaller − 02& ). For example, for a droplet grows to :;& = 5 , 02& ≈ 0IJK ≈ 0.998. With the environmental saturation ratio becoming subsaturated at a later time ( − 02& < 0), a negative 250 growth rate is expected. In RH_GC, the solute effect is not important ( 0IJK ≈ 1.0 :;& = 10 ) and therefore, the growth rate is mainly determined by the environmental saturation ratio.
In summary, the value of supercooled liquid imposes an opposite impact on the ice growth and glaciation rate: a higher LWC can generate a higher IWC but results in a slower glaciation rate. To maintain both effective ice growth and mixed-phase for a long time, the environment needs either a high RH or a source of liquid water. Generating cells with a higher LWC and RH 255 than the adjacent portion of the clouds provide ideal environments for efficient ice growth. In the meantime, solute effects from large aerosols can be important at the initial growth stage to maintain droplet growth and thus the increased LWC; this has an important implication for future cloud seeding operations to a maintain high LWC for a more active WBF process. RH_noGC. The ice fraction is defined as the ratio between ice water content and total water content.
We also examined the evolution of particle size distribution (PSD) and found that it is highly sensitive to RH and LWC available in the air parcel. In RH_noGC, the RH plays a critical role in determining the mean radius: the higher the RH, the larger the ice mean size (Figure 4(b)). At RH<95%, most droplets evaporate, and ice growth is inhibited. Given that all ice particles have the same initial size, the broadening of the spectral width is solely due to the turbulent 270 fluctuation of the supersaturation field. The width of the ice spectral peaked at RH = 90% where the supersaturation with respect to ice is close to 100% (Figure 2g). Ice particles can experience positive or negative supersaturation in their immediate vicinity as they are advected by the turbulent flow. This results in a highly variable Lagrangian growth history within the ice population and thus the widest spectra. For cases with RH < 90% where / ≪ 0 when reaching quasi-stationary, the growth of the entire ice population is inhibited; for cases with RH>90% where / ≫ 0 at quasi-stationarity, the effect of mean 275 supersaturation dominates the fluctuations, and ice growth is less diversified. This can also be seen in RH_GC where / ≫ 0 in all simulations, the impact of initial RH on the mean ice radius and spectral width is negligible, and the ice PSDs are very    (Figure 4(d)). The theoretical study from Korolev (2008) found that for a typical LWC and IWC, ice growth is significantly less sensitive to the vertical velocity in mixed-phase clouds than the droplet growth. This finding is consistent with our finding, as the vertical velocity determines the saturation ratio (i.e., RH). 280 The effective impacts of supersaturation fluctuations on spectral widening can be quantified by the relative dispersion, which is a measure of the relative width of the PSD and is defined as the standard deviation of the particle radius normalized by its mean radius. In RH_GC, the relative dispersion of the ice size spectra only differs initially during the RH adjusting to quasistationarity, and converge to the same width at the end (Figure 5a). A mildly negative relation between the initial RH and relative dispersion of droplet size spectra is found (Figure 5b), indicating that supersaturation fluctuations, which broadens the 285 size spectra, play a more important role in a low sub-saturated environment. In RH_noGC, the spectral broadening for both ice and droplets is more sensitive to the initial RH. Different from RH_GC, the relative dispersion of droplet size spectra shows a reverse trend and increases with RH ( Figure 5d). This reverse trend is because most of the droplets evaporate in the low RHs and therefore reach a nearly-zero relative dispersion. 290 Figure 5. The relative dispersion of (top row) the ice size distribution and (bottom row) the droplet size distribution in air parcels initiated with different RH. The left column shows results from RH_GC, the right column from RH_noGC.

Microphysical impact
Next, we investigate the microphysical impact, such as the initial particle number concentration and particle size, on the ice growth. In Group High_IN, the range of ice number concentrations tested is much higher than the typical number 295 concentrations of natural ice observed in mixed-phase clouds, which usually is much lower than 1 !" in the atmosphere.

Group RH_GC Group RH_noGC
However, it is well recognized that measuring the ice particles below 100 is challenging and bears a high uncertainty. Ice particles with size <<100 can be largely under detected. This means that 1 !, of ice particles > 100 could be a result of ice crystal growth at a higher concentration. Therefore, the purpose of using this higher range is (1) to examine the sensitivity of resulting large ice particle number concentration to the initial ice crystal number concentration and (2) to test the upper 300 bound effects of ice number concentration.
Changing ice concentration exerts a huge impact. At / = 100 !" , ice particles quickly consume the water vapor ( Figure   6 c, d) and droplets completely evaporate after 50 seconds. Therefore, the parcel is glaciated within one minute (Figure 6(e)).
Though both IWC and total water content (TWC) stay the highest of the three simulations, the mean ice radius could not grow beyond 10 , as too many ice crystals or ice nucleating particles "overseeds" the cloud parcel and subsequently inhibits the 305 formation of precipitating-size particles. In contrast, with / = 1 !" , the parcel maintains mixed-phase much longer than the rest of the cases and sustains ice growth through the WBF process. Therefore, the formation of precipitating ice particles and the long-term persistence of mixed-phase conditions requires a relatively low ice number concentration ( / ≤ 1 !" ) and a high LWC environment. As is already demonstrated in High_IN above that for / ≤ 1 !" the parcel favors effective ice growth in a high LWC (0.55 !" ). In Low_IN experiments, the simulations were initiated with a lower LWC environment by reducing its droplet number concentration by a factor of ten (LWC = 0.055 !" , $ = 10 !" , $ = 10 , see also Table 2) to examine 315 how ice number concentrations influence the ice growth in a low LWC. Unlike the High_IN, when / ≤ 1 !" the growth of ice and evaporation of droplets were insensitive to the ice number concentration (Figure 7a-b) compared to High_IN.
However, the higher number concentration can lead to a much higher glaciation rate (Figure 6e). The mean supersaturation only starts to diverge after droplets nearly evaporate around 40 sec (Figure 7c-d). And ice grows in general is slower than in High_IN shown in Figure 6a albeit that fewer ice particles are competing for water vapor. This implies even if the ice number 320 concentration is low (<< 1 !" ), a high LWC is necessary for rapid growth of big ice particles through diffusional growth.

Turbulent mixing impact
To further examine how turbulent-mixing-induced supersaturation fluctuations affect the ice growth and the glaciation process, 325 we compare the results of the controlled experiment with the case without the turbulent mixing effects on the thermodynamic fluctuations. By comparing the two idealized simulations we can isolate the impact of turbulent-mixing on the microphysics at local scales (≪ 1 ). except for a slightly smaller mean droplet radius in the thermodynamically-forced case. The mean radius is calculated based on the arithmetic average, it is found that the volume-weighted mean radius is not altered by the fluctuation (not shown), substantiated by the identical LWC in Figure 8  Albeit the same mean properties in both cases, the size distribution width broadening is significant. In our case, the particle distribution is monodispersed in the beginning. Therefore, the broadening is solely caused by the saturation fluctuations. We conducted another experiment without supersaturation fluctuations for comparison. It can be seen from the relative dispersion 340 that no broadening happens in the case without supersaturation fluctuations (Figure 9) because particles are exposed to a fluctuated supersaturation environment. This supersaturation-fluctuation-induced broadening can be further substantiated by looking at spatial distribution of the supersaturation. Figure 10 displays the horizontal cross-section of the supersaturation fields at T=41 second in the RH_GC for RH=80%. On average, this case manifests an active WBF process, with a mean negative supersaturation with respect to water (~-0.3%) and mean positive supersaturation with respect to ice (~15.0%) (Figure 2c-d). This leads to effective ice growth and droplet evaporation. However, Figure 10 shows turbulent mixing creates pockets of positive supersaturation with respect to 350 water (with maximum value of 2.3% at this time point). This locally high supersaturation keeps the droplets from evaporation and even favors the droplet growth as shown in the PSD evolution in Figure 4(c). To quantitatively evaluate the spatial correlation between the droplet, ice, and supersaturation, Figure 11 (a-b) shows the Pearson correlation between the particle radius and the supersaturation value at the particle's location (i.e., the Lagrangian supersaturation). For droplets, the correlation between the size and supersaturation field is always positive (with p-value < 0.001). Correspondingly, the supersaturation 355 fluctuations are strong in the first 10-20 seconds, leading to the strongest broadening. For ice particles, the correlation is positive before around 10 seconds (with p-value < 0.001) and becomes insignificant afterwards. Therefore, compared to ice particles, droplet spectral broadening is more sensitive to the supersaturation fluctuation fields. This can also be substantiated by the evolution of size distributions in Fig. 4 (c-d): the ice spectral broadening due to supersaturation fluctuations is mostly significant in the first 10 seconds while the droplet broadening is notable for the first 50 seconds. This effective broadening of 360 ice and droplet particles within such a small domain again demonstrates the importance of small-scale fluctuations that could not be resolved by the commonly used cloud models or even large eddy simulations.

Summary and discussion
In this study, we developed the first DNS model capable of simulating mixed-phase cloud microphysics at its native scales and presented the first results of such a model. Six groups of experiments with initial conditions resembling the typical 375 environments inside and outside of the cloud-top generating cells were conducted. By changing the environmental (macrophysical and turbulent) conditions and microphysical conditions, various potential factors affecting the ice growth are explored to unravel the favorable conditions for ice growth in GCs.
In summary, high LWC or high RH are critical to maintaining both effective ice growth and mixed-phase status. Generating cells with higher LWC and higher RH than the adjacent portion of the clouds provide ideal environments for efficient ice 380 growth. In the meantime, solute effects from large aerosols can be important at the initial growth stage to maintain droplet growth and thus the increased LWC, this has an important implication for future cloud seeding operations to maintain high LWC for a more active WBF process.
The results show that in high LWC conditions, the saturation field reaches quasi-stationary within 20 to 30 seconds. This quasistationarity is maintained through an unstable microphysical condition, in which the loss of water vapor from ice growth is 385 continuously replenished through the evaporation of droplets (i.e., the WBF process), leading to a quasi-stationary RH. The timescale for reaching quasi-stationary in this mixed-phase environment is much slower than in a pure liquid environment. As a result, in a mixed-phase condition, where ice, droplets, and thermodynamic fields interact concurrently, the phase-relaxation timescale is not a reliable measure of the thermodynamic response time.
Sensitivity studies on ice number concentrations show that the ice number concentration below 1 !" , a typical range in the 390 atmospheric clouds, a high LWC is needed for an efficient formation of big individual ice particles through diffusional growth.
In addition, the ice radius growth rate is less sensitive to the ice number concentrations. However, a high ice number concentration > 10 !" results in an early termination of the WBF process, fast glaciation of mixed-phase clouds, suppression of the ice growth, which does not favor the formation of hydrometeors. This finding has an applicational implication on the winter-cloud seeding. 395 By comparing the cases with and without turbulent-mixing, it is found that small-scale thermodynamic fluctuations substantially broaden the particle size spectra. It is found that droplet spectral broadening is in general more sensitive to turbulent fluctuations; ice spectral width remains very similar when it is supersaturated with respect to ice, and broadening of ice size distribution is mostly prominent when the environment is close to ice saturation. The effective broadening at a time scale as fast as a few seconds and at a length scale down to cm-scales demonstrates the importance of small-scale fluctuations 400 that could not be resolved by the commonly used cloud models or even large eddy simulations. However, the simulation also shows that the fluctuations have a negligible effect on the mean properties such as the mean radius and mean saturation ratio.
The current study represents the first attempt to use the DNS model to understand the mixed-phase cloud microphysics at local scales. This first result reveals the capability of such an ultra-fine-resolution model to address the outstanding problems in mixed-phase cloud research. It bears some limitations and requires future work to fill the gaps. Firstly, studies that consider 405 collisions between droplets and ice particles are needed to further investigate the turbulent effect on the evolution of particle size distribution through collisions. Secondly, the mixing of air parcels of different properties not only affects the fluctuations but also influences the mean properties. In addition to the turbulent mixing effects on thermodynamic fluctuations, impacts on the mean temperature and mixing ratio also need to be considered to reflect the realistic entrainment-mixing process.
Observations on the cooling rate, variation of mixing ratio at scales much larger than the DNS, and dynamical and microphysical information along the Lagrangian trajectory in a 3D high-resolution cloud model such as LES can help to better constrain the DNS for more realistic investigation on the mixed-phase processes at the microscale. Thirdly, the ice geometry information is important to determine its fall velocity and diffusional growth rate and needs to be considered to better represent the evolution of ice particle size distribution (Korolev et al. 1999;Fukuta 1969;Avramov and Harrington 2010;Chen and Lamb 1994;Jensen et al. 2017;Jensen and Harrington 2015;Morrison et al. 2020). Ongoing modification is conducted to 415 improve the model. With the current model capability and ongoing development, the future application of the DNS on ice and mixed-phase microphysics can also be promising in reproducing or comparing with the results in the laboratory experiments (Desai et al. 2019), which can provide process-level understanding of the experiments and help design new experiments and new laboratory facilities with the well-validated model (Shaw et al. 2020). 420

Acknowledgments
This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. This work is also supported by Idaho Power Company. We would like to thank Jeff French, Wojciech Grabowski, and Anders Jensen for their valuable discussion. LX thanks Paul Field for the inspirational discussion on this topic that led to the development of the mixed-phase 425 DNS capability. Computing resources were provided by the Climate Simulation Laboratory at NCAR's Computational and Information Systems Laboratory (CISL).

Code and data availability
The data of this study are available to download in Harvard Dataverse (Chen 2022). Figures were made with matplotlib on the Jupyter Notebook. Part of the software associated with this manuscript for the 430 visualization is also stored in Harvard Dataverse (Chen 2022).