Mountain Waves Analysis in the Vicinity of the Madrid-Barajas Airport Using the WRF Model

Department of Earth Physics and Astrophysics, Complutense University of Madrid, Madrid, Spain High Resolution Limited Area Model Consortium (HIRLAM), Madrid, Spain State Meteorological Agency (AEMET), Santander, Spain Institute of Interdisciplinary Mathematics (IMI), Complutense University of Madrid, Madrid, Spain Atmospheric Physics Group, IMA, University of León, León, Spain Department of Applied Mathematics, University of Valladolid, Valladolid, Spain


Introduction
Severe weather conditions, such as hail, heavy precipitation, lightning, mountain waves, icing, wind shear, or turbulence, can affect aircraft safety [1][2][3]. Particularly, turbulence and icing episodes can promote loss of control of the aircraft, which is one of the main causes of aviation accidents [3,4]. According to the National Transportation Safety Board [5], turbulence accounts for 71% of all the 446 weather-related accidents in commercial aviation in the USA for the period of 2000-2011. Concerning mountain wave turbulence, the number of accidents recorded in the USA was 68 (7, 9% of the total accidents) and 113 (13, 1% of the total) were classified as clear air turbulence (CAT) during 1987-2008 [6]. In Europe, for the period of 2012-2016, 2991 incidents were related to turbulence and 312 others to in-flight icing, with four of them resulting in fatal accidents [3]. Both phenomena, turbulence and icing, are associated with mountain wave episodes, rendering this meteorological phenomenon a considerable risk to aviation. ere are two principal atmospheric conditions required for the mountain wave generation: a strong wind flow perpendicular to an orographic barrier and a stable environment [7][8][9]. As wind is forced to climb over the barrier, if the air parcel reaching over the top encounters a slightly stable layer, it will induce a gravity wave. us, mountain waves are formed downwind and over the inducing orographic barrier. Depending on the conditions, this wave can propagate vertically and horizontally, reaching several kilometers downwind [10]. Several authors have researched the different factors by which the turbulent flow associated with mountain waves can occur at both high and low levels in the atmosphere [9,[11][12][13][14][15]. e gravity wave is necessarily associated with an alternating vertical wind speed. Depending on the frequency, this alternation can generate in-flight turbulence for an aircraft flying through the mountain wave. Bolgiani et al. [4] study a mountain wave event on the leeward side of the Guadarrama mountain range, relating a moderate-to-strong turbulence observation to a simulated vertical wind speed above 2 ms −1 . e associated turbulence is more usual during the winter season (especially during December and January in the northern hemisphere), as winds are stronger, and the atmosphere is more stable [16]. Furthermore, when there are no clouds evidencing the mountain wave, the turbulence then becomes CAT, often associated with wind shear [6]. Wind shear conditions can produce asymmetric changes in the aircraft's lifting force, which is most critical during take-off and landing operations.
In addition to turbulence, mountain waves can also produce aircraft icing. Condensation is promoted in the updrafts by collision and coalescence processes [17] which produce an increase in the liquid water content (LWC), occasionally forming cloud bands associated with the mountain wave's phases [18]. When these wave clouds develop in areas with temperatures below 0°C, LWC in the cloud can become supercooled as the process of nucleation is less efficient than condensation [19]. If these supercooled liquid droplets (SLD) impact the aircraft's surface, icing can be produced. Due to the nucleation process efficiency, SLD prevail at temperatures above −10°C, rendering icing more likely to occur between −2°C and −15°C [20]. Pilot reports register 50% of aircraft icing events between −8°C and −12°C [21] and between 1,500 and 4,000 meters above sea level (masl). Ledesma and Baleriola [2] established a temperature of −40°C as the limit where supercooled drops can be found, even without ice nuclei presence (through homogeneous nucleation processes).
Considering the LWC, Tafferner et al. [22] establish several grades of aircraft icing: LWC lower than 0.1 g kg −1 sustained for one hour does not affect the aircraft's safety, values between 0.1 and 0.5 g kg −1 can produce light icing, and LWC between 0.5 and 1.0 g kg −1 produce moderate icing. Contents greater than 1 g kg −1 can generate severe aircraft icing, meaning that the anti-icing systems equipped by the aircraft cannot cope with the amount of contamination on the wing surface. A particularly dangerous type of icing can happen when temperature is approximately −2°C, as SLD instead of freezing immediately flow over the wing slowly freezing on it. is is known as "clear ice," which is hard to remove, since it is located behind the anti-icing system. On the contrary, if this process occurs at lower temperatures, the drops will immediately freeze over leading edges (wings, propeller. . .) producing "rime ice" [2]. is is less dangerous for the aircraft safety since anti-icing systems are located in these leading edges. Both types of ice increase the weight of the aircraft and change the lift and thrust profiles. Moreover, clear ice can detach the airflow causing turbulence downstream and loss of control [23].
Several forecasting products are available for aviation users nowadays. e World Area Forecast Centers produce, among others, operational turbulence and icing forecasts tailored for aviation [24]. e Significant Weather Chart (SIGWX) and the Significant Meteorological Information (SIGMET) are adverse meteorology information reports dedicated to aircraft in a specific area. Another product is the Graphical Turbulence Guidance (GTG), which produces a forecast by using a weighted combination of several different turbulence diagnostics [25]. Nevertheless, these products are far from perfect and require updating and the use of newer and more precise methods [4,26].
Forecasting mountain waves is not an easy process because of the conjunction of several atmospheric factors such as the synoptic pattern, wind speed and direction, and stability profile. In addition, the phenomenon is strongly influenced by orography; therefore, the use of mesoscale numerical prediction models is necessary to produce an accurate icing forecast. In this regard, the weather research and forecasting (WRF) model has been proven valuable to simulate mountain wave events in several studies [4,15,27,28]. De la Torre et al. [14] evaluated the frequency of mountain wave episodes located near and over the range peaks of the Andes Mountains taking into account different thresholds of downdrafts and updrafts velocity. Furthermore, in the last years, other authors have carried out sensitivity analyses of LWC using several parametrizations of the WRF model [29]. Finally, the use of satellite remote sensing has been proven useful for nowcasting mountain waves and icing conditions [4,[30][31][32].
One example of the evidence of the risk generated by mountain waves to aviation and the need of improving the forecasting can be found in the case of flight IRC3704. is flight operated in an ATR72 aircraft encountered turbulence and icing associated with a mountain wave en route to Yasuj (Iran) in February 2018. As a result, from both phenomena, the crew lost control of the aircraft and impacted the terrain resulting in 66 casualties and the total destruction of the airframe [26]. e SIGMET and SIGWX reports underestimated the mountain wave and icing conditions [21,33] failing in the purpose of warning the pilots about the unsafe conditions. is paper is performed within the framework of the SAFEFLIGHT research project for the improvement of aviation safety related to icing and weather hazards. e objectives of this study are to evaluate if the WRF model can reproduce mountain wave episodes observed via the Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager (MSG-SEVIRI) products. Using high-resolution simulations, thirteen observed mountain waves episodes affecting the Adolfo Suarez Madrid-Barajas Airport (LEMD hereafter, as per the International Civil Aviation Organisation code) are analysed. ese are located south of the Guadarrama mountain range with prevailing northern winds and are liable to generate a hazard to the arrival and departure operations of the airport [2,8]. Several microphysics parameterizations schemes are used to simulate the selected events in order to assess the best parametrization or the need to use an ensemble for forecasting. In addition, the principal atmospheric variables conducting to mountain waves are evaluated. Finally, a case study is presented in order to analyse if the WRF is a useful tool to forecast mountain wave episodes which cannot be validated with the observation of wave clouds.
is paper is organized as follows: Section 2 describes the experimental design, including the study domain, the WRF model setup, and the MSG-SEVIRI data used. Section 3 explains the methodology used to produce and validate the results. e results and discussion (Section 4) are divided in four subsections: initial analysis, study case of mountain wave event, brightness temperature (BT) analysis, and the main features associated with mountain waves. Finally, the major conclusions are summarized in Section 5.

Study Domains.
e spatial domain for this study is designed to cover the area where mountain waves affecting LEMD are usually observed. LEMD is the busiest airport in Spain and the sixth in Europe, with more than 55 million passengers during 2018. e airport elevation is 610 masl and it is located on a plateau in the centre of the Iberian Peninsula, approximately 40 km southeast of the Guadarrama mountain range (Figure 1(b)). is mountain range, which is part of the Central System, runs from northeast to southwest with an approximate length of 80 km, with Peñalara being the highest elevation at 2, 428 masl. When northwesterly winds hit this range, mountain waves can be formed to the south. us, the spatial domain selected covers the windward side of the mountains and the leeward side beyond the location of LEMD, making a 121 × 121 km domain.

WRF.
e WRF numerical model (version 3.5.1) [34] is used for simulating the selected events.
is is a nonhydrostatic model extensively used and validated for weather research and forecasting. Initial and boundary conditions are taken from the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) operational analysis, with a horizontal spatial resolution of 0.25°and a temporal resolution of 3 hours [35]. Each episode is independently simulated in periods of 24 hours. Simulations are initialized at 00:00 UTC, allowing for the first 6 hours as spin-up time.
ree geographical domains are defined, following a two-way nesting strategy. Each domain has 121 grid points in north-south and east-west directions and 40 sigma vertical levels. Figure 1(a) depicts the outer domain (d01) including the Iberian Peninsula with 9 km of spatial resolution, the middle domain (d02) localized in the centre of Spain with 3 km resolution, and the inner domain (d03) with 1 km of spatial resolution covering the Guadarrama mountain range and LEMD.
As different physics options allow optimizing the WRF simulations depending on the meteorological event and resolution to evaluate [34], several combinations of physics parametrization schemes are used to study the selected episodes. Particularly, different microphysics, Planetary Boundary Layer (PBL), surface layer, land surface, and radiation parametrizations are combined, creating 20 experiments (Table 1). e main differences between the PBL and microphysics schemes chosen are as follows.

PBL Schemes
(i) Yonsei University (YSU) [36]: this scheme intensifies the boundary layer mixing in the thermally induced free convection regime, while reducing it in the mechanically induced forced convection. It has been tested in other studies related to cloudiness [37] and icing [38]. (ii) Mellor-Yamada-Janjic (MYJ) [39]: this scheme is optimized for deep convective regimes assuming a new parameter called cloud efficiency. It has been proven to be a good performer in similar BT [40] and icing studies [38]. (iii) Mellor-Yamada-Nakanishi-Niino (MYNN) [41]: this parametrization considers the effects of buoyancy on pressure and stability on the turbulent length scale. It has been used in mountain waves studies over the same area [4,28]. (iv) University of Washington (TKE) [42]: it involves a new moist turbulence parametrization. It is characterized by the use of moist-conserved variables and a new diagnosis and computation of turbulent kinetic energy. is PBL scheme is suitable in study areas with complex terrain [42].

Microphysics Schemes
(i) Goddard (GOD) [43]: this scheme is optimized to simulate condensation and evaporation processes; it has been proven to be a good performer in aircraft icing studies in the Iberian Peninsula, although not specifically for mountain waves [38,44]. (ii) ompson (THO) [45]: this microphysics scheme is defined to improve the forecast of aircraft icing, since it considers mixed-phase processes relative to another bulk microphysics. It has been used by Otkin et al. [46], Bormann et al. [40], and Montejo [37] for the validation of simulated clouds using BT. (iii) Morrison 2-moment (MOR) [47]: even if this microphysics scheme is not optimized for icing conditions, it has been proven to be a good performer in similar studies [27,38].

Advances in Meteorology
(iv) WRF single-moment 6-class (WSM6) [48]: this scheme is developed as an improvement of the default scheme optimized for WRF. It has been tested in other studies related to cloudiness [40] and icing [38]. (v) WRF double-moment 6-class (WDM6) [49]: this microphysics scheme based on WSM6 scheme includes a prognostic variable of cloud condensation nuclei number concentration. It is commonly used in studies about great convective events like cyclones, storms, and heavy precipitation.

MSG-SEVIRI.
e MSG-SEVIRI products are used as observational data to validate the WRF simulations through the presence of cloud bands associated with mountain waves. e MSG is a geostationary satellite operated by the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT). It is composed by 12 spectral wavelength channels; two visible channels (0.6 and 0.8 µm), a high-resolution visible (HRVIS) channel and a near-infrared channel (1.6 µm); and eight infrared (IR) channels (from 3.9 to 13.4 µm). Within these, there are two watervapor channels (6.2 and 7.3 µm). e horizontal resolution is 3 × 3 km at nadir, except for the HRVIS channel which offers images with a horizontal resolution of 1 × 1 km. e temporal resolution is one hour. Of all the channels available, only HRVIS and 10.8 µm IR are used, for the reasons explained in the following section.

