Evaluation of WRF microphysics schemes in the simulation of a squall line over IRAN using radar and reanalysis data Evaluación de esquemas de microfísica WRF en la simulación de una línea de turbonada sobre IRAN utilizando datos de radar y de reanálisis

A squall line was recorded in Dayyer port over southwest of Iran, on 19 Mar 2017. In the present paper, we have simulated the characteristic features associated with the squall line by Weather Research and Forecasting (WRF) model using five different microphysics (MP) schemes. For validating the simulated characteristics of the squall line, the latitude-height and longitude-height cross section reflectivity and precipitation value derived from observed reflectivity gathered by Doppler Weather Radar at Bushehr, synoptic weather station data at Dayyer port along with NCEP-NCAR and ERA-INTERIM reanalyzes data were used. To verify the simulated precipitation, the Fractions Skill Score (FSS) curve was calculated. Examining the simulation results for geopotential and sea level pressure show that the model simulations using different MP schemes, agree well with the verifying reanalyzes. Also, the spatial rainfall distribution of simulations and verifying observations did not show big differences. However, there are significant differences in the details of simulations such as the maximum reflectivity of the convective cells, vertical extent of the storm cells, speed and direction of the wind, rainfall values and FSS curves. Though, all of the simulations have shown convective cells over Dayyer port at the time of occurrence of the squall line, but, only the model simulation using Lin MP scheme is consistent with the corresponding radar reflectivity and vertical extent. The FSS chart showed that the skill changes with spatial scale. Results using Lin microphysics scheme crossed the FSSuniform line at lower scales when compared to other MP schemes.


