Evaluation of Cloud Microphysical Schemes for Winter Snowfall Events in Hokkaido: A Case Study of Snowfall by Winter Monsoon

The performances of bulk microphysical schemes were evaluated through comparisons with observations, including a new volume scanning video disdrometer targeting one of a typical snowfall event in Hokkaido, Japan. For the evaluation, downscaling experiments using three bulk microphysical schemes were conducted: a two-moment bulk scheme, a one-moment bulk scheme, and an improved version of the one-moment bulk scheme coupled with an identical dynamical core. The two-moment scheme successfully simulated the measured relationship between particle size and terminal velocity distribution (PSVD). On the other hand, the one-moment scheme overestimated the graupel frequency, its terminal velocity, and underestimated the particle diameter. The improved version of the one-moment bulk scheme reduced the overestimation of the terminal velocity of hydrometeors, but still overestimated the graupel frequency. We improved the overestimation of terminal velocity and the frequency of graupel in the one-moment scheme by modifying the assumption of the velocity–diameter relationship and the intercept parameter of the size distribution of graupel based on the results of the new disdrometeor. The observation of the new disdrometer would give us hints to improve the microphysics schemes in snowfall cases. (Citation: Kondo, M., Y. Sato, M. Inatsu, and Y. Katsuyama, 2021: Evaluation of cloud microphysical schemes for winter snowfall events in Hokkaido: A case study of snowfall by winter monsoon. SOLA, 17, 74−80, doi:10.2151/sola.2021-012.)


Introduction
Cloud covers broad areas of the Earth and plays important roles in its energy and water cycles. In mixed-phase and ice-phase clouds, a variety of solid particles are formed through microphysical processes such as deposition, aggregation, and riming, depending on the atmospheric conditions (Bailey and Hallett 2004;Kikuchi et al. 2013); thus, the cloud microphysical properties tend to be complex. The microphysics of solid particles are commonly characterized by mass-diameter (m-D) and terminal velocity-diameter (v-D) relationships, as well as the density and size distribution (Locatelli and Hobbs 1974;Mitchell 1996). The representation of the complex properties of solid hydrometeors in numerical models has been highly simplified in general (Tapiador et al. 2019) and remains as one of the major challenges in accurately representing mixed-phase and ice-phase clouds (Morrison et al. 2020).
Bulk microphysical schemes (e.g., Kessler 1969;Lin et al. 1983) are widely used to simulate cloud microphysics in globalscale models with explicit representation of cloud microphysics (Tomita et al. 2005;Satoh et al. 2014;Chern et al. 2016) and weather prediction models (Ferrier 1994;Morrison et al. 2005;Seifert and Beheng 2006;Saito et al. 2006; Thompson et al. 2008;Milbrandt et al. 2012). In these bulk microphysical schemes, v-D relationships are prescribed, although these relationships show large variability in the atmosphere (Locatelli and Hobbs 1974). As such, the prescribed v-D relationships are major factors in the uncertainty in cloud microphysical schemes, as reported in numerous studies (e.g., Woods et al. 2007;Iguchi et al. 2012).
This has prompted a number of studies that validated v, D and their relationships (i.e., v-D relationships) assumed by various models through comparisons with in situ observations by aircraft and ground-based measurements (e.g., Molthan and Colle 2012;Molthan et al. 2016). In particular, ground-based measurement by disdrometers such as the Particle Size Velocity disdrometer (PARSIVEL; Löffler-Mang and Joss 2000), the Laser Precipitation Monitor (LPM; Angulo-Martínez et al. 2018), and the two-dimensional video disdrometer (2DVD; Kruger and Krajewski 2002) have largely contributed to the evaluation of the v-D relationship of microphysical schemes. However, Katsuyama and Inatsu (2020) pointed out that the commonly used disdrometer introduces bias in the size-velocity distribution for solid precipitation measurements, that is rarely addressed, because the disdrometer does not observe the number concentration of precipitation particles directly. In addition, these disdrometers are usually expensive. Recently, Katsuyama and Inatsu (2021) developed a low-cost, portable volume scanning video disdrometer (VSVD) which can directly observe the number concentration, based on the work by Muramoto and Shiina (1989); the VSVD showed better performance in sizevelocity observations than the conventional disdrometers. This implies that evaluations focusing on the v-D relationship using the results of the new disdrometer provide knowledge to improve cloud microphysical schemes.
Here, we evaluated bulk cloud microphysical schemes using measurement data from the new disdrometer targeting a snowfall event in Hokkaido, Japan, in an attempt to improve the representation of cloud microphysics schemes.

