A Case Study on Assessing Cumulonimbus Induced Flight Vulnerabilities Over the Nepalese Himalayan Terrain

With a boom in the tourist industry both in India and adjoining Nepal, middle-class tourists from the Indian subcontinent have started routinely flying over the Himalayan terrain over the last decade. This influx of tourists stretches its aviation sector, in a country with some of the most vulnerable airports. The Himalayan terrain plays a key role in shaping weather systems over the region which has witnessed several aircraft crashes. Tourists in Nepal routinely fly 19-seater aircrafts such as the Beechcraft 1900D for a joy ride over the Himalayan range and are particularly vulnerable. This paper assesses the risks and vulnerabilities associated with such short duration flights starting from Kathmandu and covering parts of the Himalayan region. The region experiences deep cumulonimbus clouds which form over a period of a few days and are ubiquitous during the monsoon season (June to September). The vertical extent of such clouds ranges from 2 km upwards to up to the tropopause. This paper first assesses cumulonimbus mediated hazards along flight routes with thunderstorm activity and then details a critique on the hazard tied down to secondary effects i.e. deposition of supercooled droplets on a Beechcraft 1900D. In particular, the paper explores the role of cumulonimbus induced icing on the aircraft surfaces which can severely affect tourist flights over this part of the developing world.


Introduction
The Himalayan region experiences one of the most intense thunderstorm activities in the world (Zipser et al. 2006). During the monsoon season, warm moist air from the Indian ocean encounters the Himalayan range and is forced to rise up (orographic lifting). The subsequent cooling of this air results in the formation of large convective clouds initiating thunderstorm activity. Several flights within this Himalayan country, including scenic mountain flights flying across this terrain are severely affected (see Fig. 1).
These short-duration flights are here to stay, at least for the coming several decades boosting Nepal's economy. With some reduction in snow cover, this region is already experiencing the ravages of global warming (Bastola 2017). Rise in global temperatures may lead to more intense hydrological cycles which in turn would lead to an increase in the frequency of convective clouds and precipitation (Fowler and Hennessy 1995). Additionally, an influx of tourists is currently stretching Nepal's aviation sector (BBC 2019). These factors severely impact Nepal and, being a landlocked country, air travel is its main connection with the outside world.
Even without convective storms, there have been several crashes owing to localized extreme weather shaped by complex terrain. Some of these crashes have been analyzed by Regmi (2014a, b) and Regmi et al. (2017) using the National Oceanic and Atmospheric Administration's (NOAA) Weather Research and Forecast (WRF) model (Skamarock and Klemp 2008). It must be noted that the WRF and other Numerical Weather Prediction (NWP) models use a combination of available observations and short-term forecasts for model initialization purposes. The model's success in yielding accurate forecasts rests on the quality of data used for the purposes of initialization. However, one must also recognize that there is often a deficiency in data acquisition, in both ground-based observations as well as upper-air soundings, mainly in low-income and developing countries (Gumber et al. 2020;Ghosh et al. 2018;Strigaro et al. 2019). In addition, such observations are unevenly distributed in regions with complex topography such as in the Himalayas (Orr et al. 2017). One, therefore, needs to undertake detailed assessment studies in order to compensate for the data assimilation bias that may arise. This first detailed study, combining the use of a time tested and widely used model, i.e., the WRF, partly overcomes this shortcoming since it is able to combine facilities for yielding interpolated data variables subsuming all available observational data. Additionally, it can handle convoluted topography and has a robust cloud microphysics scheme that can include super-cooled liquid drops along with cloud ice.
This present work is not the only study concerning the use of the WRF model for aviation-related applications, but definitely a first attempt at understanding the role of cumulonimbus induced icing on the aircraft surfaces using the model. A recent paper by Stocker et al. (2019) discusses the importance of the use of a Numerical Weather Prediction (NWP) model for forecasting the impact of adverse meteorological conditions on the movement of aircrafts around regions of complex orography. However, generic weather forecasting models are not usually configured to resolve micro-scale events. Two important papers by Wong et al. (2013) and Chan and Hon (2016) discuss the development and implementation of WRF-based high-resolution NWP model (also called the Aviation Model) for aviation meteorology at the Hong Kong International Airport where airflow-mediated disturbances occur partly due to a complex orography showing how sea breeze related wind shear and headwind change around the Airport.
Nepal has 32 operational airports, many of which lie in narrow valleys or on hilltops with heights ranging from 2440 to 3050 m Above Mean Sea Level (AMSL) (Yadav 2017). This arrangement forces aircrafts to fly through uneven terrain and narrow gorges like those found over regions near the Jomsom airport (see Fig. 1). The Jomsom airport in Nepal has witnessed several crashes attributed to the diurnal variation of the wind flow over the Kali Gandaki river valley (Regmi 2014b). Intense up-valley winds rush through a gorge and small aircrafts can hardly make their way to the airport. The 9 N-ABB crash was examined by Regmi et al. (2017), who linked it to severe turbulence induced by the undulating Himalayan terrain as well as icing on the aircraft. A. Chandra et al. Pure Appl. Geophys. However, the aforementioned author did not conduct an explicit analysis of flight performance when the aircraft was subjected to turbulence and icing but rather focused on simulating weather conditions prevalent on that day to deduce probable reasons for the aircraft crash. These are just some examples of the threats that aircrafts face from erratic weather systems. In total, from 1961 to 2016, there have been 273 fatalities from 37 incidents of turboprop aircraft crashes (Yadav 2017). Two incidents involving the Beechcraft 1900D aircraft (2004Beechcraft 1900D aircraft ( and 2011 operated by the Buddha Air airlines have led to 20 fatalities. The Beechcraft 1900D is still routinely used for scenic mountain flights originating from Kathmandu and is also the chosen aircraft for our own case study presented in this paper. Figure 2 shows the selected 19 seater Beechcraft 1900D which is routinely used by its operator, Buddha Airlines, for scenic flights that loop around Mt. Everest (Fig. 3), with the first flight scheduled at 6:30 am for a one hour flight. The aircraft is capable of carrying 19 passengers and the Buddha airlines alone, conducts 5 flights a day (Buddha Airlines 2020).
Deep cumulonimbus clouds occur frequently over the Himalayan region. The formation of these convective clouds has been studied in detail by Medina et al. (2010). It has been observed that from June to September, the number of thunderstorm days increases almost linearly with increasing latitude (Manohar et al. 1999). As thunderstorms are an excellent indicator of well-developed cumulonimbus clouds over a region, one can expect a lot of convective storms during this time period. Another study by Bhat and Kumar (2015) has shown that the foothills of the Himalayas and the northern parts of India see a high frequency of Cumulonimbus towers compared to the rest of the Indian subcontinent.
Cumulonimbus clouds are characterized by the presence of intense updrafts and downdrafts laden with super-cooled water drops, hail, and, graupel, and can extend up to the tropopause where a high wind shear gives it its distinctive anvil shape. Regions with mature cumulonimbus clouds experience intense precipitation and strong gusts caused by diverging downdrafts (Cotton et al. 2011). Aircrafts generally try to avoid cumulonimbus clouds as these clouds not only subject the aircraft to severe turbulence and icing but also obscure the pilot's range of visibility. The danger posed by Cumulonimbus clouds to aircrafts is duly recognized by the aviation industry and pilots are advised to avoid such clouds by finding alternative paths. A Delta Airlines flight (flight 191) crash (2nd August 1985) was attributed to the pilot's decision of continuing approach into a Cumulonimbus cloud rather than circumnavigating it (Ghazi and Juhany 1996).
Icing effects on performance are not easily detected in-flight and have often caused fatalities Tourists disembarking from a Buddha Air Beechcraft 1900D aircraft after an ''Everest Experience'' flight (Wang 2009) Vol. 177, (2020 A Case Study on Assessing Cumulonimbus Induced Flight Vulnerabilities 5043 (Cole and Sand 1991). Icing generally leads to a decrease in lift and an increase in drag which could then lead to a stall (Mingione and Barocco 2010). Aircrafts are equipped with de-icing and anti-icing systems on critical surfaces i.e. the leading edge of wings and windshields. However, unprotected sites like the underside of the wing can also increase drag significantly. Sand et al. (1984) observed that even with trace amounts of icing (with deicing systems enabled), the increase in surface roughness of the wing's underside was enough to increase the Coefficient of Drag (C D ) by a factor of 2, thus increasing fuel consumption. Cumulonimbus clouds have high Liquid Water Content (LWC) and aircrafts are susceptible to icing above the 0°C isotherm line i.e. the upper half of the cloud (Vogel 1988;Vukits 2002). Unsurprisingly, there has been no experimental aircraft-performance studies undertaken over this region wherein an instrumented aircraft con-currently measured ice amounts deposited, along with turbulence variability-none whatsoever over Nepalese terrain. However, a study by Todd (1964) involved the passage of an aircraft through a developing cumulonimbus cloud in light wind conditions. The study focused on the developing cloud structure rather than the effects of the cloud on aircraft performance. There have been some numerical studies- Wu and Cao (2016) studied the aerodynamic flow characteristics of a NACA0012 airfoil under heavy rain and icing and concluded that icing induces heavier aerodynamic penalties. The problem of aircraft icing is very much an ongoing issue and is still recognized as a major threat. This is evidenced from published scientific research papers i.e. Cao et al. (2018), who state that icing is still a major external cause of accidents-a 39% mortality rate in 42 icing incidents between 1986 and 1996. The authors indicate the most dangerous features and discuss causes of aircraft icing in tandem with environmental variability. Another major caveat is linked to ice crystal radar invisibility and was indeed the main factor precipitating the Air France Crash (Hylton 2011). Our paper is a case study showing how the problem of icing makes small and low flying aircrafts over the Nepalese Himalayas vulnerable. Although this is a unique case study, currently such flights routinely occur over this very popular route-about 5 per day just by Buddha airlines alone-and during the peak season when tourists still continue to enjoy scenic joy rides, cumulonimbus anvils also form very frequently (Manohar et al. 1999). In this sense, a case study based on a typical flight route, operated with a typical aircraft type and apprehending cumulonimbus anvils, is a good starting point. As the case study unfolds in the subsequent sections, it is shown how a state of the art cloud microphysics module within the framework of the WRF model helps in the quantification of the amounts of accreted ice (along with the mechanistic details of its formation from supercooled cloud droplets), associated deposition patterns, and the associated risks therein. This first attempt is by no means complete-there are avenues for improvement i.e. (1) by using a more sophisticated aerosol scheme within WRF-Chem and thereafter feeding in parameterized relations into the main module is a better way to count supercooled droplet numbers, (ii) a finer meshed grid over a smaller domain better resolves the up and downdraughts. These are some pointers on the way forward for stimulating further research in this multidisciplinary area albeit an overarching emphasis on how cloud microphysical impacts precipitate aviation hazards.
Presently, both ground-based as well as on-board radars systems are used for navigating an aircraft through thunderstorms-however, these afford only limited capabilities. The strength of a radar echo, or reflectivity (measured in dBZ) is directly proportional to the concentration and the radii of droplets, and can be correlated with the intensity of precipitation. For example, when the reflectivity exceeds 30 dBZ, one observes moderate rainfall. On-board radar systems do not yield a 3D profile of the associated hydrometeors. Another limitation of the system is related to limited signal attenuation when a cell absorbs or reflects a very broad spectrum of radio signals. This 2D representation can obscure a thunderstorm cell located immediately behind another cell (FAA (Federal Aviation Administration) 2013). Additionally, Ground-based radars return near-instantaneous observations of turbulence across large spatial extents. In a recent paper by Feist et al. (2019), it is shown how a statistical assessment of (the TKE dissipation rate) in convective as well as shower clouds show an increasing turbulence profile with values ranging from 0.01 m 2 s -3 at 1 km to 0.1 m 2 s -3 at 10 km for convective clouds, and a near constant vertical profile of with values ranging between 0.02 and 0.03 m 2 s -3 for shower clouds. We show later, in Sect. 3.3, modelled values of for the thunderstorm event pertaining to our own case study where similar trends are noticeable.
With a surge of affordable high-performance computing along with a tremendous progress in remote sensing operations, it is now possible to not only predict the onset of cumulonimbus activity, but also, to simulate the dynamical and microphysical environment within it (Cracknell and Varotsos 2007). Advances in gas chemistry, especially in the upper tropospheric regime, will further increase operation accuracy. For example, by accounting for diffused gases in ice will one can expect more accurate model estimates of ice morphology and may even lead to better de-icing mechanisms (Varotsos and Zellner 2010). Pre-run simulations, as well as real time simulations, will give airlines and pilots some control over uncertainties leading to efficient resource usage and safer flights.

A cumulonimbus Based Case Study Over Nepalese Terrain
Weather stations in Kathmandu recorded heavy rainfall on the 10th of July 2018 (2 am-4 am) with the presence of cumulonimbus activity-this was modeled with a WRF simulation which we now discuss in detail. Our case study involves a flight analysis of a Beechcraft 1900D aircraft (shown earlier in Fig. 2) passing through this particular cumulonimbus cloud. Figure 3 shows a schematic representation of the main study attributes. One will notice from Fig. 3, that an interlinked study to account for the in-flight icing is also in place. The latter involves a commercial CFD code, Star-CCM?, to accurately model the ice accretion process over the wing (CD-adapco 2018). In-flight icing develops from the cloud liquid water content (LWC). Accounting for microphysical cloud attributes, calls for a scheme well suited to the study objective, i.e. accurately accounting for accreted ice on the Beechcraft's surface. The Thompson and Eidhammer (2014) aerosol-aware scheme was used for a precise cloud microphysical representation of all the processes involved. This scheme includes an explicit module for the treatment of cloud microphysical processes detailed in the treatise by Pruppacher and Klett (2010). The scheme includes an explicit module for the treatment of cloud droplet activation and ice nucleation, crucial requirements for aircraft icing assessments. Cloud-ice nucleation parameterizations for heterogeneous ice nucleation are described by Phillips et al. (2008), while immersion freezing (above water saturation) follows the formulation by DeMott et al. (2010). The scheme also accounts for the freezing of deliquescent aerosols through parameterizations suggested by (Koop 2000) (Thompson and Eidhammer 2014). Figure 3 makes a clear distinction between the warm and cold portions of a towered cloud, separated by a 0°C isotherm. Warm regions consist of droplets grown by condensation onto aerosols. This portion of the cloud has only liquid water and therefore poses a minimal icing threat to the aircraft. The leading edge of its wings, susceptible to icing, has been marked in orange. The region of the cloud with temperatures ranging between -15 and 0°C consists mainly of super-cooled liquid droplets and pose 'a significant' icing threat to the aircraft. The cloudy region situated further aloft between -15 and -40°C features mixed-phase clouds, and pose 'some' icing threat to the aircraft. Aeronautical Information Manuals point out that, warmer temperatures between 8°and 14°C can be critical for aircraft performance at lower altitudes, relevant to this case where joyrides typically fly over low altitudes (TC AIM 2019). When an aircraft accretes liquid droplets along with supercooled ones in these low altitude flights, especially soon after takeoff, the liquid mass subsequently freezes during the aircraft's ascension well into the flight. Further, since thawing times are longer than freezing times, these frozen hydrometeors last through the aircraft's last leg as it descends to eventually land. This is, of course, more critical over parts of the aircraft devoid of de-icing mechanisms. Past studies have confirmed that aircrafts traversing through this chosen Nepalese route routinely encounter icing phenomenon (Regmi et al. 2017). To summarize, the accumulated ice may change the airfoil cross section, thereby marring lift, and increasing drag and stall speeds (Bragg et al. 2000; TC AIM 2019).
Additionally, aircrafts can also be affected severely by in-cloud turbulence associated with monsoon clouds as they navigate through eddies and vertical currents. Turbulence extends from the cloudbase to the cloud-top in towering Cumulus or Cumulonimbus clouds (see Fig. 3). In summary, an aircraft performance is typically linked to not only variations in hydrometeor size distribution and LWC, but also to the strength of the updrafts.
Sections 2 deals with the selection of the appropriate spatial and temporal domains required for the WRF analysis and these are garnered from an assessment of typical flight paths. Section 3 deals with the case study simulation and its validation with observed data, whilst Sect. 4 deals with surface icing on the aircraft involving the morphology of the deposition and its impact on aircraft performance-the morphology is crucial for an assessment of the roughness lengths impacting aero-dynamical drag. The 'Wider implication' section (Sect. 5) concludes the paper and talks about further work that can be done to aid pilots or the ATC to guide planes stuck in Cumulonimbus clouds. This section also proposes an outline of a step by step alerting protocol for pilots.

Study Region
The study focuses on a section of the Nepalese Himalayan region spanning the full extent of most low lying flight routes (27.71°N, 85.32°E to 27. 98°N , 86.92°E). varying from around 1400 m over Kathmandu to 8848 m along Mt. Everest. This large topographic variability in aerodynamic roughness lengths causes turbulence within the ensuing wind characterized by a succession of up and downdrafts through which the aircraft has to negotiate. Apart from shaping the weather system over the region, the high terrain also makes it hard for pilots to escape thunderstorms. This region is frequently battered by intense thunderstorm activity featuring deep convective clouds (Bookhagen and Burbank 2010;Palazzi et al. 2013). The test case chosen corresponded to a thunderstorm event that started off during the late evening of 9th July, during a monsoon depression with rainfall activity continuing into the next day, the 10th of July 2018. The region witnessed heavy rainfall as evidenced by the INSAT-3DR satellite image (MOSDAC 2018) and is shown in Fig. 5. This was the only available timeframe from INSAT-3DR-however, the cumulonimbus activity continued until the early hours of the 10th of July. The subsequent WFR simulations covered the entire duration. Figure 5 shows the GOES precipitation index (GPI) estimated from cloud top temperature measurements. GPI is an algorithm that estimates rainfall using only IR (Infrared) measurements from satellites. IR observations provide information about the location and temperature of the cloud tops. These observations are then correlated with radar and rain gauge measurements enabling an estimation of rainfall amounts using IR data (Arkin et al. 1994). This seems a reasonable procedure applicable to deep convective clouds over Tropical regions (Ba and Gruber 2001). The elongated black loop within the precipitation contour represents the flight path (Kathmandu-Mt. Everest-Kathmandu). One can see that the flight path crossed 2 bands of GPI contouring-it left off from a region where the precipitation index was 4.5 and within minutes approached an index of 7.5. Additionally, the aircraft flew across these bands at an altitude, placed between 5 and 6 km (see Fig. 4). At these heights, some of the hydrometeors are likely to be frozen (see Fig. 9). Upon, examining Fig. 9, one observes that the aircraft contacted air masses that were at first barely frozen (0°C at 5 km altitude) and ascended to very cold regimes, with temperatures of the order of -10°C at its highest cruising altitude of * 7.6 km. The lapse rate is thus * 0.4°C/100 m which is within the range of lapse rates generally found over the Himalayan terrain (Romshoo et al. 2018). This temperature regime is broad enough to initiate heterogeneous nucleation at the lower altitudes, mediated by copious amounts of aerosol pollution blown over from continental India along with the South-Western Monsoonal flow. Thereafter, at the higher cruising altitude, the temperatures were cold enough to initiate formation of some 'mixed phase' ice. Since the flight duration along the flight trajectory was of the order of 60 min, the flight afforded long enough contact times for the accretion of frozen hydrometeors in the form of ice as well as super-cooled droplets on the underside of the airfoil which did not have icing mechanisms in place.
It is now established that the AF447 crash killing 228 had been triggered by ice formation on the speed sensors. Icing from super-cooled water droplets can occur both during climbing as well as during descent. It is appropriate therefore to qualitatively discuss some mechanistic details of the freezing of supercooled water droplets (Hylton 2011).
The Nepalese Beechcraft Airplane came into contact with a low lying cloud (the maximum cruising altitude was 7-8 km where the temperatures varied between -2.3 and -14°C (Fig. 8). It is well known that super-cooled water droplets occur over sub-zero temperatures going down to even -20°C.
New research by Kintea et al. (2014) suggests that when such droplets in low lying clouds make contact with aircraft parts exposed to oncoming air currents then a portion freezes immediately. This exposure to air currents is very evident in the present case as the Beechcraft steered itself around the periphery of the Cumulonimbus event. It is quite likely that any nonfrozen amounts trickled across the wings causing a water film that eventually froze on the sub-zero temperatures of the underlying structure.
In contrast to this Beechcraft flight in Nepal, the AF flight flew at a much higher altitude. Kintea et al. (2014) suggests an interesting mechanism when planes fly over storm cells at cruising altitude. It is quite possible that a mixture of water and ice adheres to the aircraft surfaces, much like wet sand and further reduce the temperature of the metal surface until the freezing point is reached. Finally the water film freezes on the surface on which the cold ice particles are stuck. This seems a plausible mechanism for sensors to freeze as in the AF flight.
Having described the flight path, its incursions through a region of deep convection, we now move onto reconstructing the dynamical and microphysical attributes of the case study in question through a carefully set up WRF run. 2.2. WRF Model Set-Up for the 9-10th July Cumulonimbus Case Study In this case study, we have used the Advanced Research Version 3.6.1 (ARW) dynamical solver of the WRF code (Skamarock and Klemp 2008). It is in essence a mesoscale atmospheric model used the world over for numerical weather prediction and operational forecasting applications. We had to use 2 nested model regions to fully capture the event with horizontal grid spacing of 1350 9 1269 km 2 (outer parent grid) and 396 9 369 km 2 (inner daughter grid), as shown in Fig. 6.
As stated earlier, because of the presence of a succession of ridges and clefts over the Nepalese Himalayas, it is important to be able to fully resolve the associated up and downdraughts which have a direct bearing on the microphysics. For these reasons, it deemed appropriate to use 2 domains to resolve all of the large mesoscale features over this complex study area (Jiménez and Dudhia 2013). The simulation ran from 9th July 2018 1200 h to 10th July 2018 1200 h UTC spanning 2115-2215 h over the time domain. In order to precisely simulate the weather system over Kathmandu, 9 h of spin-up time was essential. The WRF has a vertical coordinate system with g levels marked from the surface up to 20 km. We have used the Final analysis (FNL) data for configuring the model's initial and lateral boundary conditions. This dataset is available at every 6 h over a horizontal grid resolution of 1°9 1°, covering 26 pressure levels from 1000 to 10 hPa (National Centers for Environmental Prediction/National Weather Service/NOAA/U.S. Department of Commerce 2015).
Although the WRF allows a choice of several microphysical schemes, we employed the 'Thompson   Grell and Freitas (2014) was used-this is a scheme facilitating a seamless transition from subgrid scales to larger scales.
Additionally, a recent paper by Thompson et al. (2017) also corroborates some of these options but specifically recommend the 'Thompson Aerosol Aware' scheme. We have therefore adopted this double moment scheme for the present analysis.
At this stage, it also relevant to review other small aircraft studies linked to icing instabilities from elsewhere, since no such aircraft based studies have been reported for this Nepalese terrain. Cooper et al. (1984) and Sand et al. (1984) steered research flights through icing conditions (over Montana and Northern California) using a similar aircraft (Beechcraft super king Air) and found that both the LWC and the mean volume diameter (MVD impact aircraft performance. They found that the rate of climb decreased by approximately 0.45 m/s per g/cm 2 of ice accumulation with a concomitant doubling of the Coefficient of Drag (C D ) with little to no change in the Coefficient of Lift (C L ). This led them to conclude even a thin layer of accreted ice onto parts of the aircraft (where de-icing systems were not in place, such as at the underside of wings, led to a significant increase in the surface roughness. Sand et al. 1984 also confirmed that acute icing was observed when the MVD values attained a threshold of 15 lm. Studies have also categorized the extent of icing on aircrafts-Jeck (2001) classed it as trace, light, moderate and severe for ice accretion rates 0-1, 1-6, 6-12 and [ 12 g cm -3 h -1 respectively. As we shall show later, the ice accretion rates for the present study fell under the moderate category.

Microphysical Characterization of the Cumulonimbus EVENT on 10th July 2018
The main thrust of this study is to assess the vulnerability of a Beechcraft negotiating a Cumulonimbus event. This involves an accurate quantification of the ice amounts deposited on the aircraft. In order to predict the right ice amounts, one must ensure that the modeled cloud morphology is captured accurately. This calls for a prior validational analysis. To this end, it is appropriate that modeled precipitable water vapour amounts are compared with observations. We show results from an atmospheric reanalysis using the Goddard Earth Observing System Model, Version 5 (GOES-5) (also called MERRA-2) (Global Modeling and Assimilation Office (GMAO) 2008). MERRA-2 is an hourly, vertically integrated diagnostics, available at a spatial resolution of 0.5°9 0.625°. The precipitable water vapor is a precursor for all forms of hydrometeorsalthough, as we show later in Sect. 3.2, in this instance, it primarily drove warm rain microphysics. This is evident from the flight path which lies embedded majorly in the warm region as can be seen in Fig. 8. For these reasons, it was appropriate to compare total precipitable water content from both the WRF, and, the MERRA-2 observations (Fig. 7). Figure 7 indicates a remarkable degree of agreement between observations and model resultsprecipitable water vapor observed at the starting point of the flight route i.e. Kathmandu (see Fig. 7a) ranges from 51 to 59 kg m -2 . A closer look at Fig. 7b, around the same location (27.7°N, 85.3°E), also shows similar water vapor amounts varying between 45 and 53 kg m -2 . Moreover, the total precipitable water vapor obtained from the WRF and MERRA-2 observations are also comparable along the entire flight path. This comparison attests to the robustness of the 'Thompson Aerosol-Aware' microphysical scheme. One observes that over regions spanning the flight path, the Total Precipitable water (Fig. 7) ranges from 20 to 60 kg m -2 . Since the temperatures encountered here range between 15 and -14°C, one can also infer that the precipitable water amounts can be linked to both supercooled water droplets as well as droplets generated from warm rain microphysics. However, it is important also to explore regions of the flight path where the Beechcraft would be most susceptible to instability induced by liquid and frozen hydrometeors. This is discussed in the next section.

Hydrometeor Distribution Along Flight Trajectory
The atmospheric temperature profile drives cloud microphysics within moist air. Prior to a discourse on the modeled hydrometeor profile, it is crucial also to look at isotherms along the flight path. Figure 8 shows the temperature contours on the 10th of July 2018 spanning four time slots.
From Fig. 8, one observes that during the first 3 quarters of the flight (27.72°N to 27.88°N), traversing the 1400-6000 m gradient, the temperature varies from 20 to -2.3°C, a temperature band where warm rain microphysics dominated. This was to the tune of 34%. Only during the last quarter of the flight path, when the aircraft covered a flight path from 27.84°N and 27.94°N, one notices sub-zero temperatures (-2.3°C to -13.6°C), where icedriven microphysics operated. Figure 9 shows vertical cross-sections of hydrometeor distributions along a west to east track of this flight for the same time slots. During takeoff, the aircraft collected liquid water which froze as it ascended to colder regions at a cruising altitude of around 7.6 km.
Panels b and c show hydrometeor distributions at 0320 and 0340 h where one notices reduced ice and precipitable water amounts along the entire west-east track. The aircraft has a total flight time of 60 min (30 min each way), and as it nears Mt. Everest at around 0330 h (Nepal time), it penetrates into regions with freezing temperatures mediated by cold-rain  processes alone. When Super-cooled Liquid Droplets (SLD) in low lying clouds make contact with aircraft parts exposed to oncoming air currents then a portion freezes immediately. This exposure to air currents is very evident in the present case as the Beechcraft steered itself around the periphery of the Cumulonimbus event. It is quite likely that any non-frozen amounts trickled across the wings causing a water film which eventually froze on the subzero temperatures of the underlying structure . The other causes of concern relate to the poor radar detection of ice signatures and the problem of ice incursions in aircraft sensors.
As we move on to panel d (0400 h), we notice a further reduction in the hydrometeor mixing ratios as the aircraft approaches Mt Everest (27.97°N, 86.68°E ). One observes that the Beechcraft majorly traverses through regions with freezing temperatures (to the tune of -10°C), from 27.92°N to 27.98°N. This stretch is mainly dominated by high concentration of snow particles, which have the propensity to stick to the aircraft's underside and wings. Cao et al. (2018) present a detailed review on the types of aircraft icing and the associated hazards. The authors also mention the statistically most dangerous features and causes of aircraft icing. They also list out the . Also, note the positioning of the isotherms (obtained from the WRF model) shown in grey where a 0°C line can be seen at an altitude of 5 km. Warm-rain microphysics is operative below 5 km where one notices cloud liquid water (shown in red) and rain mixing ratios (shown in green). Cloud tops extend up to a maximum of 16 km (see a) and gradually lowers to a height of 14 km as the cloud dissipates (see b-d) Vol. 177, (2020) A Case Study on Assessing Cumulonimbus Induced Flight Vulnerabilities 5053 environmental conditions and causes of several aircraft crashes caused by ice accretion. Specifically, they discuss the formation of three main types of icing: rime, glaze and mixed ice. The formation of rime ice is mediated by the presence of super-cooled droplets located in regions with sub-zero temperatures. In the present Nepalese case study, we notice from Fig. 9 that the aircraft traverses through regions with below freezing temperatures during its final leg of cruise prior to landing. During these times, the Beechcraft encounters SCDs which freeze immediately upon contacting the surface of the wings. On the other hand, the formation of glaze ice occurs in warmer conditions (see isotherms below 6 km in Fig. 9) where liquid droplets flow along the wing's surface under the influence of aerodynamic forces. These overflowing droplets cluster upon freezing to form appendages on the airfoil structure (Cao et al. 2018). We indeed observe similar results in our study in Fig. 16, where one observes double-horn appendages on the Beechcraft's airfoil. Lastly, mixed ice is formed in clouds with temperatures between -20 and -10°C. Since we are concerned with a light and low-flying aircraft with a maximum cruising altitude of less than 8 km and with a temperature variation between -2.3 and -13.6°C, it is clear that over this temperature band the existence of Supercooled Liquid Droplets (SLD) would indeed prevail over mixed ice.

Cloud Water, Vertical Velocity and MVD Distribution Along Flight Trajectory
Having established that the Beechcraft primarily encountered, hydrometeors initiated from both warm and cold microphysical processes, it is worthwhile to explore other tangible indicators for ascertaining the extent of precipitation that the aircraft had intercepted. The most obvious parameter that pilots look for is the radar reflectivity (Kumar et al. 2013;LeMone and Zipser 1980;Sauvageot and Omar 1987). Radar reflectivity contours provide useful information on the microphysical properties of clouds, for instance, a radar echo (measured in dBZ) whose strength is directly proportional to the concentration and the radii of droplets, correlates directly with the intensity of precipitation. One can also distinguish between portions of the cloud with more liquid water than others with lower amounts, and for obvious reasons a flight path should operate through a chosen track comprising of the least amounts of cloud water. Matsuda and Onishi (2019) show an increase in the radar reflectivity factor over regions of high liquid water content (LWC), particularly over regions where a non-uniform distribution of droplets is observed. A large value of radar reflectivity is attributed to scattering by a higher concentration of droplets. These authors have reported a high total reflectance to the tune of 20 dBz in the upper reaches. This observation corroborates with what we have also shown for this Nepalese case study (see Fig. 10a-d). We observe similar results from Fig. 10a-d where values range between 20 and 40 dBz over a same vertical stretch. Additionally, the very recent Matsuda and Onishi paper indicates cloud liquid amounts in excess of 0.1 g m -3 , an observation again consistent with our own findings (see Fig. 11a-d). We also extracted the radar reflectivity values associated with this Cumulonimbus event for the entire duration of this flight (from takeoff to cruising altitude prior to its landing). Figure 10 shows the vertical transects of the modeled radar reflectivity (dBz) for four time slots during the Cumulonimbus activity. Ideally, one should additionally have dual Doppler/polarization measurements where one can specify crystal size and habits (Moisseev et al. 2015;Posselt et al. 2015;Sinclair et al. 2016).
Upon a closer look at Fig. 10 a, one observes moderate to heavy rainfall around 27.72°N reflected by reflectivity values ranging between 40 and 50 dBz. In fact, the region reported heavy rain on the 10th of July 2018 during the same time period. However, as we move from panels b-d (0320-0400 h), we see reduced reflectivity values signifying a state of dissipation. At 0400 h (Fig. 10d), one observes very light to moderate rainfall, corresponding to radar reflectivity values spanning between 20 and 40 dBz marked by wispy clouds left behind.
A past study on cloud reflectivity has also shown that a threshold value of 20 dBZ at an altitude of 12,000 m is a strong indicator of the presence of deep convective clouds (Bhat and Kumar 2015). From Fig. 10, we indeed observe a similar value of cloud 5054 A. Chandra et al. Pure Appl. Geophys. reflectivity at an altitude extending from 12,000 to 15,000 m. Since these comparisons are encouraging, we now proceed with further microphysical assessments for accessing flight vulnerability. We discuss first, the variability in LWC viz a viz the morphology of the updrafts and downdrafts. As stated earlier, this varying water-signature would leave a frozen imprint subsequently on the Beechcraft, sometimes on undesired parts, devoid of de-icing circuitry.
From panel a (0300 h), one observes four regions of high LWC, vertically extending up to 6000 m, and horizontally spanning between 27.72°N and 27. 88°N . As we move from panel b-d (from 0320 to 0400 h), a reduction in cloud water amounts is noticed corresponding to a dissipation of thunderstorm activity. We use analytical formulations to estimate this time. It is found that the cloud dissipates over a total duration of 60 min (Ghosh et al. 2018;Ghosh and Jonas 1998). The Aerosol Aware microphysical scheme yields cloud water autoconversion rates ranging between 10 -7 and 10 -4 g m -3 s -1 (Lee and Baik 2017). The modelled time scale for this rate (* 60 min) matches very favourably with analytical estimations of Ghosh et al. (2018).
It is observed that as the aircraft takes off from Kathmandu, it ascends through regions of intense updrafts. Updraught speeds directly correlate with the strength of the vertical mean winds-when the latter are under predicted, then once expects a concomitant under prediction in the former. Yang et al. (2013)  reported underestimation of WRF-modelled vertical winds between 600 and 1000 km AGL in both stable as well unstable atmospheric conditions. If one examines Fig. 12, one notices that the updraft velocity approaches 1 ms -1 soon after takeoff. This rising motion of air quickly initiates condensational growth and results in increased cloud water amounts. This is also evident from Fig. 11a, where one notices a zone of high LWC around the same region containing an updraft. The flight path covers four updraft transects encountering high LWC, and upon attaining cruising altitudes also encounters downdrafts, where the LWCs are much lower. This indicates that the succession of updrafts and downdrafts has a direct bearing on the extent of liquid water generation (shown in Fig. 11) modulating in turn the extent of supersaturation. However, even after some decrease in the liquid water amounts over a period of 40 min (see Fig. 11c) the third band of updraft still has significant liquid water, which could adhere to aircraft parts. As the aircraft ascends further upward to the cruising altitude of 7.6 km, the accreted liquid water would have got entirely frozen. It is in these higher reaches, that the telltale anvil shape, associated with thunderstorm activity, is expected to form. The simulations affirm this anvillike protrusion in the higher reaches quite clearly.
Since the super-saturation profile is perturbed by in-cloud turbulence, it is also interesting to look at the turbulence characteristics of the thunderstorm in question. Figure 13 is particularly revealing; it shows a modelled vertical profile of up to 100 hPa for this showering thundercloud as observed on the 9th July 2018 at 0320 h UTC. Figure 13 illustrates two distinct regimes-a 'showering' zone up to 300 hPa (10 km) and a deep convective regime spanning 300-100 hPa (10-16 km). The former, subsumes a vertical stretch of near-constant values of the order of 0.01 (m 2 s -3 attributed to a dissipating thundercloud. On the other hand, a well-marked deep convection zone is also present marked by increasing values, first starting with a minimum value of * 0.01 (m 2 s -3 at 300 hPa to a maximum value approaching 0.08 (m 2 s -3 ) at 100 hPa. This is also seen in Fig. 12b where strong updraughts at these altitudes indicate zones of intense convection. These trends are consistent with radar observations by Feist et al. (2019) strengthening the case study discussed here.
The Beechcraft although bigger than a glider aircraft can also be severely buffeted by up and down draughts caused by the Himalayan orography. In particular, shear effects between updraughts and downdraughts can cause Kelvin-Helmholtz instabilities that can severely damage an aircraft. Additionally, when there is a strong updraught below a super cell cumulonimbus, a small plane can even be sucked into the cloud. Importantly, this case study concerns a tourist circuit over one of the most impoverished regions of the world where a heavy demand exists for such scenic flights sought after by European and American Tourists (paid not in local currency, but in USD). Beechcraft used do not have sophisticated sensors on board. These small private aircrafts are generally not equipped with on-boardweather-radars. Cumulonimbus tops are invariably frozen-the vertical extent of the Cb in question varied from 2 (cloud base) to 16 (cloud top) km. In these lofty realms water droplets will invariably freeze and the associated phase change makes the cloud top very turbulent. Although, not relevant to this event, small aircrafts are often forced to ascension close to cloud tops when it gets covered with ice and the controls can freeze and remain immobile or even stuck. Following on from the preceding discussion, comparing the LWC viz a viz the morphology of the updrafts and downdrafts, we now discuss the variability of the droplet mean volume diameter across the flight path.
From Fig. 14, it is noticed that the aircraft makes contact with both small and big hydrometeors ranging between 14 and 42 lm during its ascent. These hydrometeors are majorly cloud droplets contained within warm regions of the atmosphere. Thereafter, as the aircraft ascends to its cruising altitude, it is exposed to droplets with much reduced diameters in the range of 4.75-19 lm.

Aircraft Icing
The scenic joy rides connecting Kathmandu and Mt. Everest typically hover along low altitudes allowing tourists to photographically capture a spectacular snowline. During the flight duration, the aircraft accretes liquid droplets along with some super cooled droplets as it flies through the prescribed route. These accreted droplets are left to spread out and mingle with other drops before they freeze. Further, thawing times are longer than freezing times, and these frozen hydrometeors last through the aircraft's last leg as it descends to eventually land. The structural accretion of ice destroys the smooth flow of air around the aircraft's airfoil, thereby decreasing its ability to induce a lift force requiring additional power to overcome the ensuing drag force. In this section, we present an icing simulation of an airfoil section of the 1900D Beechcraft wing in order to demonstrate how the microphysical conditions in this cumulonimbus case study could actually promote the growth of dangerous protuberances on the leading edge of the aircraft. The methodology shown in this case study should help other interested researchers to model the prevailing dynamics and the associated microphysics of a small aircraft negotiating Cumulonimbus events, and generally, to have a handle on the overall icing conditions, and specifically, to A. Chandra et al. Pure Appl. Geophys. predict the shape of the accreted structure on vulnerable aircraft parts. Patterns of ice deposition on vulnerable parts of an aircraft affect aircraft performance. In particular, any extraneous morphological additions on the airfoil shape can prove to be dangerous (Mikkelsen et al. 1985). It is therefore very important for us to explicitly explore this caveat in our study.
A computational fluid dynamics (CFD) program known as Star-CCM? was used for deriving the morphology of the deposited ice on the Beechcraft. For our analysis, a CAD model of the wing of Beech 1900D aircraft was imported into the software and then meshed (as seen in Fig. 14). In order to simulate the flow around the airfoil (NACA23018), and thereafter quantify ice-accretion over the wings, an implicit-unsteady simulation was carried out with the k-2 turbulence model with dispersed Multiphase interaction enabled.
The input conditions to the icing software were imported from the WRF simulations-these include the Temperature, Pressure, the LWC and the MVD of  the hydrometeors for each of the three flight-segments i.e. take-off, cruise, and landing (see Fig. 4 for details of the flight segments). Table 2 summarizes the chosen input configurations obtained from the WRF model. The 2D airfoil along with the mesh is shown in Fig. 14.
To visualize and predict the growth of ice accretion, a Morpher displacement plot was generated with 90 s of icing time for each of the three flight segments i.e. take off, cruise and landing. Morphing technique is employed in applications of fluid dynamics where it is essential to monitor mesh deformation caused by a transient flow of fluid or accretion of solid on to any prescribed geometry (Fig. 15). The Morpher functionality in Star-CCM allows one to include the effects of moving boundaries in a transient simulation also. The function changes the shape of the mesh (and thereby the geometry) to account for the accreted ice in the previous time steps. The updated geometry with new mesh structures is then used for subsequent iterations. The quiver plots in Fig. 16 illustrate the Morpher geometry comprising of the accreted ice shown as appendages on the Beechcraft's airfoil. The vectorial representation yields information on the depth of the accreted ice and the direction of its further growth. For instance, the red vectors in Fig. 16a (Take off) correspond to an ice-layer with a maximum thickness of 3.3 lm, with a South-westerly accretion. During cruise and landing (see Fig. 16b), a 42% thicker icelayer (4.7 lm) is noted on the Beechcraft's airfoil. It is clear that over cruising altitudes, super-cooled droplets on making contact with aircraft parts exposed to oncoming air-currents freeze immediately, thus increasing the ice thickness (evident from the displacement vectors). Figure 16 is particularly interesting; one notices a vectorial representation (Morpher plots) of ice thicknesses along the leading edge of the Beechcraft's airfoil during 'takeoff' and 'cruising'. The mean ice thickness values (for 90 s of icing simulation) during these aforementioned phases of flight trajectory are 2 lm and 3 lm, respectively. With this information, one can ascertain the time rate of accretion of ice on the aircraft's wings as a product of the mean ice thickness (lm), the density of ice (kg m -3 ) and the total icing time (s). This yields the

Figure 16
Morpher displacement during take-off, cruise and landing (ice growth) showing ''Double Horn'' appendages on the Beechcraft airfoil Vol. 177, (2020) A Case Study on Assessing Cumulonimbus Induced Flight Vulnerabilities 5061 flux of the time rate of accretion typically * 0.1 g m -2 s -1 on the Beechcraft's wings. Thawing times are large and the frozen ice on the aircraft would impact its stability due to the added aerodynamic roughness. However, with this level of sophistication, were we able to discern a double horn appendage on the airfoil as indicated in Fig. 16. This ungainly accoutrement is characteristic of glaze ice and is known to severely disrupt the airflow by affecting engine torque and the angle of attack (Cao et al. 2018;Mikkelsen et al. 1985).

Conclusions and Wider Implications
This first study, combining the use of a numerical weather prediction model (WRF) with an ice simulation procedure, (Star CCM?) yields newer insights into flight-vulnerability assessments over a highly convoluted Himalayan terrain. The unique orography of the Nepalese Himalayan terrain is conducive to the formation of cumulonimbus activity during the South-West (SW) monsoon season. Since tourism is a major industry in cash-starved Nepal, short-haul flights are used routinely to ferry tourists along lowaltitude flights using small Beechcraft airplanes, along well-marked scenic air-routes, enabling mountainscape photography. Unfortunately, the government does not issue travel advisories, unless, when it comes face-to-face with severe weather hazards. Such flights take off and land even when moderate cumulonimbus activity is forecast. This study shows for the first time, how the associated hazard arises is not only from severe in-cloud turbulence associated with cumulonimbus clouds, but also are linked to ice deposition on parts of the aircraft devoid of de-icing mechanisms. For these shorthaul flights, deposited ice is not able to thaw, and as we have shown, a double horn appendage on a Beechcraft's airfoil can cause severe instabilities.
We demonstrate these effects through a microphysically modeled cumulonimbus observed on the 10th of July 2018 during the SW monsoon season. Importantly, we show how the flight encountered liquid hydrometeors during its ascent and subsequently froze at cruising altitude gaining ice-mass.
Flight paths can be altered to completely avoid a Cumulonimbus cell, or, if it is deemed sufficiently safe, a path can even be charted through it. Lim and Zhong (2018) present a route re-planning method which suggests avoiding specific regions of a convective airspace-these zones include prohibited areas, restricted areas, and the danger areas. The authors use a cellular automation model to chart flying zones which do not intersect with the three aforementioned regions. The airspace is divided into a cell matrix where the state of each cell is dependent on its eight neighboring cells in the up, down, left, right and diagonal directions. The state of the central cell can be determined based on the spatial coverage of the radar reflectivity contours, such as the ones shown in Fig. 10 of this Nepalese case study. Regions with moderate to high reflectivity can be treated as restricted zones, whereas zones of intense reflectivity can be marked as prohibited areas. An optimal speed trajectory can thus pass through regions surrounding low to moderate radar reflectance values. Another crucial factor influencing the aircraft speed is the strength of winds which can be quite erratic inside a cumulonimbus cloud-if the airspeed is low and it suddenly drops due to a wind reversal, the aircraft can stall leading to a loss in altitude. If the airspeed is too high and the aircraft is hit by severe turbulence, the induced high structural stresses could damage the aircraft. The present study shows how a Numerical Weather Prediction Model can be used to estimate strength of the winds, particularly the convective updraughts/downdraughts inside a cumulonimbus cloud. Ground-based radar returns along with model results can yield a comprehensive picture showing zones of intense convection. This additional information along with the former route planning method can help in chalking an optimal speed profile, allowing pilots to conveniently decide their chosen flight trajectories.
A natural progression from this study should involve a remedial strategy. We provide this through a sequence of steps. It is proposed that a product be developed to help pilots and ATC personnel to tackle problems associated with flights navigating thunderstorms. The schematic diagram shown in Fig. 17 illustrates the working of such a product. The flight paths can be optimized, and a speed profile could be 5062 A. Chandra et al. Pure Appl. Geophys. created for a safe and fuel-efficient journey through thunderstorms.
For an aircraft scheduled to fly through the Cumulonimbus event on the 9th-10th July, high risks could be averted with a mere flight postponementfrom Fig. 14 a good starting time would be an hour later when the cloud would have dissipated much of its kinetic energy. Had that not been possible, an alternate route avoiding super-cooled liquid droplets could have been charted. This is quite manageable from the data presented in Figs. 9, 11, 12, 13 and 14, along with the associated CFD simulations using Star-CCM ? enabled to quantify any loss of performance.
This is a first study and to our knowledge, also a first attempt. Although we have shown how small aircrafts plying through a much loved Himalayan terrain, are very vulnerable to the problem of aircraft icing, we can also explore how such studies can be improved. At the outset, new research should be undertaken to include additional de-icing mechanisms along vulnerable aircraft parts in addition to the ones already in place. This need has been discussed in Sect. 2.1 where it was said that as aircrafts encounter (and also circumnavigate) Cumulonimbus mediated thunderstorms, particularly, over relatively low cruising altitudes (required for a panoramic view of the Himalayan snowline), liquid water as well as super-cooled droplets trickle to wings as well as to other sensors and eventually freeze. Whilst placing de-icing devices on wings of such small crafts may still be possible, equipping them within sensors still remains undone. We now discuss the second theme which can further improve related studies. This concerns pre-flight updates. Whilst, the WRF model shall still continue to be the chosen tool, it can be much better customized. We anticipate a surge in computational power within the sub-continent. First, parameterized results from WRF-Chem yielding accurate aerosol counted precipitation amounts must be used. We suggest the use of new parameterisations linking cloud auto-conversion rates to aerosol spectral properties (mainly mode radius and dispersion) to quantify cloud water amounts. Currently, the standard version of the WRF under-predicts precipitation amounts (Conrick and Mass 2019). For an exact quantification of the extent of icing, one must first accurately quantify the expected amounts of the main hydrometeors. Finally, a better future study should use a version of the WRF operating over a more localized domain (see Li et al. 2019). Although in this paper a staggered vertical grid was used, the Flowchart describing the process of aircraft path planning in real-time or in advance for any weather condition Vol. 177, (2020) A Case Study on Assessing Cumulonimbus Induced Flight Vulnerabilities 5063 horizontal resolutions were rather coarse (see Table 1). Using a newer version and Local version of the WRF shall have an additional benefit-Orography modulated up and downdraughts can be much better resolved using the embedded LEM option in this version.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.