INTRODUCCIÓN
Mesoscale convective systems that are linear or quasi-linear can generate heavy rains, hail, detrimental winds, and even sometimes tornadoes, called squall line (Meng et al. 2012). Radar images can show development, movement (Wong and Yip 2006) and arc shape (Meng et al. 2011) of squall lines. The bowing of squall lines is commonly associated with swaths of detrimental winds (Fujita 1978). Prediction of thunderstorms and squall lines are especially important to space vehicle launch operations, aviation, electricity services, to name a few. The methods for predicting thunderstorms can be arranged into two categories (Wilson et al. 1998). One method is a historical treatment of thunderstorm extrapolation techniques, first assuming no change in motion, size and intensity and second allowing for changes in size and intensity based on past trends (Wilson et al., 1998).
The second technique is consisted of using the numerical weather prediction (NWP) models. Prediction of squall line and thunderstorm is one of the hardest jobs in weather prediction, owing to small spatial and temporal scales and the non-linearity of their physics and dynamics. The insufficient treatment of sub-grid convection is broadly believed to be a main barrier for improving the poor efficiency of NWP models in precipitation forecasting (Liu and Moncrieff 2007). Commonly, for evaluation the performance of NWP models, the simulated thermodynamic or/and kinematic fields are compared against the corresponding observations (Jankov et al. 2010). For example, the parameters for comparison may include pressure, temperature, surface winds and precipitation values. Simulations of precipitation with high resolution 3-6 km grid spacing using convection configuration showed good agreement (Done et al. 2004;Trier et al. 2006;Liu et al. 2006;Moncrieff and Liu 2006;Lean et al. 2008;). There are various causes for better performance with the high resolution models. Also, many studies have shown that decreasing the horizontal grid spacing of the model may raise the model's capability to simulate thunderstorm and precipitation (McQueen et al. 1995;Katzfey 1995;Martin 1996;Doyle 1997;Colle et al. 1999;Davis and Carr 2000;Adlerman and Droegemeier 2002;Petch et al. 2002;Bryan et al. 2003;Kain at el. 2006Kain at el. , 2008Xue and Martin 2006;Schwartz at el. 2009;Weisman et al. 2008). Done et al. (2004) run the model with two configurations, one of them with 10-km and the other one with 4-km grid spacing. Their research showed that forecasts with 4-km resolution showed better results when compared with those with 10-km.  run WRF model to generate 5 and 7 day forecasts with the same boundary and initial conditions, same physical configurations, but different horizontal resolutions. The results showed that horizontal resolution of 2 km captures more details. Therefore, we used model output with 3 km grid spacing for Check the results.
Cloud microphysical operations play a significant role via influences on the latent heating (owing to condensation) and cold pool stability (owing to rainfall evaporation) (Rajeevan et al. 2010). Hence, microphysical parameterizations could be a main basis of uncertainty in convection in the predictions of NWP model simulations. Therefore, an important concern in convection models is microphysics parameterization (MP). Various researches have been conducted to reveal the cloud microphysics sensitivity of the models with in the prediction of thunderstorms, squall lines and associated precipitations over different region (e.g. McCumber 1991;Reisner et al. 1998;Gilmore et al. 2004;Naegele 2014;Tan 2016;Shrestha et al. 2017;Stergiou et al. 2017;Chawla et al. 2018;Eltahan and Magooda 2018;Gboode et al. 2018). Rajeevan et al. (2010) examined sensitivity of the WRF cloud microphysics to simulations of a severe thunderstorm event over southeast India. They simulated the thunderstorm using four different MP schemes. However, all the schemes underestimated stability and vertical extent of the updraft cores. Also, the microphysics schemes displayed problems in the downdrafts. While the Thompson scheme simulated rainfall closer to the actual rainfall, the other three MP schemes overestimated rainfall. Tao et al. (2011) investigated the impact of microphysics schemes on Katrina hurricane. In general, they found that microphysics schemes do not have an impact on simulated storm track but do have a major effect on the simulated intensity. Song and Sohn (2018) evaluated the WRF microphysical schemes for the simulation of heavy rain over the Korean peninsula, and found that the WRF Double Moment 6-class (WDM6) scheme was the best scheme for simulation of heavy rainfall. Our work is similar to the works conducted in the above mentioned references with some differences in the methods used for verification of the results including the results of simulation of vertical profile of maximum reflectivity, simulation of 24h precipitation accumulation, simulation of wind speed and direction and Fractions skill score.
New methods for verifying precipitation have been presented in recent years. These methods provide fractions skill score. Procedures have been explained by Ebert and McBride (2000), Casati et al. (2004), Davis et al. (2006). Robert and Lean (2008) have used radar data for precipitation accumulations for convective event to find the forecast skill.
In this research, we have simulated a squall line and seiche event that occurred over Dayyer port and Jam stations located in southwest of Iran on 19 Mar 2017 and searched the sensitivity of the WRF simulations to various cloud microphysics schemes: Lin, WSM6, Morrison, Thompson and Thompson aerosol-aware. In general, our work is very much similar to Tan, (2016); Shrestha et al. (2017) and Chawla et al. (2018), but the main purpose of the current paper is to investigate the ability of the current operational modeling system used in I.R. of Iran Meteorological Organization (IRIMO) for the prediction of squall lines which frequently occur over the southern Iran. It is to be mentioned that prediction of squall line and associated weather is a challenging task and sometimes are missed in the prediction. For example, the selected case in this study was not well predicted by the forecasting center at the time. Therefore, we want to research predictability of the squall lines and the model sensitivity to various microphysics. In Section 2, data and methodology are described. The Section 3 the main results of the simulations are presented and the conclusions are drawn in Section 4.