Model
The meteorological model used in this study was Scalable Computing for Advanced Library and Environment (SCALE; Nishizawa et al. 2015;Sato et al. 2015) version 5.3.3. The twomoment bulk microphysical scheme of Seiki and Nakajima (2014;SN14), the one-moment bulk scheme of Tomita (2008;T08), and the improved version of one-moment bulk scheme of Roh and Satoh (2014;RS14) were used. In this study, we propose a onemoment bulk scheme based on RS14, with a modified v-D relationship and intercept parameter for the size distribution formed by graupel, as shown in Table 1. The modification of the v-D relationship assumes that v increases gradually with D compared

Comparison of the microphysical properties between the model and observations
We compared the observed PSVD and the frequency of terminal velocity with those obtained from simulations accumulated over 3 × 3 grids around the nearest grid point to the observations in longitude-latitude grid space. The simulated data were from the same period as the measurement data.
To create PSVD from the model, mean D (D ̅ ) was calculated from the model as: where D i and n(D i ) are the diameter of i-th bin in unit of mm, the particle number of i-th bin. Using D ̅ , mean v (v ̅ ) was calculated according to v-D relationship assumed in each scheme. The frequency distribution of v used for the evaluation is given by where v k , f (v k ), and n(v k ) are v ̅ of k-th bin, the frequency distribution of the k-th bin, and particle number in k-th bin. j 0.5 and j 6.0 is bin number of D of 0.5 mm and 6.0 mm, which is minimum and maximum value of the disdrometer measurement range approximately supported by VSVD. The width of bin for D and v was same as that of measurement; 0.5 mm and 0.2 m s −1 , respectively. n(D i ) in the model was calculated based on the assumption of size to RS14. The reduction of the intercept parameter assumes that graupel is diagnosed as a low number concentration compared to RS14. The turbulence and radiation processes were calculated using the Mellor-Yamada Nakanishi-Niino level 2.5 scheme (Naka nishi and Niino 2004) and MSTRN-X (Sekiguchi and Nakajima 2008), respectively. The bucket land model and single-layer urban canopy model (Kusaka et al. 2001) were used for calculating temperature and soil moisture in non-urban and urban areas, respectively. The bulk coefficient of the surface flux was determined based on Beljaars and Holslag (1991) and Wilson (2001).

Experimental setup
Simulations were conducted 33 hours from 0000 Coordinated Universal Time (UTC) on 12 February 2018, when snowfall occurred during a north-west winter monsoon over Hokkaido. The snowfall by the winter monsoon commonly occurs off the west coast of Japan including Hokkaido (Magono et al. 1966). Thus, the knowledge obtained from this case can be applied for most of the snowfall events in winter over Japan. The calculation domain covered an area of 900 × 900 km 2 , with a grid interval of 1 km, as shown by the gray-shaded area in Fig. 1. The model top height was 19,528 m, and the model layer was divided into 57 layers. The layer thickness was increased from 40 m (at the model bottom) to 650.5 m (at the model top). The results of the model were output every 30 minutes. The pressure, horizontal wind field, temperature, specific humidity, mixing ratio of hydrometeors, soil moisture, and skin temperature of a 3-hourly Mesoscale ANLysis (MANL) provided by the Japan Meteorological Agency (JMA: JMA 2019) were used for the initial and boundary conditions. RSTAR-7 Tanaka 1986, 1988) implemented in Joint Simulator for Satellite Sensors (Joint-Simulator: Hashino et al. 2013) was used to calculate the 0.64 µm radiance and the 13.3-µm equivalent blackbody temperature (TBB) from the output of SCALE. For the calculation of Joint-Simulator, the information about the size distribution of the hydrometeor is required. The size distribution assumed in each scheme was also used in Joint-Simulator. For 1-moment schemes, the mono-disperse distributions with median diameter of 16 µm for cloud water and 80 µm for cloud ice were adopted.

