Evaluation of Microphysics Schemes in the WRF-ARW Model for Numerical Wind Forecast in José Martí International Airport †

: A sensitivity study was developed with Lin, Morrison 2-moment, weather research and forecasting (WRF) single-moment 5-class (WSM5), and WRF single-moment 6-class (WSM6) micro-physics schemes available in the weather research and forecasting-advanced research WRF (WRF-ARW) for the numerical forecast of the wind field at José Martí International Airport, in Cuba. The selection of these schemes was based on their use in numerical weather forecast systems operating in Cuba. As case studies, five storms associated with synoptic patterns that cause dangerous conditions at this aerodrome were selected. The simulations were initialized at 0000 UTC with the forecast outputs of the global forecast system (GFS) model. The schemes were evaluated according to the wind field’s representation in the region where the airport is located, the headlands, and the center of the runway. The errors observed are strongly dependent on the occurrence of convection, especially on the intensity and the factors that cause it. During the dry season (November–April), the lowest errors are observed, while the worst performance is appreciable for the rainy period (May– October). Lin and WSM6 schemes reproduce the best behavior of the wind field on the aerodrome.


Introduction
The wind field forecast is one of the most important meteorological supports for air operations.The spatial resolution of meteorological phenomena that modifies the wind field along the runway often exceeds the range of local sensors.Numerical weather prediction (NWP) models are an alternative to be used as alarm systems in aeronautics [1][2][3][4].
The weather research and forecasting (WRF) model [5] is widely used to simulate the near-surface wind for both research and operational applications.It has two dynamic cores, the advanced research WRF (ARW) and the non-hydrostatic mesoscale model (NMM), developed by the National Centre for Atmospheric Research (NCAR) and National Centers for Environmental Prediction (NCEP), respectively.Shaw et al. [2] implemented the WRF-ARW v2.2 model [6] for the Dubai International Airport aviation weather decision support system.The authors installed an operational system assimilating data from satellites, radiometers, wind profiles, radar, and surface observations.In the Hong Kong International Airport, a sub-kilometric NWP capability in capturing low-level wind shear was evaluated [4].This aviation model (AVM) [7,8] is a sub-kilometer resolution implementation of the WRF.WRF model offers multiple applications and, like most of NWPs, several physics options.For wind forecast, commonly, sensitivity boundary-layer parametrizations (PBL) studies were developed [9][10][11][12][13][14][15][16][17].In Cuba, the meteorological model WRF-ARW sensitivity to physics options was tested [18,19].The prediction systems SiSPI [20] and SPNOA [21] were developed and implemented in the Center of Atmospheric Physics of the Institute of Meteorology of Cuba, but not specific to the aviation application.For this purpose, Díaz-Zurita et al. [22] improved a numerical wind surface derived from WRF-NMM for José Martí International Airport.This airport is located near the elevations of Cacahual.The catabatic flow modifies the characteristics of the meteorological variables at the aerodrome and its vicinity [23].Furthermore, at this aerodrome, Sosa [24] points out that storms in the vicinity of the aerodrome usually originate dangerous phenomena associated with wind field variations.
Microphysics parametrizations are significant in predicting storms, as is described in many types of research [25][26][27][28][29].To provide a preliminary evaluation of numerical wind field forecast over José Martí International Airport, a sensitivity study was developed.From WRF practice recommendations [30], ARW core was used in this research.Based on results that show the influence of storms in aerodrome vicinity wind field variations, different microphysics schemes' capability to represent these changes was verified.

Study Area and Case Studies
Barcía et al. [31] divided Cuba into forecast regions according to the behavior of the meteorological variables recorded in the observational meteorological network.The regions are classified according to the extreme temperatures, the influence of the sea breeze, and the physical-geographical characteristics.José Martí International Airport is located in the inland forecast region of the Artemisa, Havana, and Mayabeque provinces (Figure 1a).This airport has a 4 km runway (Figure 1b), with a southeast-northeast orientation.The airport is surrounded by terrain with complex orography.At the aerodrome, there is a catabatic flow of moist air that favors a drop in temperatures around the runway [23].
Sosa [24] used meteorological data from observations in the José Martí International Airport during the period between 2012 and 2017 to describe the behavior of the low-level wind shear.This hazardous phenomenon for aircraft is associated with meteorological systems: cold front (18.36%), anticyclone (53.06%), and tropical wave (10.20%) [24].The author identified the following synoptic patterns as the most frequent in which low-level wind shear occurs: 1. Influence of the North Atlantic Subtropical Anticyclone with trough medium and high levels; 2. Influence of the North Atlantic Subtropical Anticyclone in the entire tropospheric column; 3. Migratory anticyclones; 4. Tropical waves into the south of western Cuba. 5. Cold fronts on western Cuba.
Sosa [24] refers to that significant wind field variations over the airport are often reported under storms.As case studies, five storms associated with synoptic patterns that cause dangerous conditions at this aerodrome were selected (Table 1).

