Abstract
Turbulence and aircraft icing associated with mountain waves are weather phenomena potentially affecting aviation safety. In this paper, these weather phenomena are analysed in the vicinity of the Adolfo Suárez Madrid-Barajas Airport (Spain). Mountain waves are formed in this area due to the proximity of the Guadarrama mountain range. Twenty different weather research and forecasting (WRF) model configurations are evaluated in an initial analysis. This shows the incompetence of some experiments to capture the phenomenon. The two experiments showing the best results are used to simulate thirteen episodes with observed mountain waves. Simulated pseudosatellite images are validated using satellite observations, and an analysis is performed through several skill scores applied to brightness temperature. Few differences are found among the different skill scores. Nevertheless, the Thompson microphysics scheme combined with the Yonsei university PBL scheme shows the best results. The simulations produced by this scheme are used to evaluate the characteristic variables of the mountain wave episodes at windward and leeward and over the mountain. The results show that north-northwest wind directions, moderate wind velocities, and neutral or slightly stable conditions are the main features for the episodes evaluated. In addition, a case study is analysed to evidence the WRF ability to properly detect turbulence and icing associated with mountain waves, even when there is no visual evidence available.
1. Introduction
Severe weather conditions, such as hail, heavy precipitation, lightning, mountain waves, icing, wind shear, or turbulence, can affect aircraft safety [1–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.
There 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–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. Thus, 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–15]. The 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. The 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. This 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]. This 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. The World Area Forecast Centers produce, among others, operational turbulence and icing forecasts tailored for aviation [24]. The 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–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. This 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]. The 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.
This paper is performed within the framework of the SAFEFLIGHT research project for the improvement of aviation safety related to icing and weather hazards. The 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. These 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.
This 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. The 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.
2. Experimental Design
2.1. Study Domains
The 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. The 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)). This 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. Thus, 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.
(a)
(b)
2.2. WRF
The WRF numerical model (version 3.5.1) [34] is used for simulating the selected events. This 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. Three 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). The main differences between the PBL and microphysics schemes chosen are as follows.
2.2.1. 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. This PBL scheme is suitable in study areas with complex terrain [42].
2.2.2. 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)Thompson (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].(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.
2.3. MSG-SEVIRI
The MSG-SEVIRI products are used as observational data to validate the WRF simulations through the presence of cloud bands associated with mountain waves. The 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 water-vapor channels (6.2 and 7.3 µm). The horizontal resolution is 3 × 3 km at nadir, except for the HRVIS channel which offers images with a horizontal resolution of 1 × 1 km. The 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.
3. Methodology
3.1. Selection of Episodes and Synoptic Situation
The 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. The 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:10 November 2017, 02 December 2017, 09 December 2017, 02 January 2018, 26 January 2018, 03 February 2018, 12 March 2018, 02 November 2018, 06 November 2018, 26 November 2018, 14 December 2018, 20 January 2019, and 27 January 2019.
In order to understand the mountain waves formation, the synoptic situations are analysed first. The 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 (Figure 2), thus creating a trough over the western Mediterranean Sea. This 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]. The 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.
3.2. Initial Analysis
The 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. Thus, the longwave IR channel (10.8 µm) is selected as observational data to evaluate the simulations. The 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. The 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. This 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(b).
Having established a proper domain, an initial analysis is carried out. The 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. The results allow us to discard some of the experiments. For the remaining, further validation examination is performed.
Based on the methodology followed by Lopez et al. [51], Loew et al. [52], and WWRP/WGNE [53], several skill scores are used for analysis.(1)BIAS where BTs is the simulated BT at a specific grid point for every time step and BTo is the observed BT at a specific grid point for every time step.(2)Root-mean-square error (RMSE) where N is the total number of time steps considered.(3)Linear product-moment correlation coefficient of pearson (R)(4)Standard deviation (SD)
The results allow us to select the two best performing experiments. Only these are used for the simulation of the 13 events and further validation.
3.3. Spatial Validation, Main Features, and Case Study
The best performing experiments are then fully evaluated. The BT frequency distribution is generated, and the aforementioned skill scores are averaged over the 91 timesteps for each grid point. Thus, 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. This 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)). The 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. The 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.
4. Results and Discussion
This section is structured as follows: first, the initial analysis is shown in order to elucidate the best parametrization schemes. Then, 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. This order allows a better understanding and evaluation of the experiment, as the reader will have a previous example to visualize.
4.1. Initial Analysis
The 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. Three 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. The 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. These 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. This 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.
The 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. This underestimation is more pronounced in the MYNN schemes, while the best score (−2.31) is obtained by YSU-THO and TKE-THO experiments. The averages of BT for 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.
These 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. The 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. The 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 Thompson [45] and WSM6 [48]. Henceforth, the experiments are named according to the microphysics scheme used (Thompson and WSM6).
4.2. 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 Thompson experiment for 26 November 2018. This is selected as per the results of the microphysics validation, presented in Section 4.3. The 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. The stability is 0.0025 K Pa−1 over N at 2,500 masl. The 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. The 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. The 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. The wave clouds are 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. These 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. This 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]. The 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]. Thus, 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.
4.3. Brightness Temperature Analysis
The 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 distribution (10.8 µm channel) of the observations and experiments for the selected events. The MSG-SEVIRI produces a curve with a long tail extending to the lower BT. The tail is a sign of the observation of clouds at several altitudes, while the platykurtic zone (270–280 K) may result from the abundance of lower clouds hiding the ground BT. The Thompson and WSM6 curves do not display such zone, although they show similar distribution and variability range of BT values. Moreover, clouds are simulated with BT lower than 265 K, proving that both experiments properly capture the clouds at various altitudes and reproducing the characteristic tail seen in the observations curve. These results are in agreement with the studies by Otkin et al. [46] and Bormann et al. [40] who obtain similar curves using the same IR channel and the Thompson 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. The 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 Thompson results (Figures 7(b) and 7(c), respectively), although the temperature differences found in the respective domains are lower than the observations. The spatial average BT for WSM6 (267 K) and Thompson (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. The observations (Figure 7(d)) present lower deviations to the northeast of the domain, while displaying higher deviations to the southwest. This suggests that the BT 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. The WSM6 and Thompson experiments (Figures 7(e) and 7(f), respectively) present the largest variations over the same area, although they generate much larger SD than observed. This 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 (Thompson), 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 Thompson (Figures 7(e) and 7(f), respectively) produce large dispersion in their data, the error is also high. The averaged RMSE is 15.27 K for the Thompson experiment and 15.57 K for the WSM6 experiment. The 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 Thompson 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 Thompson 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). This 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 Thompson microphysics can be considered the best experiment for simulating wave-induced clouds in the study domain, in line with the scheme design by Thompson et al. [45] and previous result by Bolgiani et al. [4]. Even if SD shows relatively poor results, the performances presented by the Thompson 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 Thompson simulations. Only a small improvement may be achieved by incorporating other experiments analysed here. The potential robustness that could be gained by an ensemble does not compensate the computational costs of these simulations.
4.4. Main Features of Mountain Waves
Once the Thompson 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. The data over W shows a WD range from 315° to 005° and a typical WV between 12 and 18 m s−1. This 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)). This 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)). The 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 west-northwest 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].
(a)
(b)
(c)
Table 3 presents a distribution summary of the main variables considered. Regarding the WD and WV values, the previous results are confirmed. The 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.
(a)
(b)
(c)
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). The 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 wave-related 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)). The 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.
5. 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. The 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 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. The 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 Thompson microphysics schemes are able to simulate the cloud bands associated with mountain waves. The Thompson 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. The 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. These WRF simulations using the Thompson 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.
Data Availability
Initial and boundary conditions are taken from the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) operational analysis; (National Centers for Environmental Prediction)/National Weather Service/NOAA/U.S. Department of Commerce; “Updated daily: NCEP GFS 0.25 Degree Global Forecast Grids Historical Archive”; Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, 2015, DOI: https://doi.org/10.5065/D65D8PWK (accessed 14/02/2020). Satellite data are taken from European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) (https://archive.eumetsat.int/usc/, accessed 23/07/2020).
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was partially supported by the following research projects: PID2019‐105306RB‐I00, PCIN-2014-013-C07-04, and PCIN2016-080 (UE ERA-NET Plus NEWA Project), CGL2016-78702-C2-1-R and CGL2016-78702-C2-2-R (SAFEFLIGHT Project), FEI-EU-17-16 and SPESMART and SPESVALE (ECMWF Special Projects). Special thanks go to Roberto Weigand for computer support. Javier Díaz-Fernández acknowledges the support grant from the MINECO-FPI program (BES-2017).