Observational data
Observational data from the Himawari-8 geostationary satellite (Bessho et al. 2016;Yamamoto et al. 2020;Takenaka et al. 2020), Automated Meteorological Data Acquisition System (AMeDAS) operated by JMA, the radar/raingauge-analyzed precipitation product produced by Japan Meteorological Agency (JMA) (Makihara et al. 1996), and VSVDs (Katsuyama and Inatsu 2021) installed at Sapporo (43.0830°N, 141.3367°E, ◯ point of Fig.1) and Asahikawa (43.7870°N, 142.3465°E, △ point of Fig. 1) were used for the evaluation of the model. The 0.64 µm radiance and 13.3 µm TBB of Himawari-8 was used to confirm the validity of the horizontal cloud distribution and the cloud top height simulated by the model. The total amount of precipitation from 1500 UTC on 12 February to 0900 UTC on 13 February 2018 observed at 220 AMeDAS sites (shown as × points in Fig. 1) and radar/raingauge-analyzed precipitation were used to evaluate the simulated precipitation amount. Disdrometer data was used to evaluate the number of hydrometeor particles, particle size, and terminal velocity distribution (PSVD) from Asahikawa and Sapporo from 2100 UTC 12 to 0000 UTC 13 and from 1500 UTC 12 to 0930 UTC 13 February 2018, respectively. 4.0 × 10 6 (same value of T08) Lump graupel (Locatelli and Hobbs 1974) Lump graupel (Locatelli and Hobbs 1974) 4.0 × 10 8 Graupel (Grabowski 1998) Lump graupel (Locatelli and Hobbs 1974) Fig. 1. Geographical map of the study's target area. Shading shows the domain for model calculations. Yellow and red box is the analysis area for the comparison with Himawari-8 satellite data and that of vertical mass profile among models, respectively. Cross marks denote the Automated Meteorological Data Acquisition System (AMeDAS) stations in Hokkaido operated by the Japan Meteorological Agency (JMA). Circles and triangles represent disdrometer observation sites at Sapporo and at Asahikawa installed by the authors, respectively. distribution of hydrometeor's particles in each scheme. The sum of the number concentrations of snow and graupel was used as simulated n(D i ), because only a few amounts of cloud ice were reached to ground due to its small v, and cloud water and rain were not simulated due to the low temperature near the ground. In contrast, the measured values were accumulated during the same time interval of the model output; 30 minutes.

General characteristics of simulated snow clouds
First, we validated the simulation through comparisons of the radiance, TBB, and precipitation amount with those obtained by Himawari-8 and AMeDAS. Figures 2a and 3a indicate the  horizontal distribution of 0.64 µm radiance and 13.3 µm TBB by Himawari-8. Cloud streaks formed by the winter monsoon were observed over the Sea of Japan and Hokkaido. The horizontal distribution of simulated 0.64 µm radiance and TBB indicated that the three schemes simulated the cloud streaks, regardless of the cloud microphysical scheme (Figs. 2b−d, and 3b−d), although the cloud coverage differed among models. In addition, most of the bias values of simulated precipitation amount over AMeDAS sites were smaller than 0.5 mm, and the distribution of the surface precipitation simulated by the model was similar to that by radar/ raingauge-analyzed precipitation product (Fig. S1), regardless of the cloud microphysical scheme. These results show that the three schemes reasonably simulated snow clouds and surface precipitation generated by them during the winter monsoon over land.
The reasonable representations of the cloud streaks and precipitation prompted us to conduct detailed analyses about the differences in the microphysics among individual models. SN14 well reproduced the horizontal distribution of clouds and cloud top height, based on a comparison with the observed 0.64 µm radiance (Figs. 2a and 2b) and 13.3 µm TBB (Figs. 3a and 3b), respectively. By contrast, cloud coverage by T08 was smaller than that observed; therefore, the simulated radiance (TBB) by T08 was smaller (higher) than the observed radiance (Figs. 2c and 3c). These trends were also evident in the results of RS14 (Figs. 2d and 3d). Due to the underestimation of cloud coverage, T08 and RS14 showed a smaller radiance and a higher TBB than that of observations (Figs. 2f and 3f). These trends in the simulated radiance and TBB were also seen in outgoing solar radiation (OSR) and outgoing longwave radiation (OLR) respectively (Fig. S2). As these results are calculated by the same dynamical core (SCALE), these differences originated mainly from the microphysics in the models. Next, we discuss the reason for the difference among the results of three microphysical schemes based on a comparison between the model and in situ measurements by the disdrometers.