Selection of Episodes and Synoptic Situation.
e mountain waves episodes for evaluation are selected using MSG-SEVIRI HRVIS images between 2017 and 2019. Aircraft icing conditions are more likely during winter, due to the melting level being lower; thus only the periods from November to March are used. e events are considered if two conditions are met: north-northwest winds are observed in the low troposphere over the study area, and wave cloud bands are observed on the leeward side of the Guadarrama mountain range. Using this method, 53 events are observed. From this selection, only the best observations are chosen for analysis, considering high-level clouds concealing the event, how well defined the wave is, and how far southeast the wave clouds propagate, giving priority to those events reaching LEMD. In total, 13 mountain wave episodes are selected: In order to understand the mountain waves formation, the synoptic situations are analysed first. e events selected follow very similar meteorological conditions at synoptic scale.
In 10 of the 13 events, the mid (500 hPa) and lower (700 hPa) troposphere configuration at synoptic scale was dominated by a conjunction between the Azores Anticyclone and a cyclone or relative low over central Europe  Advances in Meteorology ( Figure 2), thus creating a trough over the western Mediterranean Sea. is promotes north and northwesterly flows at 500 hPa, often accompanied by cold air advection (maritime or continental) over the Iberian Peninsula. At lower levels, a flow perpendicular to the Guadarrama mountain range can be observed. Similar synoptic configurations are described in studies of mountain waves in the study area [4,28]. e situation for 02 December 2017 is governed by an omega block, with a cut-off low east of the Iberian Peninsula and a high west of Ireland, producing similar flows to those of the previous cases. On 12 March and 06 November 2018, the synoptic configuration is dominated by a large and strong Icelandic Cyclone, which drives south the Atlantic anticyclone generating westerly flows over the Iberian Peninsula.   e MSG-SEVIRI images are also used to validate the results of the WRF simulations. As there are no observational data for the variables aloft for these events, wave clouds are used as an objective variable intrinsically related to the phenomenon. For similar studies, Otkin et al. [46] and Bormann et al. [40] already used BT to validate simulated WRF data. Accordingly, in this paper, BT is used to compare the observed and simulated wave clouds images. Considering the resolution used in the simulations, it would have been desirable to use the observation HRVIS channel and d03 data for validation. However, WRF does not generate a BT product equivalent to an HRVIS image, but it does for IR pseudosatellite images using Unified Postprocessing (UPP) System developed at the NCEP, which uses a couple radiative transfer model to compute the derived brightness temperatures. In this line, according to Otkin et al. [46], the 8.7 µm, 10.8 µm, and 12.0 µm IR channels are the most appropriate to detect cloud top and surface temperatures. Even more, 10.8 µm and 12.0 µm are coincident with the IR atmospheric "windows" central wavelength [37,40,50], making these channels particularly sensitive to the presence of clouds. us, the longwave IR channel (10.8 µm) is selected as observational data to evaluate the simulations. e assessment is performed from 09: 00 to 15:00 UTC, as during this time period wave clouds are present for every selected event.
As the MSG-SEVIRI IR resolution is 3 km, it would have been desirable to validate against d02 simulations. Nevertheless, an initial evaluation shows that this domain does not properly capture wave clouds from every event selected, while d03 does. e satellite images are cropped to the same spatial domain of d03 and the pseudo-IR-images from d03 are interpolated and regridded to upscale the resolution from 1 km to 3 km.
Consequently, both the satellite images and the pseudo-IR-images that are postprocessed match with the same domain and resolution. is allows analysing the BT from WRF and MSG-SEVIR at the same resolution (3 km) and same grid points number. To use the same temporal domain as the observations, the simulation hourly results from 09:00 to 15:00 UTC are taken, making 7 daily timesteps and a total of 91 timesteps for the events selected. With these adjustments, the observed and simulated images can be compared. In addition, to avoid any possible noise generated by the orographic clouds in the mountain wave evaluation, the final domain of study is produced by removing the windward side (in both products), as depicted in Figure 1 Having established a proper domain, an initial analysis is carried out. e 20 experiments are used to simulate three randomly chosen mountain waves events (26 January 2018, 26 November 2018, and 27 January 2019). Bormann et al. [40] describe the importance of the BT frequency distribution to establish the realism of the images, particularly in terms of the general distribution of clouds. Following their methodology, the frequency distribution of BT is plotted for every experiment and the observations. e results allow us to discard some of the experiments. For the remaining, further validation examination is performed.
where BT s is the simulated BT at a specific grid point for every time step and BT o is the observed BT at a specific grid point for every time step.
where N is the total number of time steps considered. (3) e results allow us to select the two best performing experiments. Only these are used for the simulation of the 13 events and further validation.