DATA AND METHODOLOGY
In this study, we use the Advanced Research version of WRF (ARW), version 3.9, which is a compressible, nonhydrostatic and scalar-conserving state-of-the-art atmospheric model (Skamarock et al. 2005) for the simulations of the Squall line and seiche events associated with the thunderstorm observed over Dayyer port located in southwest of IRAN on 19 Mar 2017. For model simulations, we have considered a configuration with three nested domains of 27 km, 9 km and 3 km grid spacing. The domain configurations used for model simulations is shown in Figure 1. All runs were initialized at 18:00 UTC 18 Mar 2017 and Global Forecast System (GFS) data with 0.5degree horizontal resolution were used for initial and boundary conditions. To evaluate the sensitivity of the model simulations (forecast) to cloud microphysics, five microphysical schemes were tested. Microphysics in the WRF model includes cloud, water vapor and precipitation processes. The schemes considered are Lin (Chen and Sun 2002), WSM6 (Hong et al. 2004), Thompson (Thompson et al. 2004), Morrison (Morrison et al. 2009) and Thompson aerosol aware (Thompson and Eidhammer 2014) schemes. All five schemes used, divide condensed water into cloud liquid, snow, cloud ice, rain, and graupel.
The Lin scheme is according to Lin et al. (1983), Rutledge and Hobbs (1984) and Chen and Sun (2002), modification is based on Tao et al. (1989) for ice sedimentation and saturation adjustment. The WSM6 scheme is based on Tao et al. (1989), with a different accretion calculation (Hong and Lim 2006). Split time is applied to the melting and freezing processes for increasing accuracy in the vertical heating profile. The saturation adjustment follows Dudhia (1989) and Hong et al. (1998) in separately processing ice and water saturation. In the Thompson scheme ( Thompson et al. 2004) a variable has been added to predict the concentration of cloud ice. In this scheme ice nucleation and autoconversion is calculated according to Cooper (1986) and Walko et al. (1995) respectively. For graupel category represents a gamma function. All hydrometeors in the Morrison's double-moment scheme Pinto 2005, 2006) except for cloud water are double moment.
One of the factors affecting cloud microphysics is Aerosol chemistry (rajeevan et al. 2010). Khain et al. (2005) investigated the effect of aerosols on dynamics and cloud microphysics and found a significant effect. Clouds arising under aerosol conditions produce a powerful downdrafts and strong convergence in the boundary layer. As such, aerosols can contribute to the formation of convective cells and thunderstorms. Being triggered by dynamical forcing, clouds arising in microphysically air are stronger and can produce a squall line. Because of this, we will also review a version of Thompson aerosol aware scheme (Thompson and Eidhammer 2014), that was carried out into the WRF model from 2014, and considers water-and ice-friendly aerosols.
The Yonsei University (YSU) PBL scheme (Troen and Mahrt 1986) was used to parametrize the PBL processes. For longwave and shortwave radiation, the Rapid Radiative Transfer Model (Mlawer et al. 1997) and Dudhia scheme (Dudhia 1989), were used respectively. All simulations in the ARW-WRF model were performed with 30 vertical levels. For all simulations and in the first and second domains the Betts Miller Janjic (BMJ) Cumulus parameterization scheme (Janjic 1994) was used while the convection for inner domain (3 km grid spacing) was considered zero. The first 12 hour of the simulation was considered as the spin-up time. Therefore, for discussions, we have used only the results from 00:00UTC of 19 Mar 2017 onward. All runs were initialized at 1200 UTC 18 Mar 2018 and Global Forecast System (GFS) data with 0.5degree horizontal resolution (forecast NCEP GFS data) was used for initial and boundary conditions.
The main objectives of this paper are to investigate if the WRF model with high resolution is capable of capturing the observed characteristics of the squall line, to evaluate the reactivity of simulation to varying microphysics and also to investigate predictability of the squall line by the WRF model. To this end, the vertical profile of maximum reflectivity associated with the squall line, spatial distribution of 24-h precipitation accumulation and surface wind (speed and direction), were examined first. Then the fractions skill score for precipitation verification was applied. To perform the fractions skill score method, the base reflectivity data for elevation angle 0.5° (Brandes et al. 1999) from Bushehr Doppler Weather Radar over southwest Iran and the Marshal-Palmer equation for calculating precipitation are used to estimate the forecast and radar data in a same grid. Four thresholds of 0.5, 5, 10 and 15 mm and 95 th percentile threshold (as threshold to represent the heaviest precipitation in the system) are considered and used to convert the forecast and radar precipitation into binary (I O and I M ) and the fractions skill score (FSS) were calculated, according to: Where, the reference MSE (MSE(n)ref) and MSE(n) are given by, (2) Where Nx and Ny are the number of columns and rows in the domain, O(n)i,j and M(n)i,j are observed and forecast fractions obtained from IO and IM. Then, the FSS curve is plotted for neighborhood length (n) (Refer to Robert and Lean 2008).
Also, for verifying simulated characteristics of the squall line, the synoptic data at Dayyer port and JAM stations and Doppler Weather Radar (DWR) data were used.

