Process-Tracking Scheme Based on Bulk Microphysics to Diagnose the Features of Snow Particles

We have developed a new method of diagnosing the character- istics of ice particles using a bulk microphysics model. Our model tracked the mass compositions of different classes of ice particles, using their microphysical process of origin, such as water vapor deposition and riming. The mass composition from depositional growth was further divided into six components by the temperature and humidity ranges corresponding to the typical growth habits of ice crystals. In test simulations, the new framework successfully revealed the influences of riming and depositional growths of ice particles within clouds and on surface snowfall. The new approach enables weather prediction models to provide much more information on the characteristics of ice particles regarding crystal habits and the extent of riming. ( Process-tracking scheme based on bulk microphys- ics to diagnose the features of snow particles. SOLA


Introduction
During the last decade, new microphysical models were developed to incorporate the various physical properties of solid hydrometeors. Hashino and Tripoli (2007, 2011a, 2011b developed a sophisticated spectral-bin microphysics model that represents variations in ice crystal habits, based on Chen's (1992) model. Brdar and Seifert (2018) developed a Monte-Carlo particle model to represent riming and aggregation of ice particles. Shima et al. (2019) extended their superdroplet method (Shima et al. 2009), which is a particle-based, probabilistic simulation model of cloud microphysics, to include mixed-phase microphysics.
Recent advances in computational power enable such models to work under realistic atmospheric conditions, although the computational cost for this remains unfeasibly high. Bulk microphysics models provide a cheaper alternative to spectral-bin and particle-based models. Over the last few decades, developments and improvements have been made to bulk microphysics models to better understand the mechanisms of a cloud and precipitation system, and thus improve numerical weather forecasting (Lin et al. 1983;Reisner et al. 1998;Stoelinga et al. 2003;Morrison and Milbrandt 2015).
The capability of a bulk microphysics model to represent the various properties of solid hydrometeors and mixed-phase microphysical processes improves when one increases the number of classes of solid hydrometeors and the number of prognostic variables implemented in each class. Typically, most modeling studies have adopted visually distinguishable particle types as classes and physically observable properties as variables. Two strategies have been implemented to produce more sophisticated bulk microphysics models: (i) the multiclass strategy adds a new hydrometeor class and (ii) the multivariable strategy adds one or more variables to a hydrometeor class.
Using the multiclass strategy, Lin et al. (1983) added a new class for snow to Orville and Kopp's (1977) model, which had four classes -cloud water, rain, cloud ice, and hail -and solved a time evolution equation for the mixing ratio for each class (namely a single-moment scheme). Their work established the basic system of hydrometeor classification that is still used in bulk microphysical modeling. Currently, many bulk microphysical models have five or six classes of hydrometeor.
Using the multivariable strategy, Reisner et al. (1998) solved time evolution for the number concentration in every class of solid hydrometeor in addition to the mixing ratio (namely a doublemoment scheme) to diagnose the volume-mean size of a solid hydro meteor, which is crucial to accurately calculating microphysical processes. Milbrandt and Yau (2005) proposed a triplemoment microphysics scheme that solved the time evolution of the sixth moment of the size distribution function, in addition to the mixing ratio and particle number concentration in each class, to determine a new size distribution parameter associated with the dispersion of particle size with the goal of obtaining more accurate solutions for the microphysical equations. Recalling that the number concentration and mixing ratio are essentially interpreted as the zeroth and third moments of the number size distribution function, respectively, developments using the multivariable strategy were limited in improving the diagnostic method for the size distribution parameters until the work of Milbrandt and Yau (2005). Morrison and Grabowski (2008, hereafter MG08) proposed a new concept for bulk microphysics modeling. They adopted a single-class framework of ice, instead of using the traditional classification system. The mixing ratio for ice was divided into that originating from vapor deposition and that from riming. They were tracked separately, which enables us to know quantitatively the contributions from particular microphysical processes to the predicted characteristics of ice particles. Further developments were made to their model (Morrison and Milbrandt 2015, hereafter MM15), creating new trends in the numerical simulation of a cloud and precipitation system (Morrison et al. 2015;Milbrandt and Morrison 2016).
In our new approach, we developed a bulk microphysics module in the Japan Meteorological Agency Non-Hydrostatic Model (JMA-NHM; Saito et al. 2006) to examine the characteristics of solid hydrometeors in much more detail. We adopted a multivariable strategy, by adding several new variables representing microphysical processes to hydrometeor classes, as was done in the MG08 and MM15 models, but the development enables us to diagnose ice crystal growth habits within the traditional classifications. Although, the currently diagnosed characteristics do not influence the original solution of microphysical equations (including gravity sedimentation and collision efficiency), the new approach potentially expands the applicability of the bulk microphysics model in research collaborations with microphysical observations. The results from our test simulations using the improved JMA-NHM demonstrate the effectiveness of our approach.

Numerical model
We applied the JMA-NHM to our numerical experiments, using the original configuration with the following exceptions.
Here, S is defined as (Q c + Q v )/Q vsi , where Q c , Q v , and Q vsi are the mixing ratios of cloud water, water vapor, and saturated water vapor with respect to ice, respectively. By using the expression P xdep in Ikawa and Saito (1991), the PRD is defined for the depositional growth as follows: where x = i, s, and g. For ice nucleation, PRD is defined by the same stratification as shown above. Figure 1 shows these subclasses with colors representing the class of cloud ice or snow. The effects of ice nucleation on the mass of solid hydrometeors are quite small compared with the effects of riming and deposition. Therefore, for the sake of simplicity, in this study, we will show only the results associated with riming and deposition. The classification of the addressed riming and depositional growth processes is summarized in Table 1.

Numerical simulations
Test simulations were conducted in a two-dimensional domain. The domain measured 600 km horizontally and 12.6 km vertically (Fig. 2a). The domain contained a 1500-m-high and 150-km-wide bell-shaped mountain centered at x = 400 km. Horizontal grid spacing was 1 km, and the vertical grid spacing was 300 m. We used 48 vertical layers in a terrain-following coordinate system. The timestep was 4 s. We performed 13 simulations. The initial times are listed in Table 2. The integration time was 60 h, which included a period of snowfall events observed in the Echigo mountains (37°03′19.9″N, 139°00′11.4″E) during January 1995 (Harimara and Nakai 1999). Initial conditions for the simulations, air temperature, and relative humidity profiles for the middle of Sea of Japan [40°N, 135°E] are taken from the JRA-55 reanalysis data (Kobayashi et al. 2015) for each timepoint listed in Table 2. A horizontal wind of 10 m s −1 toward the right (Fig. 2a) was applied across the domain at the beginning of every simulation. It was First, we adopted a double-moment scheme for all classes of solid hydrometeor, such as cloud ice, snow, and graupel, while the original configuration adopts a double-moment scheme only for cloud ice and a single-moment scheme for the other hydrometeors such as snow and graupel (Ikawa and Saito 1991). For liquid hydrometeors (such as cloud water and rain), we adopted a single-moment scheme to predict only the mixing ratio of particles. Second, we removed the original configuration's ice-saturation adjustment scheme (Tao et al. 1989) to avoid the unrealistic formation of ice clouds in the upper troposphere (Hashimoto et al. 2007). Third, we removed a cumulus parameterization. Finally, we implemented new formulations to track microphysical parameters representing features of ice particles, as described in the next section.

Formulation
When a hydrometeor class x that has mixing ratio Q x (kg kg −1 ) acquires water mass through a microphysical growth process p, we define Q p, x (kg kg −1 ) to represent the partial mass originating from the process p. In this study, we address three microphysical growth processes for cloud ice (x = i ), snow (x = s), and graupel (x = g) due to riming ( p = Rim), water vapor deposition ( p = Dep), and ice nucleation ( p = Inu) processes. The Q p, x changes in dynamical transport and microphysical conversions with water vapor or other hydrometeors. The prognostic equation of Q p, x is expressed as follows; where ADV and DIF represent the time tendencies of the advection and diffusion processes, respectively. For the sedimentation term SED, the partial mass Q p, x is assumed to move downward at the same decent velocity as the total mass Q x . EXC x, y represents the conversion rate of Q p, x between the categories x and y through a microphysical transition.
where ( ¶Q x / ¶t ) x to y represents the mass transition from categories x to y, which is determined in accordance with the microphysical processes implemented in our model as shown in Fig. B-11-1 of Ikawa and Saito (1991). Equation (2) means that the mass conversion from one category to another is assumed to be proportional to the mass fraction of each microphysical component. The PRD represents a production rate for the riming, deposition or ice nucleation processes.
where P i.iacw and P g.iacw are the increasing rates in the mixing ratio of cloud ice and graupel, respectively, due to accretion of cloud droplets by cloud ice; P s.sacw and P g.sacw refer to the increasing rates in the mixing ratio of snow and graupel, respectively, due to accretion of cloud droplets by snow; and P gacw refers to the riming growth rate of graupel through the accretion of cloud droplets by graupel itself. These specific forms are described in Ikawa and Saito (1991). Because ice crystal growth habit depends on air temperature and humidity during depositional growth (Nakaya 1954;Kobayashi 1961;Bailey and Hallet 2004, 2012Takahashi 2014), we stratify the production term of depositional growth with air temperature T (°C) and supersaturation ratio S (%) with respect to ice, following the Nakaya's T as diagram (Nakaya 1954). Fig. 1. Color scale representing the hydrometeor classes or the subclasses for different growth modes. Navy blue shows cloud water or rain, and red shows graupel. These are the original hydrometeor classes in the JMA-NHM. Orange shows the riming growth mode for cloud ice or snow. Other colors show the depositional growth modes of cloud ice or snow: blue, irregular needle; sky blue, needle/scroll/cup/column; light green, sector/ plate; green, dendritic; and yellow, column. Gray represents all the others including irregular types, which appear at low temperatures. applied at the upwind boundary throughout each simulation. The surface fluxes of heat and humidity were given according to surface states (sea, blue; land, green; and snow cover, white in Fig.  2a). Output data were available every 10 min from the surface and every 60 min from the atmosphere. Figure 2a shows an example of the results on accumulated precipitation amount from the forecast time (FT) of 36 to 60 h Table 1. Variables representing the partial mixing ratios for cloud ice, snow and graupel originating from the riming process or from the deposition process at different ranges of temperature T and supersaturation ratio with respect to ice S. In addition to the 21 variables listed in the table, 18 variables were implemented for representing the compositions of ice nucleation.   Table 2). Colors show the contributions to precipitation from the hydrometeor classes or the subclasses in cloud ice/ snow in accordance with Fig. 1. (b) Vertical profiles of the mass contributions of the hydrometeor classes or the subclasses in cloud ice/snow at five points in the domain: x = 180 km, 280 km, 370 km, 400 km, and 430 km, which are averaged during FT = 36−60 h in the same simulation as in (a). Solid and broken curves show the total (Qt ) and cloud water (Qc ) mixing ratios, respectively, in g kg −1 . in the simulation with an initial time of 15:00 JST on 17 January 1995. Graupel (red) shows a large contribution to the total amount of precipitation over the sea and near the coast. This is consistent with the statistical analysis of Mizuno (1992). Moving inland, the contribution of graupel decreases, and components of riming (orange) and deposition (blue, sky blue, light green, or green) growth of snow increase. This means that partially rimed snow becomes the major contributor to snowfall. On the lee side of the mountain, most of the snowfall originates from depositional growth. The depositional growth at temperatures from −20°C to −10°C (green or light green) shows the largest contribution over the top or lee side of the mountain. Figure 2b shows the averaged in-cloud characteristics for the same period as in Fig. 2a. Graupel (red) and the component of riming growth of cloud ice and snow (Q Rim, i + Q Rim, s , orange) show a large contribution among solid hydrometeors over the sea (x = 180 and 280 km) and the upwind side of the mountain (x = 370 km). The columnar deposition over both warm (Q Dep-4, i + Q Dep-4, s , light blue) and cold (Q Dep-20, i + Q Dep-20, s , yellow) temperature ranges has a part of the mass of cloud ice and snow over the upwind side of the mountain (x = 370 km). The dendritic growth mode (Q Dep-14, i + Q Dep-14, s , green) shows a large contribution around −15°C over the upslope (x = 370 km) where updraft and, consequently, ice supersaturation increases. Implementation of the formulations described in Section 2.2 successfully revealed the microphysical characteristics of surface snowfall and solid hydrometeors in snow clouds. In general, this kind of data has not been produced using numerical weather prediction models before. One of the practical benefits of this improvement is to reveal the temporal-spatial distribution of particular types of snow particles that increase the risk of snow avalanche (Hirashima et al. 2018).

Discussion
In this section, we confirm the validity of the microphysical features simulated by the new microphysical data. Harimaya and Nakai (1999) used data obtained during seven snowfall events in January 1995 (Fig. 3a) to show the relationship between snowfall intensity and riming proportion (the ratio of riming mass to total mass constituting a solid hydrometeor, Supplementary material S1). They stated that:

Riming proportion in surface snowfall
Although the riming proportion varies greatly when the snowfall intensity is weak, the lower limit of the riming proportion tends to increase, and the blank space without solid points tends to expand with an increase in snowfall intensity. Therefore, stronger snowfall intensity does not occur when there is a low riming proportion. In other words, weak riming suppresses large snowfall rates. Figure 3b shows scatter plots of the same parameters but obtained from simulations. The lower limit of the simulated riming proportion increases with increasing snowfall intensity, as in the observed result. However, the upward variability of the simulated riming proportion (Fig. 3b) is smaller than the observed result (Fig.  3a), especially when there is intense snowfall.
These results show that the model can reproduce the observed feature: weak riming suppresses large snowfall rates, but the riming rate is underpredicted in the model as snowfall intensifies. The modeled riming rate is quite sensitive to the rate of ice nucleation. A large ice nucleation rate tends to suppress the riming rate because the generated ice crystals consume cloud droplets through the Bergeron-Findeisen process. Currently, the ice nucleation process is not well understood. More study through a combination of field observation, laboratory experimentation, and modeling is required to understand its impact.

Depositional growth mode of ice particles
The microphysical characteristics of winter orographic snow clouds were observed by a hydrometeor videosonde (HYVIS, Orikasa and Murakami 1997) launched from the upwind side of the Echigo mountains in Japan for 9 years from 1994 to 2002. The shape, maximum dimension, and water content of hydrometeors were derived from detected images (Supplementary material S2). Figures 4a, 4b, and 4c show the observed relationships between temperature and the mass content of different crystal habits. A columnar crystal (Fig. 4a) shows double peaks (−5°C and −25°C) in the mass content profile. A dendritic crystal (Fig. 4b) shows a maximum from −15°C to −10°C, and a planar crystal (Fig.  4c) shows a relatively flat profile. Figures 4d, 4e, and 4f provides the corresponding simulation results at x = 370 km, which mimics the HYVIS launching site. The simulated mass content of the columnar growth mode (Fig. 4d) shows double peaks (−9°C and −23°C) as in the observation result. The peak at the colder temperature originates from crystal growth in the temperature range from −36°C to −20°C (yellow in Fig. 1). The peak at the warmer temperature originates from crystal growth in the temperature range from −10°C to −4°C (sky blue in Fig. 1). Profiles for the dendritic (Fig. 4e) and planar (Fig. 4f) growth modes show almost the same characteristics as in the observed profiles at temperatures warmer than −20°C. The peak around −15°C in Fig. 4e comes from the crystal growth under the condition of −17°C < T < −14°C and S > 7% (green in Fig. 1). In contrast, for temperatures below −20°C, these profiles show a sharp decrease with increasing altitude (Figs. 4e and 4f), while the observed profiles show a gentle decrease (Figs. 4b and 4c). The classification of partial mass originating from depositional growth in our model currently follows Nakaya's T as diagram (Nakaya 1954), where columnar crystals are predicted below −20°C. Bailey and Hallet (2004, 2012 reviewed crystal habits in the temperature range from −20°C to −70°C based on their laboratory experiment and in-situ observation. They gave examples showing that many plate-like crystals were observed below −20°C, in contrast to Nakaya's diagram. In light of this, the discrepancy between the observation and simulation may indicate the necessity of updating our classification method. To renew this model in the near future, we will consider newer work, such as Bailey and Hallet (2004, 2012 and Takahashi (2014Takahashi ( , 2018.

Conclusion
In this article, we proposed a new bulk microphysics framework to examine the microphysical properties of ice particles in the atmosphere. We demonstrated its effectiveness in an investigation of the microphysical features of a cloud and precipitation  Harimaya and Nakai (1999). The black curve represents the lower limit of the riming proportion. Numerals in the upper row show the ratios (%) of snow particles with a rime mass of over 50% to the total snow particles. (b) Same as (a) but for the simulation results at a point on the upwind side of the mountain (x = 370 km). Parameters from the 13 simulations are plotted every 10 min during 24 h from FT = 36 to 60 h. The broken curve in (b) is the same as the solid black curve in (a).
system. In the simulations, the new framework successfully revealed the influence of the riming and depositional growth processes not only on the mass of ice particles suspended in the atmosphere, but also on snowfall at the ground surface. The simulations reproduced the relationship between the riming proportion and snowfall: that a larger snowfall rate on the ground requires more active riming in a cloud. However, the simulations failed to reproduce a large variation of riming proportions. One of the explanations for this is the uncertainty regarding the parameterization of ice nucleation. The simulated mass content profiles for the different depositional growth modes of ice crystals essentially concurred with those obtained by balloon-borne measurements. For thoroughness of modeling, our model needs the implementation of microphysical feedback reflecting the diagnosed characteristics of hydrometeors into the original equations. The proposed approach provides a new perspective on the validation of a bulk microphysics model using microphysical observations.