Spatial Validation, Main Features, and Case Study.
e best performing experiments are then fully evaluated. e BT frequency distribution is generated, and the aforementioned skill scores are averaged over the 91 timesteps for each grid point. us, spatial distributions of each skill score can be plotted. To summarize this information into a single value, averages of the spatial patterns are also computed.
Once the optimum experiment is established, the relevant atmospheric variables for mountain wave development are studied for every selected event. To establish an altitude for considering these variables, the LWC is analysed, finding its maximum at approximately 2,500 masl. is altitude is used to obtain the data for relative humidity (RH), LWC, wind direction (WD), wind velocity (WV), and static stability (ST), as these govern the mountain wave and related icing. For this purpose, three grid points are established: the Navacerrada mountain pass (N) and the leeward (L) and windward (W) sides (Figure 1(b)). e latter two are located 20 km from Navacerrada along a line perpendicular to the direction of the mountain range and aligned with LEMD. WD, WV, ST, and RH are used to obtain the most frequent values. A frequency distribution is generated for each variable at each study point and a summary table is computed  for 10, 20, 50, 80, and 90 percentiles. Finally, a case study is examined in order to evaluate the turbulence and aircraft icing conditions associated with mountain waves. Pseudo-IR-images, LWC, and vertical wind speed are horizontally plotted for that purpose. In addition, temperature and LWC are plotted as a cross section following the plane defined by W, N, L, and LEMD points. e same case is also evaluated at a time when no wave cloud is present to assess the ability of WRF to detect mountain waves when no visual evidence exists.