Details of the squall line
The squall line convective phenomenon and the seiche event occurred in the morning of 19 March 2017 over southwest of Iran and were recorded in Dayyer port and Jam synoptic meteorological stations. shows the spatial and time distribution of maximum reflectivity observed by the DWR at Bushehr station which represents the genesis, growth and extension of the squall line event. The start of the event was observed over west side of the Persian Gulf around 0000 UTC 19 March 2017, then the convective cells grew and moved east (Figure 2a, b). At 0300 UTC (Figure 2c), in the center of the Persian Gulf, they formed a squall line. This squall line could be seen around Dayyer station at 0330 UTC ( Figure 2d) and 0400 UTC ( Figure 2e) and also, could be seen around Jam station at 0430 UTC (Figure 2f). The squall line started dissipating around 0500 UTC. Due to the passage of the squall line, strong wind and cold air subsidence took place over the Persian Gulf that generated a standing wave and seiche event over the Dayyer port and Jam.     Figure 6d with the strong negative anomaly located in the region of the trough. The contour map of precipitatble water (Figure 6e) shows a region of deep humidity with a value more than 50 (kgm −2 ) over southwest of Iran. The contour map of the precipitable water anomaly shows at the maximum positive anomaly of 15 (kgm −2 ) which is again located in the region of low pressure area (Figure 6f). Results of our past experience in this region (Arkian and Karimkhani;2014) show that the value of 9 kgm −2 for precipitate water is a critical value, above which the probability of precipitation is high.

Comparison of simulated geo potential fields with reanalysis data
To verify the model results, the simulated fields including 500-hPa geopotential and mean sea level pressure are compared against the verifying ERA-INTERIM reanalysis dataset (Jankov et al. 2010). The simulated results for the outer domain (27 km resolution) along with the verifying reanalysis contour map for 500hPa geopotential heights are presented in Figure 7. The show results in Figure 7b are for model run with the Lin scheme as microphysics. It is to be noted that the main features of the synoptic patterns shown in 500-hPa and sea level pressure maps are almost the same for all five aforementioned microphysics schemes used. As seen in Figure 7b the simulated 500-hPa geopotential shows a through with axis in the southeast direction over southwest of Iran and agree well with the verifying reanalysis ( Figure 7a). The same results hold for sea level pressure (Fig is not presented). and lower values for the Morrison scheme. It is seen that the highest reflectivity occurs at the altitude 0f 13000 m. This indicates that there are stronger updrafts for the Lin scheme at 04:00, which is consistent with the corresponding radar reflectivity (  Figure 5b). Also the pattern of the multicell convective systems simulated by the model using Lin MP scheme is closer to the radar reflectivity ( Figure 5b).

Simulation of 24-h precipitation accumulation
24-h precipitation accumulations (mm) at 06UTC 18 to 19 Mar 2017 simulated by the model over the inner nest with five difference MP physics are presented in Figure 9 For five simulations with different MP schemes, the simulated spatial distribution of the precipitation amounts over southwest of IRAN was consistent with the corresponding observed precipitation accumulations (Figure 4). Though, the spatial distribution and areal coverage for five simulated patterns are similar but they differ from each other in terms of value. As such, the highest and lowest values of precipitation were simulated using Lin and Thompson schemes respectively.