Microphysical properties
The observed PSVD and frequency of v observed by the disdrometers (Fig. 4) enabled further clarification of the microphysical properties. The v value of most of the measured particles was ~0.8 m s −1 , and the median D was ~1.5 mm. In addition, the black box plot was close to the dashed line and the observed PSVD was mainly distributed along the dashed line, which represents the v-D relationship of snow in Figs. 4a, 4c, and 4e (also Table 2). Thus, most of the measured particles here corresponded to snow in the model.
In SN14, the simulated PSVD was close to the observed PSVD, in which v was well reproduced (Fig. 4a). In addition, the frequency distribution of v indicated that snow was dominant at the surface (blue dotted line in Fig. 4b), similar to the results observed from in situ measurements (black line in Fig. 4b). Thus, SN14 reasonably simulated the v-D relationship of solid hydrometeors measured at the surface. As well as the surface, snow was dominant in the upper layer, based on the vertical profile of the mixing ratios of graupel (q g ), snow (q s ), and total hydrometeors (q hyd ), which corresponds to the sum of cloud water, rain, cloud ice, snow, and graupel, as shown in Fig. 5.
By contrast, the PSVD by T08 (Fig. 4c) indicated that D and v was underestimated and overestimated, respectively. The frequency distribution of T08 revealed a higher graupel frequency, in which v exceeded 3 m s −1 (red dotted line in Fig. 4d), although graupel was not measured. The reason of the higher frequency of the graupel was clarified in the conversion rate from snow to graupel by the accretion of snow by graupel (P acc,s ; blue line of Fig. 5e). P acc,s in T08 was larger than that of SN14. The large P acc,s in T08 resulted in the higher frequency of the graupel with fast terminal velocity, and overestimation of v. Due to the overestimation of v, hydrometeors in the atmosphere tended to fall faster to the ground by sedimentation (figure not shown); therefore, q hyd , q s , and q g of T08 were smaller than those of SN14 (Figs. 5a−5c). The small q hyd resulted in the small cloud coverage of T08 (Figs.  2c and 3c).
The overestimation of v was reduced in the results of RS14 Fig. 4. (Left panels) Comparisons of particle size-velocity distribution (PSVD) between disdrometer observations and simulations with cloud microphysics schemes of (a) SN14, (c) T08, (e) standard RS14, or (g) modified RS14. The shaded area of (a) indicates the frequency of the observed PSVD and two-dimensional box-whisker plots indicate quartiles and the range between the 1-tile and 99-tile of the projected data onto (horizontal direction) diameter (mm) and (vertical direction) terminal velocity (m s −1 ). The shaded area and black plot are based on all observed data by the disdrometer at Asahikawa and at Sapporo throughout the observation period; the blue plot is based on all simulated snow data at the nearest gridpoint to each observational site and eight gridpoints surrounding the nearest gridpoint. The red plot is based on all simulated graupel data. The green plot is based on sum of all simulated snow and graupel data. The solid and dashed lines represent the v-D relationships of graupel and snow, respectively, assumed in a cloud microphysics scheme. (Right panels) Comparison of normalized frequency between (black) disdrometer observations and (green) simulations with (b) SN14, (d) T08, (f) standard RS14 and (h) modified RS14. The simulated normalized frequency is separated into (blue dots) snow and (red dots) graupel. The v data was divided into 20 bins from 0 to 4 m s −1 using a uniform width.
( Figs. 4e and 4f). The reduction of overestimated v was originated from the large interception parameter of RS14 (Table 1), which corresponds to the assumption that graupel has high number concentration, smaller size, and smaller v. In addition, as the RS14 turned off the accretion from snow to graupel (Roh and Satoh 2014), P acc,s was also reduced in RS14 (Fig. 5e). Due to the small P acc,s , the vertical profiles of q hyd , q s , and q g by RS14 were similar to those of SN14 (Figs. 5a−5c), in which q s was larger than q g . However, the overestimation of v and the graupel frequency, and the underestimation of D, were not still improved using RS14 (Fig. 4f); the overestimations of v and the graupel frequency also suppressed the broadening of cloud coverage, as with T08 (Figs. 2d and 3d).