Model and Domain Configuration
The WRF v.3.9 model was used with ARW dynamic core.The simulations comprised 12-4km two-way nested domains (Figure 2) and 34 verticals levels.Briefly, the setup includes for both domains: Rapid Radiative Transfer Model longwave radiation parametrization [32], Dudhia shortwave radiation scheme [33], Unified Noah land-surface model [34], Grell-Freitas ensemble cumulus parametrization [35], Mellor-Yamada-Janjić planetary boundary layer [36].The model was initialized and forced at the boundaries by a 0.500 GFS forecast, with every 3 hourly updated boundary conditions.

Experimental Design
Storms simulations were performed with four selected microphysics schemes: Lin [36], , WSM5 [38], and WSM6 [39].Table 2 shows the main species of prognostic variables in these schemes.The selection was based on their use and performance in numerical weather forecast systems operating in Cuba [18,19].The forecasts were for 54 h, started from the storm observation date at 0000 UTC.  1a).The model was initialized with the GFS forecast, which is freely available at https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p50.pl(accessed on 15 January 2020).On the other hand, radars products were obtained from the Key West, Florida, USA (KBYX) doppler radar (available online at https://www.ncdc.noaa.gov/nexradinv/chooseday.jsp?id=kbyx (accessed on 15 Janaury 2020)).In addition, the software NOAA Weather and Climate Toolkit v.4.5.0 (free available at https://www.ncdc.noaa.gov/wct/install.php(accessed on 12 Febraury 2020)) was utilized to analyze the radar data.
In the domain of 4 km (d02) of the resolution, the density of the nodes in the neighborhood of the airport is low.For this reason, the rectangular grid developed by Díaz-Zurita et al. [22] was used.This grid takes into account the orientation (60 from the horizon) and the length of the runway (4 km).Five points are matched: one in the center of the runway (MID), one at each headland, and the other two points at 1 km from the center.The grid with a longitudinal resolution of 0.87 km and a latitudinal resolution of 0.5 km is shown in Figure 3.
Based on the results of Díaz-Zurita et al. [22], for wind interpolation, the natural neighbor method was used.In addition, following Díaz-Zurita et al.'s [22] recommendations, a correction to the interpolated wind field is applied with a consistent mass model.This model is based on the equation of continuity for an incompressible air mass moving in a two-dimensional domain, Ω, with a velocity field  ⃗ , ,  : If the constant air density is considered for the entire domain, the equation becomes: which joins the impenetrability condition on the ground  , thus constituting the boundary condition: From conditions ( 2) and ( 3), the consistent mass models pose a least-squares problem with the velocities to adjust  ⃗ , ,  from the observed  ⃗  ,  ,  in the Ω domain, according to the functional: where { , ,  ,  , ,  , and  , ,  , } are the wind components calculated by the model through fit; { , ,  ,  , ,  , and  , ,  , } are the components of the initial field, interpolated from the observations, and  ,  ,  are the Gaussian precision modules [40].Considering  and  identical, for horizontal directions, the functional to minimize ( 4) is: The search field  ⃗ , ,  will be the solution to the problem: Find  ⃗ ∈ K such that, This problem is equivalent to finding the saddle point at ( ⃗, ) of the Lagrangian: The technique of Lagrange multipliers allows obtaining the saddle point of the expression (8),   ⃗,  ≤   ⃗,  ≤ ) such that the solution field is obtained from the Euler-Lagrange equations: where Φ is the Lagrange multiplier and  =  ,  ,  is the diagonal transmission tensor: If  ,  are considered constant throughout the domain, the variational formulation leads to an elliptic equation defined in Φ.Indeed, substituting Equation ( 8) in (2) results: which is completed by the null Dirichlet condition at permeable boundaries (vertical domain boundaries) and Neumann's condition in the raincoats (terrain and upper border) Considering  and  constant, Equation (11) becomes: eliminating the vertical component (two dimensions) was obtained: This methodology guarantees the conservation of wind direction due to the impenetrability conditions.
To evaluate the radar's basic features and structure simulation of storms, verticals profiles of reflectivity and relative humidity were generated.Relative humidity was determined using Clausius-Clapeyron [41] (Equation 1): where: • T = temperature (K); • p = pressure (Pa); • q = specific humidity or the mass mixing ratio of water vapor to total air (dimensionless); • T0 = reference temperature (typically 273.16 K) (K).

Evaluation
The model output was compared with meteorological observations.In this work, the forecast verification was calculated some statistical metrics: mean systematic error or BIAS, mean absolute error (MAE), root-mean-squared error (RMSE), and Pearson's correlation coefficient (rp).To obtain the best microphysics scheme to reproduce the surface wind properties derived from storms, the analysis was developed pre-, in-and after-the occurrence of the event [29].

Convective Storm Analysis
The late afternoon ground heat flux, the inland breeze convergence, and high-pressure levels diffluence origins the storm of 29 June 2012.Particularly, this storm caused the highest wind speeds records on the runway.At 2:00 pm local time (18:00 UTC), the storm had a maximum height of 13.91 km and a core of maximum reflectivity of 60.5 dBZ at 4.37 km.The spatial pattern of maximum reflectivity, vertical profile mixing ratio of hydrometeor particles, and postprocessed wind field are presented for the tested microphysics schemes when the storm occurred (Figure 4).WSM5 (Figure 4g) and WSM6 (Figure 4j) simulated high maximum reflectivity (50-55 dBZ) in the nearest location as radar did (22.978, −82.345).Both Lin and WSM6 schemes produced more than one core cloud echo.That behavior of WSM6 was founded previously by Sari, Pulung, and Sukma [29] for a hail event study case in Surabaya, Indonesia.
The postprocessed wind field for each scheme is presented; similar predicted wind speeds and directions were shown by WSM5 (Figure 4i) and WSM6 (Figure 4l).The authors reported that the microphysics schemes were not sensitive to surface properties (main wind flow in particular) and suggested a non-sensitivity of wind direction to microphysics parameterization [28,29].However, convective storms develop were sensitive to the microphysics schemes [42][43][44][45][46]. Figure 4 shows the differences in the postprocessed wind fields.The WSM5 and WSM6 wind fields are quite similar, which can be linked with the position and size of the simulated storm.Figure 4 also shows vertical simulated hydrometeor profiles; perceptible variations were observed.All schemes predicted mixing ratios of hydrometeors only for warm cloud processes, which is likely to occur in Cuba.WSM6 (Figure 4k) scheme shows the best radars basics features among the schemes.In the first place, WSM6 produces more graupels than the other schemes; second, WSM6 presented the greatest decrease in the vapor mixing ratio near 7 km, when the rest of the diagrams show it over 5 km; and finally, WSM6 also had the highest rain mixing ratio, which decreased with increasing height, from the surface to 6 km.

Wind Field Simulation
Figures 5 and 6 show skills metrics for forecasting wind speed and direction pre-, inand after-the occurrence of storms selected as case studies.In the rainy season storms, for the wind speed and wind direction forecasts, the highest biases are shown in the presence of storms (Table 3).In the case studies, the occurrence of several mesoscale convective cells caused the heterogeneous distribution of the wind speed and direction.In addition, during storms, the biggest differences in skills metrics for wind direction were obtained.This variable was underestimated.On the runway, the error measurements were more dispersed, both for speed (Figure 5d) and direction (Figure 6d).This may have been a consequence of the different simulated storm positions and characteristics from each tested microphysics.To consider the best microphysical scheme, this analysis remains difficult.However, WSM6 and Lin occasionally exhibit better scores.Previously research for this study area [25] report that WSM6 has a skill for wind forecast.For the dry season, the lowest biases in the wind speed forecast were observed during the occurrence of storms, as summarized in Table 4.In the case of wind direction, the highest biases are seen at this time.This behavior was similar to that observed during the rainy season.Skills metrics indicate a more accurate forecast than in the rainy season.The fundamental difference between the cases was in the origin of the storms.
In the dry season, Pearson's correlation in the region reached values of up to 0.6, as it is shown in Figure 7. On the runway, the Pearson's correlations ranged between 0.7 and 0.99 (Figure 8d), which suggests the model's ability to represent changes in wind direction, probably due to mass consistent correction applications.Once more, the schemes show fewer differences among themselves.The forecast errors may also be attributed to a lateral boundary condition, which was obtained from 0.50 × 0.50 horizontal resolution of the GFS forecast data.

Conclusions
A sensitivity study was developed with Lin, Morrison 2-moment, WSM5, and WSM6 microphysics schemes for the wind field's numerical forecast at José Martí International Airport.As case studies, five storms associated with synoptic patterns that cause dangerous conditions at this aerodrome were selected.The sensitivity of the microphysics scheme using the WRF model on 29 June 2012, was discussed.It is important to note that WSM6 showed the most realistic storm radar features and vertical profile hydrometeors, also, it represented distributions according to warm cloud processes, which are likely to occur in Cuba.In general , we obtain that the wind field was modified by the position and size of the simulated storms.
Furthermore, skills metrics pre-, in-and after-storms were obtained.Between seasons was observed pronounced differences in errors, possibly linked to dissimilarities in storms genesis conditions.Skills metrics indicated a more accurate forecast in dry season storms than in the rainy season ones.Additionally, the consistent mass correction application may cause higher correlations between wind directions forecast at the runway.
To consider the best microphysical scheme, this analysis remains difficult.However, WSM6 showed better scores in the major criterion of the study developed.Nevertheless, these results are the first attempt to obtain the best configuration of the WRF model for numerical storm forecasts in the airport.Ongoing work will, therefore, include other sensitivity meteorological field analyses.

Figure 5 .
Figure 5.Taylor diagrams for rainy season wind speed forecast.Panels from left to right inland region errors measurements and runway error measurements: (a,b) pre-; (c,d) in-; (e,f) after-storm.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Taylor diagrams for rainy season wind direction forecast, panel from left to right inland region errors measurements and runway error measurements: (a,b) pre-; (c,d) in-; (e,f) after-storm.

Table 1 .
Case studies (Synoptic patterns as the numbered list).