Results and Discussion
is section is structured as follows: first, the initial analysis is shown in order to elucidate the best parametrization schemes. en, a case study is presented, followed by the brightness temperature analysis. Finally, in the last subsection, the main features of mountain waves can be found. is order allows a better understanding and evaluation of the experiment, as the reader will have a previous example to visualize.

Initial Analysis.
e initial analysis is carried out with the aim to evaluate the ability of the 20 experiments to simulate cloud bands associated with mountain waves, as these wave clouds will be later used to assess the model. ree events are randomly chosen among the 13 days previously selected. Figure 3 depicts the BT frequency distribution of the 20 WRF's experiments. Two different behaviours can be clearly distinguished. e experiments using the THO and WSM6 microphysics are able to simulate BT lower than 260 K. On the other hand, the MOR, GOD, and WDM6 microphysics simulate BT only greater than 260 K.
ese results seem to be unrealistic, as it must be considered that these experiments are generating no clouds.
To verify this, the LWC is plotted and evaluated at each timestep (not shown) finding that, in effect, the WDM6, GOD, and MOR experiments do not simulate any clouds in the domain of study for the 3 episodes evaluated. is incompetence may be due to the fact that these microphysics schemes are developed to optimize the formation and evolution of convective systems and trailing stratiform precipitation [43,47,49]. As wave clouds are not connected with convection, present shallow vertical developments, and rarely produce precipitation, these schemes are not adequate for their simulation. Anyhow, this feature immediately invalidates these parametrizations for this study. Only 8 experiments (those using THO and WSM6 microphysics) are suitable to forecast wave cloud with the WRF model. e skill scores for the 8 aforementioned experiments are shown in Table 2. Differences can be appreciated between the PBL schemes. Experiments with YSU and TKE schemes obtain better RMSE (lower than 6.40) and R (greater than 0.50) than the MYJ and MYNN schemes. BIAS is negative in all experiments, which reveals that the BT is systematically underestimated. is underestimation is more pronounced in the MYNN schemes, while the best score (−2.31) is obtained by YSU-THO and TKE-THO experiments. e averages of BTfor every experiment are slightly lower than the MSG-SEVIRI average. Even if TKE presents better SD scores, an overall assessment makes YSU the most suitable PBL scheme to simulate wave clouds. ese results may be originated in the fact that, while the MYNN, MYJ, and TKE parametrization considers the effects of buoyancy on pressure, deep convective regimes, and stability on the turbulent length scale, the YSU PBL scheme intensifies the boundary layer mixing in the thermally induced free convection regime, while reducing it in the mechanically induced forced convection [36,39,41,42]. YSU considers the mechanical forced convection, involved in the generation of mountain waves. e results also reveal that the microphysics parameterization is more relevant than the PBL parameterization for the simulation of the selected events.
In consequence, and attending computational cost, only the YSU-THO and YSU-WSM6 experiments are used for the full analysis of this paper. e final configuration of the WRF uses Dudhia shortwave scheme [54] and RRTM longwave scheme [55] for radiation, Unified Noah [56] land surface scheme, revised MM5 [57] surface layer scheme, and Yonsei University (YSU) [36] PBL scheme, all combined with the microphysics schemes ompson [45] and WSM6 [48]. Henceforth, the experiments are named according to the microphysics scheme used ( ompson and WSM6).