Further improvement of RS14 based on in situ measurement results
To improve RS14 further, SN14, which showed good performance, was used as a reference value. Fig. 5b indicate that q g profile of SN14 and RS14 was similar to each other, but the number concentration of graupel in RS14 was higher than that of SN14 (Fig. 5d). Based on the results, the high frequency and fast v of graupel had to be modified. For these modifications, the intercept parameter of the size distribution and the v-D relationship of graupel were changed as listed in Table 1. The modified RS14 simulated lower number concentration of graupel than that by RS14 (Fig. 5d). As a result, the overestimated frequency and fast v of graupel by RS14 were improved (Fig. 4h), although the cloud coverage remained small (Figs. 2e and 3e) and D was still underestimated (Fig. 4g). These results suggest that the assumptions pertaining to graupel, such as the intercept parameter and the v-D relationship, offer room for some variability. The cloud coverage was not improved by the modification of the graupel parameter (Figs. 2 and 3). For further improvement of the 1-moment scheme, the assumption of the size distribution of snow as well as graupel is required to be modified as we discussed in supplementary material ( Fig. S3 and Text S1). Fig. 5. Vertical profile of the mass mixing ratios of (a) total hydrometeors q hyd , (b), graupel q g , (c) snow q s , (d) number concentration of graupel and snow, and (e) half-an-hourly conversion rate to graupel of unit q g (P acc,g ) averaged over the domain enclosed by 141.518°E−143.007°E and 43.115°N−45.0071°N at 0000 UTC on 13 February 2018. Dash-dotted, dotted, dash, and solid lines indicate the results of SN14, T08, standard RS14, and modified RS14, respectively. In (d), red and blue lines indicate number concentration of graupel and snow, respectively, and the number concentration was averaged over the grids in which mass mixing ratio of graupel and snow exceeded 10 −5 g kg −1 . In (e), red, blue, black, and purple lines indicate conversion rate of unit graupel q g by deposition of graupel, accretion of snow by graupel, that of cloud ice by graupel, and that of liquid hydrometeor by graupel. Table 2. v-D and m-D relationship for snow and graupel assumed by each cloud microphysical model.

Conclusion and discussion
This study evaluated three microphysical schemes through a comparison with observations, including in situ measurements by the new disdrometer; VSVDs targeting snowfall events in Hokkaido, Japan. The comparison indicated that SN14 well-reproduced the characteristics of the observations. On the other hand, T08 overestimated v due to the overestimation of graupel. The overestimation of v was small in RS14 compared to that in T08. An experiment using RS14 with modifications indicated that the overestimation of graupel frequency was improved by changing the intercept parameter and the v-D relationship of graupel.
However, further improvements with respect to the underestimation of D and small cloud coverage are required, and we examined the problem by the RS14 scheme with considering the growth process of graupel. In general, graupel grows through the riming, the deposition, and the accretion of snow and cloud ice by graupel (Knight and Knight 1973;Harimaya 1976). Figure 5e showed that the dominant process of graupel growth was different between RS14 and SN14. In SN14, the riming, the deposition, and the accretion of snow and cloud ice by graupel were all contributed as the graupel producing process. By contrast, RS14 represented only the two graupel-producing process; the riming (dashed purple line in Fig. 5e) and the deposition (dashed red line in Fig. 5e) by turning off the accretion of snow and cloud ice by graupel. Based on Harimaya (1976), who reported that the graupel contained snow by the accretion in winter Hokkaido, we imply that the accretion process is need to be included into the model. As well as the accretion, a more accurate representation of the growth process of graupel is required for simulating solid hydrometeors coexisting snow and graupel such as sea-effect snow cloud. Yasushi Fuji yoshi, Dr. Masayuki Kawashima, and Ms. Seika Tanji of Hokkaido University for their support for the disdrometeor observation, and Dr. Tempei Hashino of Kochi University of Technology for his comments about Joint-simulator. AMeDAS and Himawari-8 data were supplied by JMA and Center for Environmental Remote Sensing, Chiba University. SCALE and Joint-Simulator was developed by Team-SCALE of RIKEN, and the Earth Clouds, Aerosols and Radiation Explorer project of the Japan Aerospace Exploration Agency. Figure S1. Geographical distribution of precipitation by model and radar/raingauge-analyzed precipitation product. Figure S2. Geographical distribution of outgoing shortwave/ longwave radiation flux. Figure S3. Geographical distribution of effective radius of snow.

Supplementary materials
Text S1. Discussion about the further improvement of the 1-moment scheme