Fractions skill score of precipitation accumulations
Due to the complexity of rainfall patterns, their prediction is difficult, but when the visual assessment of the precipitation forecast is similar to the observation, FSS is also consistent with this result (Robert and Lean 2008). To evaluate the simulation results quantitatively, the FSS method of Robert and Lean were used. Figure 11 shows the observed and simulated 12-h accumulated rainfall from 0000 to 1200 UTC 19 Mar 2017. Fractions skill score charts are shown in Figure 12 for all the five forecasts using different MP schemes. To check the spatial accuracy of the forecasts, the 95th percentile threshold (about 14 mm as threshold to represent the heaviest precipitation in the system) was selected. The precipitation forecast for Morrison and Thompson aerosol-aware MP schemes showed relatively poor results in terms of precipitation area (Figure 11d, e). Figure 12 supports this view by showing that the forecasts using Morrison and Thompson aerosol-aware MP schemes are less skillful compared to a random forecast up to 450 and 250 km scales respectively.
The WRF model with Lin MP scheme predicted the spatial rainfall distribution well (Figure 12b), while the maximum rainfall value using WSM6 MP Scheme match well with the verifying observation ( Figure  12c). FSS diagram also shows that for the Lin and WSM6 MP Schemes, the results for scales between 160 km and ~120 km were better compared to the random forecasts.  Chart of FSS for accumulation rainfall thresholds of 0.5, 5, 10 and 15 mm are displayed in Figure 13. For the smallest threshold (0.5 mm) the FSS values, using five different microphysics schemes, are greater than FSSrandom and are significantly skillful over all scales (Figure 13(a)). For 5.0-mm threshold, FSS curves using Lin and WSM6 MP Schemes crossed the FSSuniform line at 235 and 215 km respectively. While FSS curves using other three different schemes doesn't cross the FSSuniform line (Figure 13(b)).
Generally, the FSS curves using Morrison and Thompson aerosol-aware microphysics schemes don't cross the FSSuniform except for the lowest threshold of 0.5 mm and thus show weaker performance compared to other schemes. For 10 mm threshold (Figure 13c), the FSS curves using Lin, WSM6 and Thompson MP schemes, cross the FSSuniform at 165, 175 and 205 km respectively. Two forecasts using Lin and WSM6 MP Schemes have the same FSS values for 15 mm threshold (Figure 13d).

CONCLUSIONS
In this paper, a squall line observed over Dayyer port (Iran) on 19 Mar 2017, was simulated using WRF model with 3 km grid spacing. The aim was to examine if the WRF model is capable to simulate the main characteristics of the squall line and whether different microphysics schemes influence to squall line simulation. For this purpose, five different microphysics schemes in the WRF model (Lin, WSM6, Morrison, Thompson and Thompson aerosol aware) were selected. The simulations showed similar results in terms of occurrence of the squall line. In general, the model using all microphysics relatively faithfully simulated the synoptic scale patterns of geopotential and sea level pressure over the southwest of Iran. Both simulated and reanalyzed fields of 500hPa geopotential height map show a trough over the southwest of Iran. Also, the model simulated the formation of the convective cells as close to the verifying observations. All the simulations show the passage of convective cells over the southwest of Iran, as observed by Doppler Weather Radar. Moreover, the spatial rainfall distribution is almost the same in all the simulations.
In general, it is found that the WRF model was capable to simulate many large scale characteristics of the squall line successfully. But examining the high resolutions reveals some significant difference between the simulations and the verifying observations. Significant differences are observed in the simulations of the number of convective cells, the vertical extent of convective cells, horizontal wind speed and direction. Also, the FSS curve for different simulations differed from each other. In brief, among the five MP schemes considered, the Lin scheme simulations were closer to observation. These results are derived from the simulation of the single squall line. We propose to do simulations for more number of squall line observed over Iran to generalize the results obtained. This paper has found the problems associated with the simulation of characteristics of the squall line and its sensitivity to five different MP schemes.
It is interesting to attend that the Thompson scheme originally designed to improve the mid-latitude winter precipitation (Rajeevan et al. 2010). The Morrison scheme is a double moment scheme, in which the concentration of cloud ice is treated explicitly and the Thompson aerosol-aware MP scheme considers water-and ice-friendly aerosols. However, the above mentioned MP schemes failed to simulate the features associated with the squall line successfully. It is important to understand why the microphysics schemes at the WRF model had a problem in simulating the vertical extent and strength of convective cells associated with the squall line. Rajeeavan et al. (2010) stated that the reasons could be associated to the sensitivity of cumulus parameterization schemes, differences in the simulations of graupel, changes in land surface properties and initial conditions.