Case Study of Mountain Wave Events.
To evaluate the ability of the WRF simulations to reproduce the mountain waves, a case study is presented, namely, the ompson experiment for 26 November 2018. is is selected as per the results of the microphysics validation, presented in Section 4.3. e atmospheric conditions at 13:00 UTC for this day favoured the generation of mountain waves. A WV of 14 m s −1 , WD of north-northwest, and a temperature value of −6°C are simulated on W at 2,500 m. e stability is 0.0025 K Pa −1 over N at 2,500 masl. e temperature is −3.7°C and the RH is 80% over L. Two times are selected to evaluate the turbulence and the aircraft icing conditions: at 13:00 UTC when wave clouds are simulated and observed in the MSG-SEVIRI product and at 22:00 UTC when no clouds evidenced the wave.
In the first place, the BT images in the 10.8 µm band for MSG-SEVIRI and for WRF are displayed in Figure 4 at both times evaluated. e observations (Figures 4(a) and 4(c)) present a poor quality due to the 3 km spatial resolution. Nevertheless, wave clouds are clearly present at 13:00 UTC, reaching beyond LEMD location with BT approximately 265 K on cloud tops. At least six phases can be differentiated with a decrease in the cloud extension as the wave propagates southeast. At 22:00 UTC no wave clouds are appreciated. It is worth noting that at this time (night local time) the visual channels of MSG-SEVIRI are not useful, which is a clear disadvantage in the use of this product for wave forecasting.
e pseudosatellite images (Figure 4(b) and 4(d)) depict a similar pattern at both times, with a much better quality as per the 1 km resolution. e wave clouds are Advances in Meteorology clearly simulated at 13:00 UTC, reaching LEMD and presenting the same number of phases although the cloud extension may be smaller than observations. BT on cloud tops is very close to the observations. At 22:00 UTC no clouds are represented in the leeward side in accordance with observations.
With the objective of evaluating if the WRF is able to capture mountain waves when no visual evidence exists, LWC and vertical wind speed are assessed. ese two variables are highly dependent on mountain waves, although vertical wind speed will always be present and LWC may not. Figure 5 shows several simulated outputs obtained for d03. As expected, LWC (Figures 5(a) and 5(d)) confirms the same results from the pseudosatellite images (Figures 4(b) and 4(d)). When the cross section is analysed at 13:00 UTC ( Figure 5(b)), it is evident that the LWC is creating the wave clouds according to the wave phases. It can also be appreciated that the wave is affecting the temperature in the lower atmosphere up to 5,000 masl. is is evidence that mountain waves disturb the troposphere much beyond the extent of the orographic barrier. Another important result is that LWC is almost completely present at temperature below 0°C; thus, it can be considered that these are SLD. Moreover, the maximum LWC is located at a height of 2,500 masl approximately, coincident with the observations by Bolgiani et al. [4]. e cross section at 22:00 UTC ( Figure 5(e)) presents no LWC but still depicts the disturbance of the temperature produced by the wave, showing weaker oscillations. Vertical wind speed results (Figures 5(c) and 5(f )) show several alternating downdraft/updraft bands. Vertical speeds of ±2 ms −1 are observed, resulting from the mountain waves. At 13:00 UTC a stronger wave can be seen extending along the complete leeward domain and presenting larger wavelengths. At 22:00 UTC the event is weaker, and the pattern is dissipating; nevertheless, the mountain wave is still evident and reaches LEMD.
As stated by the existent literature, turbulence and icing conditions are related to mountain waves [4,22,28]. us, a risk of moderate aircraft icing (LWC above 0.5 g kg −1 ) and moderate turbulence (vertical speed ± 2 m s −1 ) can be assumed at 13:00 UTC as well as risk of low-to-moderate CAT at 22:00 UTC over the LEMD airport area.

Brightness Temperature Analysis.
e previous results show that the WRF model can properly capture mountain wave events and the associated wave clouds. Nevertheless, an objective physical variable is required to validate the parametrizations. As LWC and vertical wind speed observations are not available, the analysis of BT is performed as described in Section 3. Figure 6 shows the BT frequency   [46] and Bormann et al. [40] who obtain similar curves using the same IR channel and the ompson microphysics. However, both experiments overestimate slightly the frequencies associated with the mode value, presenting leptokurtic distributions with respect to the observed frequency distribution.
To perform a more specific BT analysis, the skill scores are evaluated. e MSG-SEVIRI average BT (Figure 7(a)) presents spatial distribution with a clear influence of the orography. Low BT are observed over the higher elevations and conversely high BT are observed over the lower lands. A similar pattern can be appreciated in the WSM6 and ompson results (Figures 7(b) and 7(c), respectively), although the temperature differences found in the respective domains are lower than the observations. e spatial average BT for WSM6 (267 K) and ompson (268 K) present results of 5 K and 4 K below the observations, respectively. Regarding the spatial distribution of SD, important differences can be found. e observations (Figure 7(d)   produces larger variations throughout the temporal domain in the southwest corner of the domain; in consequence, it can be considered that this area has the most variable cloudiness. e WSM6 and ompson experiments (Figures 7(e) and 7(f ), respectively) present the largest variations over the same area, although they generate much larger SD than observed. is is an evidence of the model propagating clouds in a much more variable way than observed. Both experiments overestimate the averaged SD, with values of 14.1 K (WSM6) and 15.3 K ( ompson), while the observation is 9.9 K.
RMSE results (Figures 8(a) and 8(b)) are in accordance with the SD results. As WSM6 and ompson (Figures 7(e) and 7(f ), respectively) produce large dispersion in their data, the error is also high. e averaged RMSE is 15.27 K for the ompson experiment and 15.57 K for the WSM6 experiment. e spatial distribution of the average BIAS presents differences depending on the zone over study area. Near the Guadarrama mountain range, where the wave clouds should be more defined and present during a longer time, BIAS is close to 0 (Figures 8(c) and 8(d)). However, BIAS is smaller to the south of the domain, which means that BT is underestimated in this side. Also, the BIAS spatial distribution of both experiments is quite similar. Finally, when R is considered, substantial results are found. While the WSM6 experiment yields averaged results of 0.30, the ompson experiment clearly outperforms it, with an average value of 0.51 (statistically significant at a 99% level), which is a considerable correlation with the observations. R spatial distributions show that the worst performer is WSM6 (Figure 8(e)), presenting correlations between 0 and 0.1 in the northeast of the domain. On the contrary, the ompson experiment (Figure 8(f )) depicts the best R values over the wave propagation area, which is consistent with the more detailed ice-phase and other cold-phase processes of this scheme in comparison with the default microphysics scheme (WSM6) optimized for WRF [45,48,58]. In Figures 8(a) and 8(b) it is noteworthy that in the SW corner there is an area with high RMSE which is coincident with high values of R, shown in Figures 8(e) and 8(f ). is fact is due to a large error in this area which constant over time, and consequently the correlation is high. Finally, it is noteworthy that these skill scores are very different from those of the initial analysis, which may be due to the large variability of this type of atmospheric phenomenon. However, the YSU-THO configuration remains the best parametrization in both the initial analysis and the total one.
Considering the frequency distributions ( Figure 6) and the skill score results (Figures 7 and 8), the ompson microphysics can be considered the best experiment for simulating wave-induced clouds in the study domain, in line with the scheme design by ompson et al. [45] and previous result by Bolgiani et al. [4]. Even if SD shows relatively poor results, the performances presented by the ompson simulations in RMSE, BIAS, R, and frequency distribution of BT make it the best choice possible. An additional conclusion that can be obtained from these results is that the development of an ensemble to forecast these mountain wave events may add only marginal information to the ompson simulations. Only a small improvement may be achieved by incorporating other  experiments analysed here. e potential robustness that could be gained by an ensemble does not compensate the computational costs of these simulations.

Main Features of Mountain Waves.
Once the ompson experiment has been established as the best option for simulating mountain waves in the study area, the next objective of this paper is to summarize the values of several atmospheric variables that can induce these events. As on initial evaluation the highest simulated LWC values are found at 2,500 masl (e.g., Figure 5(b)), this level is set for the analysis of the variables at the three selected points of evaluation (see Figure 1(b)). Figure 9(a) shows the prevailing WD and WV at 2,500 masl for each point of analysis. When considering the wind variables, it is worth noting that the most important value for mountain wave formation is recorded at W, as the wave largely depends on the flow reaching the barrier. In addition, the orography will perturbate the wind over N and L, rendering these values less reliable. e data over W shows a WD range from 315°to 005°and a typical WV between 12 and 18 m s −1 .
is is consistent with the prevailing WD (from 315°to 023°) climatology data for the proximate surface station located in Navacerrada [59]. When reaching N, the WD range narrows and the WV increases, most probably due to orographic effects. In almost 35% of the time, the prevailing WD is 340°, reaching WV between 18 and 24 m s −1 (Figure 9(b)).
is is consistent since N is located over a mountain pass where the wind flow narrows, increasing the WV by Venturi effect. At L, the WD spreads again, ranging from 300°to 360°, with WV similar to the W data (Figure 9(c)). e wind observations (21.1 m s −1 of WV and west-northwest WD) by Bolgiani et al. [4] on the same area are within the wind rose range for L and the climatology for LEMD airport surface station also presents a westnorthwest mean direction during the winter [60]. Overall, the prevailing WD in the three points analysed is northwest to north for the selected events, particularly in W point vicinity, where WD matches with other studies in winter months [61,62]. Table 3 presents a distribution summary of the main variables considered. Regarding the WD and WV values, the previous results are confirmed. e prevailing direction of north-northwest (between 315°and 005°) can be considered the characteristic WD for the mountain wave events evaluated. In addition, it is evident that a minimum WV of 12 m s −1 is required in 80% of the cases for the event to appear, with typical values under 17 m s −1 . Evaluating the probability density function for this variable (Figure 10(a)), a characteristic WV around 15 m s −1 (≈P50) can be established.
Considering the ST results, the most important value is taken at N, as the variable may be affected by an orographic dipole (formation of meso-high on W and a meso-low on L). Table 3 presents a narrow distribution, especially concentrated at N, as depicted by the leptokurtic curve in Figure 10(b). e typical ST value over N is neutral, ranging from −0.013 to 0.012 K Pa −1 (P20 and P80, respectively). On W and L the data shows that the distributions tend to slightly stable conditions, consistent with the required situation for mountain wave generation [28,63].
Concerning the RH percentiles for mountain waves, priority has to be given to the values over L, as the waverelated icing conditions will only be present on the leeward side. P80 values are 100% over W and N (Table 3), which is consistent with the orographic clouds over the barrier associated with mountain waves, most noticeable in the density function for N (Figure 10(c)). e data over  L presents a more evenly distributed variable, generating the lowest RH value for P20 (42.3%) and a median value of 65%. Based on these results, the percentiles can be used to define a characterization of mountain waves and icing conditions. Minimum and typical values for wind over W and ST over N would define the generation of the events in the area of study. RH values associated with temperatures below 0°C would set the conditions for icing. Both characterisations would be of interest; however, they would require the examination of a much larger number of cases, which is beyond the scope of this paper.

Summary and Conclusions
Mountain wave and icing episodes are adverse meteorological phenomena that can affect aviation safety and air traffic management. In this paper, 13 mountain wave events are selected among 53 observations on the Guadarrama mountain range area from 2017 to 2019. e events are simulated using the WRF model, and an analysis with several parameterizations schemes is carried out. Five microphysics schemes are selected in this paper in order to assess the requirement of an ensemble and to choose the best parameterization schemes to Table 3: 10th, 20th, 50th, 80th, and 90th percentile values of the WD, WV, ST, and RH distributions. Data are taken over W (windward), N (Navacerrada), and L (leeward) at 2,500 masl for the total mountain wave events selected.  forecast mountain waves. Pseudosatellite images simulated by WRF are validated using the observed BT from the MSG-SEVIRI. Different scores are used to assess the skills of each selected parameterization scheme in simulating the BT. e best parametrization scheme is used to evaluate the main features in mountain waves. Finally, the WRF ability to detect turbulence and icing associated with mountain waves is evidenced in an example case study.
It is concluded that it is possible to detect mountain waves using WRF IR pseudosatellite representation of the wave clouds, as could be done using MSG-SEVIRI observations. However, only the experiments using WSM6 and ompson microphysics schemes are able to simulate the cloud bands associated with mountain waves. e ompson parametrization combined with the YSU PBL scheme is found to be the most suitable experiment among the different WRF schemes used to detect mountain waves in the vicinity of the LEMD airport. e microphysics parametrization is found to be more relevant in the mountain wave cloud simulations than the four PBL schemes evaluated. Based on 80 and 20 percentiles, wind direction between 315°-005°and wind velocity between 12 and 17 m s −1 on the windward side as well as a stability range within ±0.012 K Pa −1 in Navacerrada are the typical conditions for the selected mountain wave events. Turbulence and icing risks associated with mountain waves can be properly detected using vertical wind speeds and LWC simulated. ese WRF simulations using the ompson and YSU parameterizations can be considered a useful forecasting tool for mountain wave events south of the Guadarrama mountain range. Even when no observations can be made, the model properly represents the atmospheric wave.
In summary, the WRF model is considered for the forecasting of icing conditions and turbulence connected to mountain waves as a valid tool to improve aviation safety. In fact, even when no nowcasting can be made using satellite products, the WRF simulations can still represent the event. It would also be interesting to use more simulated events to produce a characterization of the mountain waves for the area, which can even lead to the development of a forecasting algorithm.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.