*Article* **Analysis of Mountain Wave E**ff**ects on a Hard Landing Incident in Pico Aerodrome Using the AROME Model and Airborne Observations**

#### **Jin Maruhashi 1, Pedro Serrão <sup>2</sup> and Margarida Belo-Pereira 3,4,\***


Received: 22 May 2019; Accepted: 24 June 2019; Published: 26 June 2019

**Abstract:** A hard landing incident in Pico Aerodrome (LPPI) involving an Airbus A320-200 aircraft is investigated using airborne observations and forecasts of the AROME (Applications of Research to Operations at Mesoscale) model. A second flight is also analyzed. The severity of the wind shear during both flights is quantified using the intensity factor "I" that is based on aerial data and recommended by ICAO (International Civil Aviation Organization). During Flight 1, 36% of the landing phase (below 2100 ft) occurred under "severe" wind shear conditions and 16% occurred under "strong" conditions. Upstream characteristics included southwest winds, stable stratification and a Froude number close to 1. According to the AROME model, these circumstances triggered the development of vertically propagating mountain waves, with maximum vertical velocities above 400 ft/min and exceeding 200 ft/min in the flight path. These conditions, together with the severe wind shear, may have caused the incident. During the second flight, a wake with lee vortices and reversed flow developed in the region of the flight path, which is consistent with a low upstream Froude number and/or with the flow regime diagram of previous studies. During the approach phase of this flight, "severe" wind shear conditions were absent, with "strong" ones occurring 4% of the time. It predominantly displayed "light" conditions during 68% of this phase. As a result of the comparison between "I" and the AROME turbulence indicators, preliminary thresholds are proposed for these indexes. Lastly, this study provides an objective verification of AROME wind forecasts, showing a good agreement with airborne observations for wind speeds above 10 kt, but a poor skill for weaker winds.

**Keywords:** mountain waves; wake; AROME model; turbulence indicators; Azores Island; Froude number; model objective verification; aircraft observations

#### **1. Introduction**

In a stably stratified atmosphere, the airflow towards a topographic obstacle may trigger gravity (or buoyancy) waves downstream of the obstacle, known as mountain waves [1]. The flow response to the obstacle depends on several factors like the mountain shape and the Froude number [2]. At a low Froude number, several phenomena have been observed: the flow may be blocked upwind [3] or it may also pass around the obstacle rather than over it (flow splitting phenomenon) [2,4]. In association with the flow splitting phenomenon, a low-level wake may develop on the leeside of the obstacle [4–7]. The wake can be expressed as a quasi-steady pair of counter-rotating vortices circulating about vertical

axes or can present an unstable pattern in which vortices of alternating sign are periodically shed downstream to form a vortex street [4,5]. Also at a low Froude number, wave breaking, indicating the breakdown of laminar flow, may occur above the obstacle [2,4]. On the other hand, at high Froude numbers, airflow will overcome the topographic obstacle and vertically propagating mountain waves may exist depending on the relation between stability, wind and mountain width [8].

The study of such orographic phenomena is of paramount importance in the realm of aeronautics as they have been known to cause a variety of aviation accidents worldwide [9,10] and are especially hazardous to light aircraft during their takeoff and landing phases. Therefore, several authors analyzed aircraft accidents and incidents in various parts of the world, with special attention given to airports or aerodromes directly exposed to orographic flow interference. Aircraft observations have often been compared to NWP (numerical weather prediction) models in an attempt to not only establish their reliability as forecasting tools but to also better comprehend the main causes behind these incidents. Keller et al. [9], for instance, investigated a Continental Airlines aircraft that caught fire after veering violently off of the runway in Denver International Airport (DIA) during takeoff. High-resolution numerical simulations enabled them to conclude that lee waves were the probable cause behind the strong crosswinds that caused the accident. Another study published in 2013 by Parker and Lane [10] revealed that turbulence engendered by orographic flow resulted in a fatal lightweight aircraft accident in Australia.

In Portugal, the behavior of Madeira's wake has been the subject of several studies [5,11,12]. However, mountain waves associated to Pico Mountain (the highest of Portugal), to the authors' best knowledge, has only been addressed previously by Barata et al. [13] in a characterization of the surrounding wind regime and validation of the existence of leeside vortices through a scaled physical model of the island placed within a wind tunnel.

Pico Island is one of nine constituents of the Azorean archipelago, situated in the Atlantic Ocean 1500 km from Lisbon (Figure 1). Pico Aerodrome (LPPI) lies on the north coast of the island (Figure S1), north of Pico Mountain with a maximum height of 2351 m (7710 ft). In the presence of southwesterly winds in the southern region of the Island, which are frequent [13], the likelihood that takeoff and landing phases will be affected by a plethora of orographic phenomena ranging from severe mountain waves, leeside vortices and overall LLWS (low-level wind shear) is increased [13,14].

**Figure 1.** Macroscopic view of Pico Island's location relative to Mainland Portugal.

Forecasts of mountain wave induced turbulence rely on algorithms based on outputs of numerical weather prediction (NWP) models. Many of these algorithms rely on wind, temperature and turbulent kinetic energy [15]. Therefore, the accuracy of these fields is a key issue. However, an accurate simulation of such atmospheric variables is particularly challenging in the atmospheric boundary layer (ABL) due to its unsteady and turbulent nature [16]. Several factors may influence the model accuracy, namely, unresolved topographic features [17], turbulence closure schemes [18,19] as well as the representation of turbulence mixing [20,21]. The dynamical core also plays a key role. Ricard et al. [22] for instance, in analyzing the kinetic energy spectra of AROME and Meso-NH models (which have the

same physical parameterizations), concluded that a model using a Semi-Implicit Semi-Lagrangian (SISL) scheme has a coarser effective resolution than a model using an Eulerian explicit scheme (less efficient).

The current article focuses on the orographic flow associated to Pico Mountain in two different synoptic backgrounds, using forecasts from the operational AROME (Applications of Research to Operations at Mesoscale) model and airborne measurements collected by SATA (SATA Internacional–Azores Airlines, S.A.) Airbus A320-200 aircraft. Two flights were considered, where the first one represents the hard landing incident (henceforth referred to as Flight 1). Upon touchdown, the aircraft bounced and registered a vertical acceleration above 2 G. No fatalities or serious injuries were reported. The second one, labeled Flight 2, will serve as a control parameter that is representative of a flight in which no incidents were reported.

This study has three goals. The first goal is to characterize the nature of the orographic phenomenon involved during the two flights. The second is to provide an assessment of the agreement between AROME wind forecasts and airborne observations. The third objective is to establish a correspondence between AROME turbulence indicators and the wind shear severity registered by the aircraft. For this purpose, five turbulence indicators (Brown, Ellrod TI1, Ellrod TI2, CAT1 and EDR or Eddy Dissipation Rate) are computed from AROME forecasts. The wind shear perceived by the aircraft is quantified using the wind shear intensity factor "I", as originally suggested by Woodfield and Woods [23]. This assessment is also of great importance because the Portuguese Institute for Sea and Atmosphere (IPMA) routinely provides AROME forecasts of the wind and EDR for Pico and Faial Islands to SATA on an experimental basis. This test phase started in November 2017.

#### **2. Data and Methodology**

This section characterizes the wind from aerial data gathered by SATA A320-200 aircraft during both flights. The method of calculation of the wind shear intensity factor "I", which is a function of these airborne measurements, is established. A brief description of the AROME model is also provided. The Froude number and turbulence indicators, which are based on forecasts, are then defined. Lastly, definitions for the RMSE (root mean square error) and MAE (mean absolute error) concerning wind magnitude and direction are included.

#### *2.1. SATA Airborne Measurements*

For each flight, wind intensity and direction were recorded approximately every 4 s along with the corresponding altitude and geographical coordinates (latitude and longitude). Only the approach and landing phases (approximately final 10 min of flight or altitudes below 5000 ft) were perused since these are the most pertinent to understanding the incident. The average trajectory for both flights along with Pico's orography in the AROME model are mapped in Figure 2.

**Figure 2.** Approximate flight path (blue dashed line) and AROME orography. Shading shows elevation in meters. Pink marker shows the aerodrome location. "AB" and "CD" are endpoints of cross-sections analyzed in subsequent sections.

#### *2.2. Wind Shear Intensity Factor "*I*"*

The severity of wind shear intensity, as perceived by the aircraft, will be measured according to the wind shear intensity factor "I", detailed in ICAO's (International Civil Aviation Organization) Manual on Low-level Wind Shear [24]. Following Guan and Yong's study [25], wind shear is classified as "light", "moderate", "strong" or "severe" based on the variation of the airspeed with time *dVTAS dt* and on the airspeed proportion <sup>Δ</sup>*VTAS VTAS* . The intensity factor "I" itself may be expressed as follows:

$$I = \frac{dV\_{TAS}}{dt} \cdot \left(\frac{\Delta V\_{TAS}}{V\_{TAS}}\right)^2\tag{1}$$

The computation of the true airspeed, *VTAS*, is based on the following relation:

$$GS = V\_{TAS} + V\_{HT\_2} \tag{2}$$

where *GS* represents the ground speed of the aircraft and was computed by employing the three-dimensional distance formula between any two points described by their altitude and geographic coordinates (latitude and longitude). Note that Δ*VTAS* is simply the difference between true airspeeds at the *n*th and (*n* + 1)th instants. The head or tailwind component (see Equation S2.1 and Figure S2, from Supplementary Materials) is then subtracted from the ground speed to obtain the true airspeed.

Once *dVTAS dt* and <sup>Δ</sup>*VTAS VTAS* are known, the wind shear intensity may be classified by consulting the proper boundaries specified in Figure S3 from the Supplementary Materials. In both flights, this was performed for altitudes approximately below 2000 ft (610 m), which according to Guan and Yong [25], represents the appropriate range for low-level wind shear considerations.

#### *2.3. Description of the AROME Model*

The AROME model that runs operationally at IPMA uses a horizontal grid spacing of 2.5 km. It uses a SISL scheme and employs a three-class ice scheme with ECMWF (European Center for Medium-Range Weather Forecasts) radiation parameterizations, as was described by Seity et al. [26]. The representation of turbulence in the ABL is based on a prognostic TKE (Turbulent Kinetic Energy) equation [27], combined with a diagnostic mixing length. AROME initial and boundary conditions emanate from the ARPEGE (Action de Recherche Petite Échelle Grande Échelle) model. At the time of Flight 1, the operational AROME model was integrated with 46 levels, while at the time of Flight 2, additional vertical levels were incorporated, leading to a total of 60 levels.

#### *2.4. Froude Number Definitions*

The Froude number is a dimensionless parameter that may be utilized to describe the interaction mechanism between kinetic and potential energies as air flows past orography, according to Sutherland [28]. The Froude number can be interpreted as the ratio of the inertial to gravity forces in the flow (Lynch and Cassano [29]). This ratio may also be interpreted physically as the ratio between the mean flow velocity and the shallow water gravity wave speed (Holton [30]). The Froude number has also been applied to infer the type of flow phenomenon that may occur in the lee of a mountain when paired with a flow regime diagram, as was done by Sheridan and Vosper [31].

Three definitions of the Froude number were considered in this study. The first is the classical form, the second is based on Sheppard's dividing-streamline concept [32] and the third is applicable when a temperature inversion is present in the vicinity of the upstream region of the mountain. The exact upstream location at which Froude numbers were computed in both flights is indicated by point C in Figure 2.

#### 2.4.1. Classical Froude Number, Fr

The classical Froude number, Fr, may be defined as:

$$\text{Fr} = \frac{\overline{\mathcal{U}}}{\mathcal{N}h\_m} \text{\textquotedblleft} \tag{3}$$

where *U* is upstream flow speed, *N* is the Brunt-Väisälä frequency and *hm* is the mountain height [6,7,33]. Vosper et al. [34] provide an additional definition, *FL*, that is based on a horizontal length scale *L*, such that *FL* = *U*/*NL*. In this study, *hm* is 1200 m, which is the maximum height of Pico Mountain according to AROME topography and *L* = 21 km. Some authors [2,4,35], alternatively, refer to the non-dimensional mountain height *M* as:

$$M = \frac{1}{\text{Fr}} = \frac{Nl\_m}{lI\_{nm}},\tag{4}$$

where *N* would now represent an average value for the Brunt-Väisälä frequency below *hm* and *Unm* a cross-mountain wind speed average also under *hm*, as is done similarly by Jiang et al. [35].

Several authors summarize their results from laboratory, numerical and theoretical studies in a schematic regime diagram showing that for Fr > 1 mostly small-amplitude gravity waves are expected, while for Fr < 1, flow splitting, lee vortices and/or wave breaking may occur depending also on the horizontal aspect ratio of the obstacle. This is seen in Figure 2 from [4] and in Figure 13 from [2], the latter of which is shown here as Figure A1a.

For the parameter *FL*, whenever *FL* ≤ 1 vertically propagating waves are expected, but are evanescent if *FL* > 1 [34].

#### 2.4.2. Froude Number (Frh) in Terms of the Dividing-Streamline Height *hc*

Under stable conditions and with varying wind speed and stratification, according to Sheppard [32] and Heinze et al. [6], it is possible to determine the dividing-streamline height *hc* (also called the critical height) as follows:

$$\frac{1}{2}\overline{V}\_h^2(h\_\varepsilon) = \int\_{h\_\varepsilon}^{h\_m} N^2(z)(h\_m - z)dz,\tag{5}$$

where Heinze et al. [6] define *Vh* as the mean horizontal velocity as a function of the horizontal velocity components *<sup>u</sup>* and *<sup>v</sup>*: *Vh* <sup>=</sup> <sup>√</sup> *u*<sup>2</sup> + *v*2. It is assumed in Equation (5) that all the kinetic energy of the flow is converted into potential energy [6]. As was done in [6], Equation (5) will be solved iteratively for *hc* via the trapezoidal rule to numerically compute the integral.

As noted by Heinze et al. [6], if Fr > 1, all fluid parcels, even near the surface, can flow over the obstacle. On the other hand, when Fr < 1, the fluid may stagnate in the upstream region below a certain critical height *hc*, passing primarily around the orographic obstacle rather than over it (flow splitting phenomenon). Above *hc*, the flow would be able to ascend the obstacle [6]. An illustration of flow splitting and of the critical height *hc* is provided in Figure 3.

So, when flow splitting occurs, an expression for the Froude number in terms of the dividing-streamline height may then be obtained, as is noted by Trombetti and Tampieri [36]:

$$\text{Fr}\_{\text{h}} = 1 - \frac{h\_{\text{c}}}{h\_{\text{m}}}.\tag{6}$$

According to Wooldridge et al. [37], Equation (6) may also be used to find the critical height *hc* if Frh is replaced by Fr (defined by Equation (3)), which is applicable when Fr < 1.

**Figure 3.** Theoretical illustration of the flow splitting phenomenon with green arrows depicting airflow. Potential temperature profile is shown in blue. The critical and inversion heights, *hc* and *zi* respectively, are also illustrated.

#### 2.4.3. Inversion Froude Number (Fri)

The inversion Froude number may be computed in the presence of an inversion by adopting Vosper's [38] definition:

$$\mathbf{F}\mathbf{r}\_{i} = \frac{\mathcal{U}\_{nm}}{\sqrt{\mathbf{g}'\mathbf{z}\_{i}}}, \ \mathbf{g}' = \mathbf{g}\frac{\Delta\theta}{\theta\_{m}}.\tag{7}$$

where *zi* is the inversion height, Δθ is the magnitude of the potential temperature difference across the inversion, θ*<sup>m</sup>* is the mean potential temperature below the inversion base and *Unm* is the cross-mountain wind speed (averaged below *zi*). See Figure 3 for a visual representation of some of these parameters.

In the presence of a temperature inversion, several authors [5,38,39] provide a regime diagram for a two-layer flow over mountains as a function of the inversion Froude number (Equation (7)) and of the ratio hm zi , as is shown in Figure A1b from Appendix A.

#### *2.5. Turbulence Indicators*

Five turbulence indicators commonly used in aviation meteorology applications [40,41] are used in this study: Brown, Ellrod TI1, Ellrod TI2, CAT1 and EDR. These parameters were calculated using AROME forecasts.

#### 2.5.1. Brown Index Φ

According to Gill and Buchanan [40] and Sharman et al. [41], the Brown Index Φ may be defined as:

$$\Phi = \sqrt{0.3\zeta\_a^2 + \left(\frac{\partial v}{\partial \mathbf{x}} + \frac{\partial u}{\partial y}\right)^2 + \left(\frac{\partial u}{\partial \mathbf{x}} - \frac{\partial v}{\partial y}\right)^2}, \quad \zeta\_a = \left(\frac{\partial v}{\partial \mathbf{x}} - \frac{\partial u}{\partial y}\right) + f \tag{8}$$

with *u* and *v* denoting the zonal and meridional components of the wind respectively. The vertical component of the absolute vorticity is represented by ζ*a*, which is the sum of the vertical component of the relative vorticity and the Coriolis frequency *f*.

The form applied in this work, Φε, was based on Equation (8) and is expressed as energy dissipation in the following manner:

$$\Phi\_{\ell} = \frac{1}{24} \Phi S\_{V'}^2 \quad S\_V = \sqrt{\left(\frac{\partial u}{\partial z}\right)^2 + \left(\frac{\partial v}{\partial z}\right)^2} \tag{9}$$

where *SV* is the vertical wind shear.

#### 2.5.2. Ellrod TI Indexes

As is mentioned in [40], the Ellrod indicators can be used to quantify turbulence. The first version, Ellrod TI1, is defined as:

$$\text{Ellrod T1} = S\_{\text{V}} \sqrt{\left(\frac{\partial u}{\partial \mathbf{x}} - \frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial v}{\partial \mathbf{x}} + \frac{\partial u}{\partial y}\right)^2},\tag{10}$$

where *Sv* is the vertical wind shear. The Ellrod TI2 index is similar to the Ellrod TI1, except that it incorporates a convergence term in the form of − ∂*u* <sup>∂</sup>*<sup>x</sup>* <sup>+</sup> <sup>∂</sup>*<sup>v</sup>* ∂*y* , as is shown in Equation (11):

$$\text{Ellrod T12} = \text{S}\_{\text{V}} \left[ \sqrt{\left(\frac{\partial u}{\partial \mathbf{x}} - \frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial v}{\partial \mathbf{x}} + \frac{\partial u}{\partial y}\right)^2} - \left(\frac{\partial u}{\partial \mathbf{x}} + \frac{\partial v}{\partial y}\right) \right]. \tag{11}$$

#### 2.5.3. CAT1 Indicator

The CAT1 turbulence indicator, as is denoted in this study, is referred to as "MOS (Model Output Statistics) probability forecast" by Sharman et al. [41] and is defined as:

$$\text{CAT1} = V\_h \sqrt{\left(\frac{\partial u}{\partial \mathbf{x}} - \frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial v}{\partial \mathbf{x}} + \frac{\partial u}{\partial y}\right)^2},\tag{12}$$

where *Vh* is the horizontal wind speed.

#### 2.5.4. EDR Indicator

The Eddy Dissipation Rate, ε, may also be used as a quantifier for turbulence. The method of calculation for this parameter in the current study followed the empirical formula delineated in Frech et al. [42]:

$$\kappa = \frac{\text{TKE}^{1.5}}{L\_{\text{eff}}},\tag{13}$$

where *Le* = 311 m, which is the chosen length scale [42]. TKE represents the turbulent kinetic energy, which is a prognostic variable from the AROME model. The EDR indicator itself, however, was defined as the cube root of ε according to ICAO [43]:

$$\text{EDR} = \varepsilon^{\frac{1}{3}}.\tag{14}$$

This study presents EDR because it is the standard measure to quantify turbulence in aviation applications. Nevertheless, for other scientific communities, it is relevant to show TKE forecasts, which are presented in the supplementary materials.

#### *2.6. RMSE and MAE for Wind Speed and Direction*

The quantitative comparison between observed (from aircraft) and forecast (from AROME) wind data is established using the root mean square error (RMSE) and the mean absolute error (MAE). The calculation of the wind speed error, RMSEWSPD, is performed following the method employed by Grubiši´c et al. [5] for a dataset with *n* points, as follows:

$$\text{RMSE}\_{\text{WSPD}} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left[ V\_h(i)\_{obs} - V\_h(i)\_{sim} \right]^2}. \tag{15}$$

MAEWSPD is similarly defined as follows:

$$\text{MAE}\_{\text{WSPD}} = \frac{1}{n} \sum\_{i=1}^{n} \left| V\_h(i)\_{\text{obs}} - V\_h(i)\_{\text{sim}} \right|. \tag{16}$$

The root mean square error and the mean absolute error for the wind direction, RMSEWDIR and MAEWDIR respectively, are determined using the method described in Jiménez et al. [44]:

$$\text{RMSE}\_{\text{WDIR}} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left[ \Delta\_{\text{WDIR}}(i) \right]^2},\tag{17}$$

$$\text{MAE}\_{\text{WDIR}} = \frac{1}{n} \sum\_{i=1}^{n} |\Delta\_{\text{WDIR}}(i)|\_{\text{\textquotedblleft}} \tag{18}$$

where ΔWDIR(*i*) in degrees is given by:

$$\Delta\_{\text{WDIR}}(i) = \begin{pmatrix} \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i); & \left| \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i) \right| \le 180\\ \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i) - 360; & \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i) > 180\\ \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i) + 360; & \text{WDIR}\_{\text{sim}}(i) - \text{WDIR}\_{\text{obs}}(i) < -180 \end{pmatrix} \tag{19}$$

The comparison between model and observation was done according to the following two-step approach in order to find corresponding values for a given airborne observation:


#### **3. Results and Discussion**

This section begins by illustrating results based on aerial observations. The wind shear intensity factor "I" is then analyzed. Next, outcomes based primarily on AROME forecasts like Froude numbers and turbulence indicators are discussed. Finally, the accuracy of AROME wind predictions is assessed by comparing them with data that was collected by SATA aircraft, using objective scores.

#### *3.1. Observations*

#### 3.1.1. Aerodrome Wind Data

At the LPPI aerodrome, wind was predominantly from the southwest at two instants close to the time of both flights (Table 1). A headwind (*WHD*) of 12/13 kt (6/7 ms−1) and a crosswind (*WCH*) of 14/16 kt (7/8 ms<sup>−</sup>1) were registered on the day of Flight 1. In regard to gusts, the aircraft experienced a headwind (*GWHD*) of 18/20 kt (9/10 ms<sup>−</sup>1) and a crosswind (*GWCH*) of 21/24 kt (11/12 ms<sup>−</sup>1) incoming from the left.


**Table 1.** Wind data (velocity in knots) based on METAR (Meteorological Aerodrome Reports) at Pico aerodrome at the time of flights. The crosswind (*WCH*) is from the left if negative. Negative values of *WHD* indicate headwinds. T1 and T2 are approximate landing times for Flights 1 and 2.

#### 3.1.2. Characterization of Wind Profiles

The wind speed profile based on airborne observations along its path (see Figure 2 and Figure S1) is presented in Figure 4, for Flights 1 and 2. During Flight 1, the aircraft experienced strong winds oscillating between 10 kt (5 ms<sup>−</sup>1) and nearly 40 kt (21 ms−1). At 440 ft (107 m), for instance, the aircraft experienced a wind of 38 kt (19.5 ms<sup>−</sup>1). During Flight 2, in the approach phase, the wind speed reaches its maximum of nearly 22 kt (11 ms−1) at 600 ft (183 m), at a distance of 1.67 NM of the aerodrome. Closer to the aerodrome, the wind decreases to 5 kt (3 ms<sup>−</sup>1) near the surface. Above 1000 ft (305 m) at a distance of nearly 4 NM of the aerodrome, the wind is weak (below 8 kt or 4 ms), due to the wake in the lee of Pico Mountain, as will be shown in the next sections. During both flights, the wind was predominantly from the southwest below 1500 ft or 457 m (Figure S4).

**Figure 4.** Variation of wind speed with altitude along the flight path, below 5000 ft (1524 m), for: (**a**) Flight 1; (**b**) Flight 2.

#### 3.1.3. Wind Shear Intensity Factor "I"

Figure 5a,b detail the relative frequency distributions in both flights for altitudes only below 2100 ft (640 m) in order to focus on the effects of low-level wind shear as is suggested in [25]. During Flight 1, about 36% of the approach phase was carried out under "severe" wind shear conditions and over 15% under "strong". In other words, more than 50% of it experienced "strong" to "severe" wind shear, which is consistent with the occurrence of the hard landing incident. During Flight 2, "severe" wind shear conditions were absent, with only 4% being classified as "strong" and nearly 29% as "moderate" (Figure 5b). Table 2 lists wind shear levels of the first three lowest levels as functions of *dV dt* and <sup>Δ</sup>*<sup>V</sup> <sup>V</sup>* for both flights. By comparing them, it is clear that within the three closest levels to the ground, the hard landing incident (Flight 1) exclusively experiences "severe" wind shear, whereas in Flight 2, a "light" intensity is prevalent.

**Figure 5.** Relative frequency per wind shear intensity category for: (**a**) Flight 1; and (**b**) Flight 2.


#### *3.2. Synoptic Analysis*

During Flight 1, the flow in Pico Island was influenced by the presence of a depression located northwest of the Azores and by a subtropical anticyclone to its west/southeast. Consequently, southwesterly winds of around 20 kt (10 ms−1) were visible south of the island, according to an ECMWF analysis (Figure 6a). For Flight 2, an anticyclone was present south of the Azores, causing southwesterly winds of approximately 10 kt (5 ms<sup>−</sup>1) south of Pico Island (Figure 6b).

#### *3.3. Mesoscale Analysis*

During Flight 1, according to AROME forecasts, the upstream conditions are characterized by a stable stratification, with *N* varying from 0.012 to 0.015 s−<sup>1</sup> and by the absence of an inversion (Figure S5a). The Froude number varies between 0.7 and 0.9 below 400 m, reaching a maximum value of 1.2 at about 1600 m (Figure 7), which coincides with a wind speed of nearly 18 ms−<sup>1</sup> (Figure S5c). Above this level, *Fr* varies between 1 and 1.2. Consistently, *Fr* and *Frh* are 0.95 and 0.9, respectively (Table 3). Therefore, according to previous studies [2,4], no wave breaking, upstream blocking or lee vortex formation are expected (see Figure A1a). The critical height defined in Section 2 is 107 m or 60 m, depending on the method. Consequently, it is expected that air will flow nearly horizontally around the mountain below this level, but will overcome the obstacle above it. Moreover, the parameter *FL* is always below 0.08 (Figure S5c) and so, according to the analysis of *FL* and the Froude number, vertically propagating waves are expected [34].

**Figure 6.** Mean sea level pressure (in hPa) and 10 m wind (kt) based on ECMWF analysis for: (**a**) Flight 1; (**b**) Flight 2.

**Table 3.** Froude number and related results by flight for Point C in Figure 2 with *hc*1, *hc*2, *zi* in meters.


**Figure 7.** Vertical profile of Froude number (from AROME forecasts) for both flights.

During Flight 2, the upstream conditions are characterized by the presence of a short inversion close to 3000 ft (872 m), where *N* varies between 0.015 and 0.018 s−<sup>1</sup> (Figure S5b). The Froude number has values smaller than 0.6 below 2 km (Figure 7). Consistently, *Fr* and *Frh* are 0.42 and 0.31, respectively (Table 3). According to the flow regime diagram presented in Figure A1b, flow splitting and lee vortex formation with reversed flow are expected.

#### 3.3.1. Characterization of Orographic Flow during Flight 1

The horizontal wind and vertical component of the absolute vorticity predicted by the AROME model are depicted in Figure 8 for Flight 1. At 30 m, the flow is partially deflected laterally (mostly on the west side), which is denounced by the curvature of the wind vectors. This deflection is still present at 200 m but is barely visible at 300 m (not shown), while at 500 m (Figure 8b) it is no longer discernible at all. This implies that the critical altitude (*hc*) should lie somewhere between 200 and 300 m. From Table 3, it may be seen that *hc*<sup>1</sup> = 107 m and *hc*<sup>2</sup> = 60 m, revealing that the estimate based on Equation (5) is more accurate. However, it underestimates *hc*, possibly due to an insufficient vertical resolution of the AROME model.

Until 2800 ft (853 m), a pair of positive and negative vertical vorticities is visible over the mountain and on the leeward region, illustrated in Figure 8a for a height of 30 m. However, no signs of reversed flow or vortex formation are present, which is consistent with the results of Bauer et al. [2] and Epifanio [4].

The horizontal wind speed at 680 ft (207 m) shows a stagnation area on the upstream region of Pico Mountain and a maximum wind speed (>40 kt or 21 ms−1) over it with weak winds on the leeward side (Figure 9a). This pattern is coherent with previous studies like in Bauer et al. [2] with their Figure 2.a2. Near Pico Aerodrome, AROME forecast winds with speeds above 35 kt (18 ms−1). Ergo, during the approach to the aerodrome, the aircraft passes a region of elevated wind shear. Figure 9b presents the maximum and minimum vertical velocities in the layer from the surface up to 6000 ft (1829 m), illustrating a wave pattern over Pico and Faial Islands. Figure 10a displays the

vertical cross-section through the flight path in the approach phase (line A–B in Figure 2), showing that AROME forecasts the upward vertical velocity in the order of 100 ft/min (0.5 ms<sup>−</sup>1) near 2000 ft (610 m) and a downward vertical velocity of 100–150 ft/min (0.50–0.76 ms<sup>−</sup>1) near the aerodrome. This is associated with pronounced vertically propagating waves, which are illustrated in Figure 10b (up to 10000 ft or 3048 m) and in Figure S6a (up to 42000 ft or 12802 m), where the highest values exceed 400 ft/min (2 ms<sup>−</sup>1), and then propagate to the vicinity of the aerodrome.

**Figure 8.** AROME absolute vorticity and wind vectors for Flight 1 at an altitude of: (**a**) 30 m; (**b**) 500 m. Dark arrows represent wind vectors.

**Figure 9.** (**a**) Horizontal wind at 680 ft (207 m) for Flight 1; (**b**) Maximum vertical velocity under 6000 ft (1829 m) for Flight 1. Both fields are from AROME forecasts. Shading in (**a**) represents the wind speed in kt and the pink dot identifies the approximate aircraft position at this level. Pink dashed line identifies the approximate aircraft trajectory. Shading in (**b**) represents intensity in ft/min, positive values indicate upward vertical velocity and negative values represent downward vertical velocity. The pink dots are the points C and D shown in Figure 2, which are the endpoints of the cross-section (dashed pink line) used in the next figures.

During Flight 1, the mean wind speed at point C (from Figure 2) is 12.76 ms−<sup>1</sup> (*U*) and the average value of the buoyancy frequency is 0.0129 s−<sup>1</sup> (*N*). Since the horizontal scale (*L*) for Pico Mountain is of the order of 21 km, the condition for vertically propagating gravity waves (*N*<sup>2</sup> > *U*2*k*2) applies [1,29]. See also Figure 13.3 in Cassano and Lynch [29], which illustrates the alternating pattern between upward and downward motion associated with vertically propagating gravity waves.

**Figure 10.** Vertical cross-section of the vertical velocity of AROME through (**a**) line A–B on Figure 2 and (**b**) line C–D on Figure 2, for Flight 1. Shading represents intensity in ft/min, positive values indicate upward vertical velocity and negative values represent downward vertical velocity. The altitude scale of (a) is half of the corresponding scale of (**b**).

#### 3.3.2. Characterization of Orographic Flow during Flight 2

The horizontal wind and vertical component of the absolute vorticity predicted by the AROME model at several heights are presented in Figure 11 for Flight 2. Figure 11a depicts a stagnation area in the windward direction of Pico Mountain and the occurrence of a flow splitting phenomenon in addition to a wake with reversed flow in the lee of the hill. This wake (at 40 m) consists of a pair of vortices, one with a clockwise rotation and negative absolute vorticity and another with a counter-clockwise rotation and positive vorticity. The wake structure exemplified in Figure 11a, which is visible below 600 m (not shown), resembles the results obtained by Smith and Grubiši´c [45] (p. 3747) in their study of Hawaii's wake, which is described as comprising "two large, nearly steady vortices". A similar structure with a pair of counter-rotating vortices downstream of the obstacle was also documented by Smolarkiewicz and Rotunno [7] (see their Figure 1c) and by Epifanio [4] (their Figure 3c). Moreover, the regime flow diagram of Figure A1b in Appendix A, which was used by Schär and Smith [39] and Grubiši´c et al. [5], accurately estimates this situation in which a wake displays reversed flow.

At higher levels, between 600 m and 850 m (Figure 11b,c), no lee vortices are present and the dipole of positive and negative vorticities is weaker and shifted to the east relative to lower levels, mainly at 850 m. At 1100 m, the incoming flow is mainly from the west and overcomes the mountain (Figure 11d). Finally, the analysis of Figure 11 suggests that the critical height in Flight 2 should be approximately 850 m, which is close to the estimate of the dividing streamline concept (see Table 3).

The horizontal wind speed at 40 m is depicted in Figure 12a, which shows that near the aircraft position, the wind speed varies between 15 and 20 kt (8–10 ms<sup>−</sup>1). It is also visible from this figure that the aircraft during the approach phase (pink dashed line in Figure 12a) has passed through the wake. Figures 12a and 13 also show signs of a vertically propagating mountain wave, however, less energetic than the one seen on the day of Flight 1 (compared to Figure 10) and propagating to lower levels (up to 25000 ft or 7620 m), as is illustrated in Figure S6b. During Flight 2, the maximum vertical velocity on the leeward side of Pico Mountain (Figure 13b) is less than 150 ft/min (0.8 ms<sup>−</sup>1), while during Flight 1 the maximum vertical velocity was above 400 ft/min (2 ms<sup>−</sup>1). Upon propagation of this wave to the aerodrome (Figure 13a), it is seen that only a weak downdraft (less than 100 ft/min or 0.5 ms−1) lies atop it with no hints of updrafts along the flight path. In Figure 13b, the reversed flow associated with the lee vortices is also visible until about 2000 ft (610 m).

*Atmosphere* **2019**, *10*, 350

**Figure 11.** AROME absolute vorticity and wind vectors for Flight 2 at an altitude of: (**a**) 40 m; (**b**) 640 m; (**c**) 850 m; (**d**) 1100 m. Dark arrows represent wind vectors.

**Figure 12.** (**a**) Horizontal wind at 135 ft (40 m) for Flight 2.; (**b**) Maximum vertical velocity component under 6000 ft (1829 m) for Flight 2. Both fields are from AROME forecasts. Shading in (**a**) represents the wind speed in kt and the pink dot identifies the approximate aircraft position at this level. Pink dashed line identifies the approximate aircraft trajectory. Shading in (**b**) represents intensity in ft/min, positive values indicate upward vertical velocity and negative values represent downward vertical velocity. The pink dots are the points C and D shown in Figure 2, which are the endpoints of cross-section (dashed pink line) used in the next figures.

**Figure 13.** Vertical cross-section of the vertical velocity of AROME through (**a**) line A-B in Figure 2 and (**b**) line C–D in Figure 2, for Flight 2. Shading represents intensity in ft/min, positive values indicate upward vertical velocity and negative values represent downward vertical velocity. The altitude scale of (**a**) is half of the corresponding scale of (**b**).

#### *3.4. Turbulence Indicators*

Turbulence indicators Ellrod TI2 and CAT1 at 66 m are presented in Figure 14a,b, for Flight 1. Both indicators show maximum values leeward of Pico mountain, with TI2 exceeding 18 <sup>×</sup> 10−<sup>6</sup> s−<sup>2</sup> and CAT1 exceeding 8 <sup>×</sup> 10−<sup>3</sup> ms−<sup>2</sup> in the proximity of the aerodrome where the incident occurred. The cross-section along the flight path shows that for a distance less than 5 NM (9.3 km) of the runway, between 300 and 1600 ft, TI2 and CAT1 reach values above 18 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>s</sup>−<sup>2</sup> and 8 <sup>×</sup> <sup>10</sup>−<sup>3</sup> ms<sup>−</sup>2, respectively (Figure 14c,d). Comparing the indicators with the wind shear Intensity Factor "I" (Figure 15a,b), the aircraft experienced severe wind shear at these levels, suggesting that these values for TI2 and CAT1 could be tested as thresholds indicative of severe wind shear conditions. Ellrod TI1 (not shown in Figure 14) has a similar spatial pattern to TI2, but with lower values (Figure 15a). The Brown index is maximum over land near the LPPI aerodrome (not shown in Figure 14), exceeding 20 <sup>×</sup> 10−<sup>9</sup> s−<sup>3</sup> near the surface and decreasing drastically with height, reaching a value of 13 <sup>×</sup> 10−<sup>9</sup> s−<sup>3</sup> at 600 ft (Figure 15c).

The EDR indicator shows a different pattern when compared to others. EDR shows higher values over land than over sea, depicting a maximum value above 0.4 m2/<sup>3</sup> s−<sup>1</sup> near the top of Pico Mountain (Figure 16a). The higher values over land are consistent with higher roughness lengths (than over sea), which contribute to a higher mechanical production of turbulence [16]. At 66 m, EDR reaches a value between 0.35 and 0.4 m2/<sup>3</sup> s−<sup>1</sup> near the aerodrome (corresponding to moderate turbulence). Figure 16b shows that along the flight path, EDR decreases with height, varying from 0.35–0.4 m2/<sup>3</sup> s−<sup>1</sup> near the surface to 0.15–0.2 m2/<sup>3</sup> s−<sup>1</sup> at 2000 ft (610 m). It is interesting to mention that the thresholds applied in Figure 16 are lower than those presented in ICAO [43] but are those proposed recently [46], after a study by Sharman et al. [47] that showed that the thresholds in reference [43] were overestimated. According to reference [46], turbulence is severe when EDR <sup>≥</sup> 0.45 m2/<sup>3</sup> s−<sup>1</sup> and is moderate for EDR between 0.2 and 0.45 m2/<sup>3</sup> s<sup>−</sup>1. Comparing EDR to the wind shear Intensity Factor "I" (Figure 15d), it is clear that in the presence of severe wind shear conditions, EDR has values between 0.2 and 0.25 m2/<sup>3</sup> s−<sup>1</sup> above 700 ft. For heights under 200 ft (61 m), EDR is close to 0.3, which suggests an underestimation of the severe wind shear conditions.

**Figure 14.** Horizontal variation of (**a**) Ellrod TI2 and (**b**) CAT1 at 66 m for Flight 1. Dotted pink line shows aircraft trajectory and pink bar locates the aerodrome. The circled "x" mark identifies the approximate aircraft position at 66 m. Cross-sectional variation of (**c**) Ellrod TI2 and (**d**) CAT1.

**Figure 15.** AROME turbulence indicators (bottom x-axis) vs. the wind shear intensity factor "I" (top x-axis), along the flight path during Flight 1. Ellrod TI1, TI2 and Brown indicators are shown in (**a**) and (**c**), respectively. Solid red lines in (**b**) and (**d**) indicate CAT1 and EDR values, respectively.

**Figure 16.** (**a**) Horizontal representation of EDR during Flight 1. The circled "X" mark identifies aircraft position at 66 m. **(b**) Cross-sectional representation of EDR.

During Flight 2, the index Ellrod TI2 shows maximum values in the leeward region of the mountain and exhibits a value of 9 <sup>×</sup> 10−<sup>6</sup> s−<sup>2</sup> near the aircraft position (Figure 17a). The index CAT1 shows maximum values above 5 <sup>×</sup> 10–3 ms−<sup>2</sup> downstream of Pico Aerodrome, close to the flight path (Figure 17b). Along the flight path itself, Ellrod TI2 reaches values above 9 <sup>×</sup> 10−<sup>6</sup> s−<sup>2</sup> below 500 ft (152 m), while CAT1 attains values between 4 <sup>×</sup> 10−<sup>3</sup> ms−<sup>2</sup> and 8 <sup>×</sup> 10−<sup>3</sup> ms−<sup>2</sup> below 900 ft or 274 m (Figure 17c,d and Figure S7), which are predominantly "moderate" wind shear conditions. The TI1 indicator is very similar to TI2 and the Brown index decreases again drastically with height, exhibiting values below 4 <sup>×</sup>10−<sup>9</sup> <sup>s</sup>−<sup>3</sup> above 300 ft (Figure S7).

**Figure 17.** Horizontal variation of (**a**) Ellrod TI2 and (**b**) CAT1 at 66 m for Flight 2. Dotted pink line shows the aircraft trajectory and pink bar indicates the aerodrome. The circled "X" mark identifies the approximate aircraft position at 66 m. Cross-sectional variation of (**c**) Ellrod TI2 and (**d**) CAT1.

Analysis of EDR depicts the occurrence of "light" turbulence near the aerodrome for Flight 2 (Figure 18). The "I" factor is mostly consistent with this categorization as wind shear conditions reach intensities that are predominantly "light" for heights between 200–500 ft (61–152 m). Nevertheless, in some instants the wind shear was classified as "moderate", while the EDR values above 400 ft are below 0.1 m2/<sup>3</sup> s−<sup>1</sup> suggesting the absence of turbulence or wind shear (Figure S7d).

**Figure 18.** (**a**) Horizontal representation of EDR during Flight 2. The circled "X" mark identifies approximate aircraft position at 66 m. (**b**) Cross-sectional representation of EDR.

This preliminary scrutiny of turbulence indicators suggests that CAT1, Ellrod TI2 and EDR could be applied together to arrive at a classification for the hazard level of turbulence. Moreover, preliminary thresholds for these indicators are proposed (Table A1), which should be subject to further validation, using more flight data.

As expected from their similarity in mathematical construction, both TKE and EDR evolve analogously during both flights (Figure S8a,b). In Flight 1, when the EDR indicator attains values of around 0.3 m2/<sup>3</sup> s−<sup>1</sup> (coinciding with "severe" wind shear conditions according to the "I" factor), TKE has values close to 4 m2s−2. During Flight 2, an EDR hovering near 0.1 m2/<sup>3</sup> s−<sup>1</sup> translated to a TKE value of about 0.5 m2s−2.

In order to certify that the turbulence was indeed predominantly greater for Flight 1 than for Flight 2, independently of the discrepant mean flow speeds in each case, the dimensionless turbulence intensity parameter (*TI* <sup>=</sup> <sup>√</sup> *TKE*/*Vh*), defined by Stull [16], was also plotted (Figure S9). It is clear from these results that the *TI* of Flight 1, with two maxima slightly above 0.25 (near 200 ft or 61 m and 1200 ft or 366 m), mostly exceeds that of Flight 2, with a maximum slightly above 0.15 (near 300 ft or 91 m).

#### *3.5. Objective Verification of Wind Forecasts from AROME*

The accuracy of AROME wind predictions is assessed by comparing them with data that was collected by SATA aircraft during Flights 1 and 2. Absolute errors of wind speed and direction are also mapped for each observation and forecast pair in Figure 19. The color coding for wind magnitude and its direction is based on ICAO guidelines [43] for operationally acceptable levels of forecast errors. For wind speed, an absolute value error of up to 5 kt (2.6 ms<sup>−</sup>1) was deemed operationally desirable. For the direction, this value is 20◦. For Flight 1, there is a good agreement between observations and forecasts.

**Figure 19.** Absolute value differences for wind speed and direction for: (**a**) Flight 1 and (**b**) Flight 2.

Values for the root mean square error (RMSE) and mean absolute error (MAE) regarding wind speed and direction are displayed in Table 4, showing that wind speed RMSE for Flight 1 is only slightly above operational recommendations (about 30% larger) and the RMSE for wind direction is 11.4◦, well within the acceptable limits. The Flight 1 MAE for wind speed shows that the 2.6 ms−<sup>1</sup> recommendation is exceeded by only 7% and the wind direction MAE is again within limits with a value of 8.1◦.

**Table 4.** Summary of RMSE and MAE for wind speed and direction for both flights.


During Flight 2, the RMSE and MAE for wind speed are 4.43 ms−<sup>1</sup> and 3.77 ms−<sup>1</sup> respectively, slightly over 30% larger than Flight 1 errors. For Flight 2 wind direction errors, the RMSE and MAE are around 83.3◦ and 65.4◦. It is important to mention that when the wind speed is above 10 kt, the wind direction errors are inferior to 45◦ (see Figure S10). However, during Flight 2 at altitudes above 1200 ft, the wind is below 10 kt (Figure S4b) and at these levels, large errors (>120◦) are found in the wind direction (Figure 4b, Figure S10). The limitations of NWP models in accurately forecasting the wind direction for weak winds was previously documented by Jiménez et al. [44], who showed that RMSEWDIR increases with decreasing wind speeds, reaching values over 70◦ for speeds under 2 kt (1 ms<sup>−</sup>1). It is also important to note that during Flight 2, a wake developed in the lee of Pico mountain in the vicinity of the flight path. In the wake region, the wind varies strongly over short distances (see Figure 12), implying that a small phase error may cause large errors in wind speed and/or wind direction. According to Ricard et al. [22], the effective resolution of the AROME model is about 9 times the grid spacing (i.e., 22 km), which may be insufficient to correctly resolve the wind variations within the lee wake.

#### **4. Summary and Conclusions**

The present article focuses on the orographic flow associated with Pico Mountain in the Azores in two different synoptic environments, using SATA airborne measurements and AROME forecasts. Aerial measurements were used to compute the wind shear intensity factor "I" according to Guan and Yong [25] and to categorize its severity according to four classes: "light", "moderate", "strong" and "severe".

During the first scenario, a hard landing incident (Flight 1) occurred in Pico's aerodrome. On this day, according to the AROME model, the upstream conditions were characterized by southwest winds with speeds of up to 40 kt (20 ms<sup>−</sup>1) in the ABL and stable stratification, which translated to a Froude number slightly above 1 near the mountain top. These conditions triggered the development of vertically propagating mountain waves, which created an undulation pattern of upward and downward motion downstream of Pico Mountain, leading to maximum vertical velocities above 400 ft/min (2 ms<sup>−</sup>1) and above 200 ft/min in the flight path. Based on ICAO Annex 3 guidelines [43], a mountain wave is deemed "severe" whenever downdrafts exceed 600 ft/min (3 ms<sup>−</sup>1) and "moderate" whenever they vary between 350 to 600 ft/min (1.75–3.0 ms<sup>−</sup>1). However, it is important to note that the AROME model may have underestimated the intensity of the vertical velocities, so this pattern of oscillating updraft/downdraft may have strongly contributed to the hard landing event. Furthermore, large wind speeds (38 kt or 20 ms<sup>−</sup>1) were registered close to 440 ft (107 m) and the wind shear intensity factor "I" based on airborne data showed that 36% of the approach phase (below 2100 ft or 640 m) occurred under "severe" wind shear conditions and 16% under "strong". At the aerodrome, METAR observations showed southwest winds (18–21 kt or 9–11 ms−1) with gusts between 28 and 31 kt (14 and 16 ms<sup>−</sup>1).

During the second flight, the upstream conditions were characterized by southwest winds with speeds below 20 kt (10 ms<sup>−</sup>1) in the PBL and a temperature inversion in the layer between 860 to 1100 m. In this range, the Froude number is less than 0.5. According to previous studies, these conditions promote upstream stagnation and flow splitting with the formation of lee vortices in low levels with reversed flow, which was predicted by the AROME model. AROME forecasts for Flight 2 also showed a pattern of vertically propagating waves, but with smaller vertical velocities and shorter vertical propagation compared to Flight 1. According to airborne observations, Flight 2 occurred under wind shear conditions predominantly rated as "light" (68%) or "moderate" (29%), with its most intense class being "strong" (only during less than 4% of the landing phase).

Moreover, the pattern of several turbulence indicators derived from AROME forecasts was analyzed for both flights. All these indicators forecast higher values for Flight 1 than for Flight 2, which is consistent with the wind shear intensity factor "I" derived from airborne observations. As a result of the comparison between "I" and the AROME turbulence indicators, preliminary thresholds have been proposed for three indexes. In the future, however, more flights should be perused in order to verify the validity of such thresholds. The use of other wind shear and turbulence indicators derived from airborne data should also be tested.

Lastly, this study provides an objective verification of AROME wind forecasts, showing a good agreement with airborne observations for wind speeds above 10 kt, but a poor skill for weaker winds. For Flight 1, characterized by wind speeds over 30 kt (15 ms−1) during several instants, AROME forecasts show a good agreement with observations in terms of wind speed and direction, with an RMSE of 3.4 ms−<sup>1</sup> and 11.4◦, respectively. The MAE score is 2.8 ms−<sup>1</sup> and 8.1◦ for wind speed and direction, respectively. During Flight 2, in the approach phase, the maximum wind speed is 22 kt (11 ms<sup>−</sup>1). In this case, AROME has an RMSE of 4.4 ms−<sup>1</sup> and 83.3◦ for wind speed and wind direction. MAE revealed values of 3.8 ms−<sup>1</sup> and 65.4◦ for the wind speed and direction. In particular, when the wind speed is below 10 kt (5 ms<sup>−</sup>1), large errors (>120◦) in the wind direction can be found. This is a limitation of NWP models that was previously documented by Jiménez et al. [44] for the WRF model.

Strong ageostrophic winds often develop in coastal gaps or channels when an along-gap pressure gradient is created [48]. A gap flow may also contribute to reinforce the winds downstream of the strait between Faial and Pico islands. This hypothesis should be subject to future research.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4433/10/7/350/s1, Figure S1: (a) Enlarged view of Pico Mountain and its Aerodrome (LPPI), (b) Location of Pico Aerodrome (LPPI) relative to Pico Mountain, Figure S2: Cross and headwind components, Figure S3: Boundaries for the classification of wind shear intensity, Figure S4: Wind profile along flight path according to aerial observations and AROME for (a) Flight 1, and (b) Flight 2, Figure S5: Absolute temperature and N profiles at point C for (a) Flight 1 and (b) Flight 2, Figure S6: (a) Vertically propagating mountain waves predicted by AROME during Flight 1, and (b)

Flight 2, Figure S7: AROME turbulence indicators vs. the wind shear intensity factor *I* along the flight path during Flight 2. Figure S8: Variation of Turbulent Kinetic Energy (TKE) and Eddy Dissipation Rate (EDR) with Height for (a): Flight 1 and (b) Flight 2. Figure S9: Turbulence Intensity for both Flights. Figure S10: Absolute Value Wind Direction Errors versus Wind Speed for both Flights.

**Author Contributions:** Conceptualization, J.M., M.B.-P. and P.S.; software, J.M. and M.B.-P.; Analysis, J.M., M.B.-P. and P.S.; writing—original draft preparation, J.M. and M.B.-P.; writing—review and editing, J.M., M.B.-P. and P.S.; Figures, J.M. and M.B.-P.; supervision, M.B.-P. and P.S.

**Funding:** FCT - Portuguese Foundation for Science and Technology, under the project UID/AGR/04033/2019.

**Acknowledgments:** This work is supported by National Funds by FCT - Portuguese Foundation for Science and Technology, under the project UID/AGR/04033/2019. The authors are also thankful to SATA International – Azores Airlines, S.A. for providing the airborne data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** (**a**) Regime flow diagram adapted from Bauer et al. [2] for conditions of Flight 1; (**b**) Regime flow diagram adapted from Grubiši´c et al. [5] for conditions of Flight 2. Red mark indicates predicted states according to results from this study. Pico's orography is approximately circular (unitary horizontal aspect ratio).

**Table A1.** Suggested preliminary thresholds for three turbulence indicators. The thresholds values of EDR are only valid for heights below 500 ft (152 m).


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Hurricane Boundary Layer Height Relative to Storm Motion from GPS Dropsonde Composites**

**Yifang Ren 1, Jun A. Zhang 2,3,\*, Stephen R. Guimond 4,5 and Xiang Wang <sup>6</sup>**


Received: 9 May 2019; Accepted: 10 June 2019; Published: 21 June 2019

**Abstract:** This study investigates the asymmetric distribution of hurricane boundary layer height scales in a storm-motion-relative framework using global positioning system (GPS) dropsonde observations. Data from a total of 1916 dropsondes collected within four times the radius of maximum wind speed of 37 named hurricanes over the Atlantic basin from 1998 to 2015 are analyzed in the composite framework. Motion-relative quadrant mean composite analyses show that both the kinematic and thermodynamic boundary layer height scales tend to increase with increasing radius in all four motion-relative quadrants. It is also found that the thermodynamic mixed layer depth and height of maximum tangential wind speed are within the inflow layer in all motion-relative quadrants. The inflow layer depth and height of the maximum tangential wind are both found to be deeper in the two front quadrants, and they are largest in the right-front quadrant. The difference in the thermodynamic mixed layer depth between the front and back quadrants is smaller than that in the kinematic boundary layer height. The thermodynamic mixed layer is shallowest in the right-rear quadrant, which may be due to the cold wake phenomena. The boundary layer height derived using the critical Richardson number (*Ric*) method shows a similar front-back asymmetry as the kinematic boundary layer height.

**Keywords:** atmospheric boundary layer; tropical cyclone; storm motion; asymmetry; hurricane; aircraft; dropsonde

#### **1. Introduction**

The atmospheric boundary layer (ABL) is the turbulent layer close to the Earth's surface that is influenced by surface friction. Physical processes in the ABL play an important role in regulating the atmosphere at both weather and climate scales, including hurricanes. The top of the ABL (i.e., the boundary layer height) is a key parameter that determines the vertical distribution of turbulent mixing in numerical models that require a boundary layer parameterization. In numerical models where turbulent fluxes are parameterized using a first-order K-profile method, the boundary layer height is usually defined as the height at which the bulk Richardson number (Ric) reaches a threshold (typically 0.25), where Ric is defined as the ratio of buoyancy to shear forcing [1–3]. Note that this Richardson number method was also previously used in observational studies to determine the boundary layer height in non-hurricane conditions [4–8].

Under non-hurricane conditions when vertical profiles of turbulent intensity and/or flux are measured, the boundary layer height is usually taken as the height where the magnitude of the turbulence parameter is much smaller (~95%) than that in the surface layer. Note that the surface layer height is usually taken as 10% of the ABL height [9]. Without turbulence data, the height of the temperature inversion layer is generally used to define the boundary layer height, in particular, under unstable conditions [9–12]. Since passive scalars are accumulated in the ABL, large gradients of temperature and water vapor occur at the inversion layer capping the ABL [13,14]. Above the inversion layer, a rising air parcel typically becomes neutrally buoyant.

The top of the stable or neutral boundary layer is more difficult to determine than the unstable convective boundary layer. The boundary layer height is typically defined using the vertical gradient of virtual potential temperature with a threshold [4,15–17], when flux data are not available. This method has been widely used when radiosonde data are available. The parcel method developed by Holzworth [18] was also used in previous studies to estimate the boundary layer height using radiosonde data [15,19,20]. Remote sensing data (e.g., lidars, radars, sodas, and wind profilers), although sometimes ambiguous, can also be used to determine the boundary layer height when other types of data are not available [21–24].

In hurricane conditions, it is even more difficult to acquire direct turbulence measurements due to safety issues and measurement constraints. Multi-level turbulent intensity and flux data that can be used to estimate the boundary layer height are scarce in hurricane conditions. The only in situ data of this kind was collected in regions either 100 km away from the hurricane eyewall or in weak tropical storms [25–28]. Note that a new airborne Doppler radar dataset from the Imaging Wind and Rain Airborne Profiler (IWRAP) is able to provide three-dimensional winds, deep into the hurricane boundary layer, at high resolution (~200 m horizontal and 30 m vertical), which will be useful for characterizing the boundary layer height using dynamical metrics [29,30].

Dropsonde data have been used to estimate the hurricane boundary layer height scales in previous studies [31]. Zhang et al. [31] pointed out differences between the kinematic boundary layer height denoted by the height of the maximum tangential wind speed (*h*vtmax) or inflow layer depth (*hi*nflow), and the thermodynamic boundary layer height denoted by the mixed layer depth (*z*i). Of note, *hi*nflow is taken as the height of strong inflow layer where the inflow strength is equal to 10% of the maximum inflow [32]. Furthermore, *z*<sup>i</sup> is taken as the height where the vertical gradient of the virtual potential temperature (theta-v) is equal to 3 K km<sup>−</sup>1, a widely used approach in non-hurricane ABL studies [4]. Ming et al. [33] confirmed the result of Zhang et al. [31] by analyzing a smaller number of dropsonde data in Pacific typhoons.

Previous modeling studies on air-sea coupling processes of individual tropical cyclones (TCs) showed an asymmetric distribution of the boundary layer structure relative to the storm motion, suggesting that the boundary layer is more stable in the right-rear quadrant [34,35]. Although a significant number of observational studies have documented the mean boundary layer structure of individual hurricanes [27,36–41], the variability of the ABL structure as a function of azimuth remains to be explored. The objective of the present study is to investigate the asymmetry of the boundary layer structure in a storm-motion-relative framework with a focus on the boundary layer height scales. This study is a first attempt to analyze the motion-induced asymmetry of the boundary layer heights in a climatological sense by compositing dropsonde data collected in multiple hurricanes. The results of this study will be useful for improving our understanding of the low-level structure of hurricanes and providing composite analyses that can be used for model evaluation and physics improvements [42].

This paper is organized as follows. Section 2 describes the dropsonde data and the composite analysis methodology. The dropsonde composites used to compute the boundary layer heights are shown in Section 3, which is followed by the discussion and conclusion in Section 4.

#### **2. Data and Composite Methodology**

Dropsondes are typically deployed from research or reconnaissance aircraft at a height of ~3 km (i.e., 700 hPa), with a descending rate of ~ 10 m s−<sup>1</sup> and a vertical sampling resolution of ~7 m. During the process of descent, atmospheric profiles of pressure, temperature, relative humidity, and horizontal winds are collected by the dropsonde. The accuracy of pressure, temperature, relative humidity and horizontal wind speed are <sup>±</sup>1.0 hPa, <sup>±</sup>0.2 ◦C, <sup>±</sup>5%, and <sup>±</sup>0.5 ms<sup>−</sup>1, respectively [43]. The dropsonde data used in this study are from the Hurricane Research Division (HRD) archives and have been quality controlled using the Editsonde and/or ASPEN programs.

In this study, a total of 7326 GPS dropsonde profiles collected by research aircraft in 37 TCs from 1998 to 2015 are analyzed. The storm intensity is obtained from the hurricane best track database produced by the National Hurricane Center (NHC). The storm tracks processed by HRD, which combine the best track storm locations with aircraft center fixes, are used in the analysis. We linearly interpolate the storm intensity, center and motion from the track data to the dropsonde time. For the composite analysis of the inner-core structure, the dropsonde profiles must meet the following requirements: (i) They have measurements of wind speed, temperature, and humidity from flight-level all the way to the surface with no data gaps, (ii) the maximum sustained wind speed (i.e., storm intensity) from NHC's best track is larger than 64 kt (Cat 1 strength), and (iii) the sondes are deployed within four times the radius of maximum wind speed (RMW). A total of 1916 dropsonde profiles are used in the final composite analysis after this screening. Note that the RMW data used in this study are calculated based on the flight-level data following [44]. Table 1 summarizes the storm information and numbers of dropsondes used in this work.


**Table 1.** Storm information and number of dropsondes (1916).

The dropsonde distribution relative to the storm center is shown in Figure 1, with observation locations rotated with respect to the storm motion direction to the top of the figure. Figure 1 shows nearly evenly distributed data at motion-relative azimuths at all radii, although more dropsonde data are located close to the RMW (i.e., r/RMW = 1). The dropsonde data are composited in four motion-relative quadrants defined clockwise from the storm motion direction, as right-front, right-rear, left-rear and left-front (Figure 2). The storm characteristics in terms of frequency distribution, including the storm intensity *Vmax*, RMW, storm translational direction, and storm speed are presented in Figure 3. Storm intensities range from 65–150 kt, sizes in terms of RMW range from 10–65 km, and translational speed ranges from 2–20 ms−1. The mean storm intensity for the whole sample is 105.5 kt, the mean RMW is 28.5 km, and the mean storm translational speed is 5.27 ms<sup>−</sup>1.

**Figure 1.** Storm-relative two-dimensional distribution of dropsonde surface observation locations. Cross- and along-track positions are normalized by the radius of maximum wind at the time of observation. The arrow indicates the storm motion direction.

**Figure 2.** Radial distribution of dropsonde counts per bin as a function of normalized distance by radius of maximum wind speed (RMW).

In this study, we use the same composite methodology as Zhang et al. [31] to construct the radial-vertical profiles of virtual temperature (θ*v*), tangential wind, radial wind, and Richardson number in four quadrants relative to the storm motion. When compositing the data, the radial bin width is 0.2 *r*∗ (*r*\* = *r*/RMW) for the inner core (*r*∗ < 2), and it is 0.4 *r*∗ for the outer part. The data are also bin averaged vertically at 10 m resolution. The final averaged data are also smoothed with three passes of a 1-2-1 filter, instead of five passes as in Zhang et al. [31]. The data sampling sizes for different quadrants are displayed in Figure 2 as a function of normalized radius and height. The largest sample size is located in the vicinity of the eyewall as expected. Figure 2 also indicates that the data samples for the four motion-relative quadrants are similar.

**Figure 3.** Frequency distribution of dropsondes according to the corresponding (**a**) storm intensity, (**b**) radius of maximum wind speed (RMW), (**c**) storm speed, and (**d**) storm direction rotated clockwise with 0◦ pointing to the north.

#### **3. Results**

The storm-motion-relative normalized radius-height representation of the tangential wind speed is displayed in Figure 4. The tangential wind first increases with height to a maximum value then decreases with height in all four motion-relative quadrants. This tangential wind maximum has been called the boundary layer jet [45,46]. The maximum tangential wind speed is largest in the right-front quadrant and it is smallest at the right-rear quadrant. Specifically, the peak values of the tangential wind speed are 53.5, 55.8, 52.6, and 44.3 ms−<sup>1</sup> for left-front, right-front, left-rear and right-rear quadrant, respectively. The fact that the front quadrants have stronger tangential wind speeds than the rear quadrants in a storm-relative framework is consistent with previous theoretical studies [45–48].

It is also evident from Figure 4 that the jet height, *h*vtmax, has a trend of decreasing with smaller r\* values in all four motion-relative quadrants. This result further supports the findings of Zhang et al. [31] in terms of the radial variation of the hurricane boundary layer height, which may be applied in the asymmetric sense. Interestingly, *h*vtmax is higher in the front quadrants than in the rear quadrants, which may be associated with the stronger tangential wind speed in the front quadrants. Specifically, *h*vtmax increases from ~600 m at a radius of *r*<sup>∗</sup> = 1 to ~1200 m at a radius of *r*<sup>∗</sup> = 2, then slightly increases with increasing radius in the front two quadrants. In the left-rear quadrant, *h*vtmax is ~100 m smaller than in the front two quadrants at r\* = 1, and increases with radius to ~800 m at r\* = 2, then increases to ~900 m at r\* = 4. In the right-rear quadrant, *h*vtmax is the lowest among all four quadrants, and it

gradually increases from ~400 m at radius of *r*∗ = 1 to ~600 m at radius of *r*∗ = 2, and increases to ~800 m at r\* = 4.

The normalized radius-height representation of the radial velocity for the four motion-relative quadrants is displayed in Figure 5. The peak values of the radial wind speed are −18.5, −18.7, −13.5 and <sup>−</sup>15.4 ms−<sup>1</sup> for left-front, right-front, left-rear and right-rear quadrant, respectively. These peak inflow values are located at ~100 m altitude between r\* = 1 − 2 and closer to r\* = 1, which is consistent with the axisymmetric structure documented by Zhang et al. [31]. In addition, a pronounced outflow of 5–10 m s−<sup>1</sup> can be seen in all the four quadrants above the inflow layer. The inflow is stronger in the front quadrants compared to the rear quadrants, which is consistent with previous theoretical and numerical studies [46–48]. Observational studies of individual hurricanes also showed similar front-back wind asymmetry [37,38,49].

The inflow layer depth, *hi*nflow, which is depicted by the white line in Figure 5, shows a decreasing trend with decreasing radius in all four motion-relative quadrants. It is also evident from Figure 5 that, *hi*nflow is larger in the front two quadrants than in the rear two quadrants. Specifically, *hi*nflow evolves similarly in left-front, right-front and right-rear quadrants from an altitude of ~600 m to an altitude of 900–1200 m with highest heights in the right-front quadrant and lowest heights in the right-rear quadrant. However, *hi*nflow remains nearly constant (~500 m) with a slightly negative trend in the left-rear quadrant. This result suggests that the asymmetric distribution of *hi*nflow follows that of *h*vtmax. From Figures 4 and 5, the maximum tangential wind speed in the eyewall region is close to the top of the inflow layer in all quadrants.

The thermodynamic mixed layer depth, zi, is depicted by the white line in Figure 6 that shows the normalized radius-height representation of the vertical gradient of theta-v (*d*θ*v*/*dz*). In all quadrants, the boundary layer is unstable near the surface as indicated by the negative value of *d*θ*v*/*dz*. Above this very shallow unstable layer, the boundary layer becomes nearly neutral up to the mixed layer depth and then becomes stable. There is a strong stable layer inside the RMW above the mixed layer at heights of 600–2000 m with *d*θ*v*/*dz* > 5 K/km. This strong stable layer is shallowest at the right-rear quadrant among the four motion-relative quadrants, which may be tied to the cold wake phenomena observed in the sea surface temperature (SST) field [50–52]. On average when the storm motion of a hurricane is ~6 ms−1, which is the average value of our data, the typical SST close to the cold wake region at the right-rear quadrant is of the order of 1–2 K smaller than the front quadrants based on in-situ data [53]. The smaller SST in the right-rear quadrant would stabilize the low-level boundary layer due to a reduction in surface enthalpy fluxes.

In all four quadrants, zi decreases with distance toward the storm center, in a similar manner as the kinematic boundary layer heights. This result again supports that of Zhang et al. [31] and Zhang et al. [54]. From Figures 4–6, it is clear that zi is much smaller than *hi*nflow and *h*vtmax in all quadrants, the kinematic and thermodynamic boundary layer heights largely depart from each other, as noted by Zhang et al. [31]. This structure is different from that of the ABL in non-hurricane conditions. Close to the eyewall region (r\* < 1.5), zi is nearly symmetric with a value of ~200 m. Outside the eyewall, the left-front quadrant has the largest zi while the other three quadrants have similar magnitudes of zi. In the outer radii (r\* > 3), zi is largest in the left-front quadrant and is smallest in the right-rear quadrant. There is a weak front-back asymmetry in zi, which is similar to the asymmetric distribution of the kinematic boundary layer height, but the mixed layer depth difference is less than ~100 m.

The storm-motion-relative normalized radius-height representation of the bulk Richardson number is displayed in Figure 7. Here the height of *Ric* = 0.25 is taken as the top of the boundary layer (hRic), which is depicted by the solid white line in Figure 7. The front-back difference in hRic is clearly shown in Figure 7 with the front two quadrants displaying a deeper boundary layer, consistent with the other height scales investigated earlier. Combining Figures 4–7, it appears that hRic lies between the thermodynamic mixed layer depth and the kinematic boundary layer height, which agrees with the axisymmetric structure documented by Zhang et al. [31]. The left-front quadrant has the deepest boundary layer while the right-rear quadrant has the shallowest boundary layer, in terms of hRic.

**Figure 4.** Composite analysis result of the relative tangential wind velocity as a function of altitude and the normalized radius to the storm center for the four quadrants relative to the motion direction. The panels show the left-front (**a**), right-front (**b**), left-rear (**c**) and right-rear (**d**) quadrants. The white dashed line in each panel depicts the height of the maximum tangential wind speed varying with radius.

**Figure 5.** Same as in Figure 4 but the results are for the relative radial wind velocity as a function of altitude and the normalized radius to the storm center for the four quadrants relative to the shear direction. The panels show the left-front (**a**), right-front (**b**), left-rear (**c**) and right-rear (**d**) quadrants. The white line in each panel represents the height of 10% peak inflow.

**Figure 6.** Same as in Figure 4 but the results show the lapse rate of the virtual potential temperature. The panels show the left-front (**a**), right-front (**b**), left-rear (**c**) and right-rear (**d**) quadrants. The thick white line denotes the contour. The contour denotes the constant contour of *d*θ*v*/*dz* = 3 K km−1.

**Figure 7.** Same as in Figure 4 but the results are for the Richardson numbers as a function of altitude and the normalized radius to the storm center. The panels show the left-front (**a**), right-front (**b**), left-rear (**c**) and right-rear (**d**) quadrants. The white line shows the contour of 0.25.

#### **4. Discussion and Conclusions**

The ABL plays an important role in the energy transport processes related to hurricane intensification and maximum intensity [32,55–60] and it is essential to understand the ABL structure. This paper analyzes a total of 1916 GPS dropsondes within four times the RMW distance from 37 hurricanes over the tropical Atlantic basin from 1998 to 2015 to study the characteristic height scales of the ABL with respect to the storm-motion direction. Figure 8 is a schematic diagram that summarizes the height scales investigated in this study based on the composite dropsonde analysis. The results show a clear departure between thermodynamic and kinematic boundary layer heights with the thermodynamic boundary layer height much shallower than the kinematic boundary layer height. Supporting the findings of Zhang et al. [31] based on the symmetric analysis of the dropsonde data, our results show that the hurricane boundary layer height increases with increasing radius in a storm-relative framework. This observed variation of boundary layer height with radius supports the theoretical scaling of a rotating boundary layer [46,59,61–64] in an axisymmetric framework. Our results indicate that this scaling also holds in a motion-relative asymmetric distribution of the boundary layer height.

**Figure 8.** Schematic diagram of the characteristic height scales of the hurricane boundary layer for the four quadrants relative to the storm motion. The height scales are based on the composite analysis of the dropsonde data. *h*inflow is the inflow layer depth (cyan dotted line); zi is the mixed layer depth (green dashed line); *h*vtmax is the height of the maximum tangential wind speed (blue solid line) and; hRic is the height of the bulk Richardson number value of 0.25.

In the eyewall region, all height scales show relatively small asymmetric distribution. The weak asymmetry is found in the kinematic boundary layer heights (*h*vtmax and *hi*nflow) that are slightly smaller in the left-rear quadrant than in other quadrants. All height scales demonstrate similar front-back asymmetry outside the eyewall region, in that the two front motion-relative quadrants have a deeper boundary layer. The front-back difference in the thermodynamic boundary layer height is smaller than

that in the kinematic boundary layer height. The mixed layer depth being smallest in the right-rear quadrant may be due to the SST cooling effect in this quadrant where the cold wake is located, following the argument of previous modeling studies [35,52].

The boundary layer height derived using the critical Richardson number (*Ric*) method shows a similar asymmetry as the kinematic boundary layer height, confirming that turbulence in the hurricane boundary layer is mainly shear driven. Consistent with Zhang et al. [31], our results suggest that the hurricane research community should use the kinematic height scales to represent the top of the boundary layer instead of the thermodynamic mixed layer depth typically used in non-hurricane conditions. Kepert [65] and Kepert et al. [66] discussed the limitation of using thermodynamic mixed layer depth to represent the boundary layer height in numerical models, which is in agreement with our observational findings here.

Our results in terms of the asymmetry of maximum tangential and radial wind speeds in the boundary layer generally agree with previous theoretical and numerical studies [46–48] Our result in terms of front-back asymmetry in the velocity fields also agrees with previous observational case studies such as the ABL structure in Hurricanes Mitch (1998), George (1998), Isabel (2003) and Daniel (2008) as shown by Kepert [37,38]. Although previous studies did not focus on the boundary layer height asymmetry, their results indicate a similar structure of the inflow and inflow layer depth as well as the height of the boundary layer jet. The surface wind asymmetry is consistent with Uhlhorn et al. [67] who analyzed stepped frequency microwave radiometer (SFMR) data, as well as Klotz and Jiang [68] who analyzed satellite data. The inflow asymmetry near the surface is consistent with the result of Zhang and Uhlhorn [69], who studied the characteristics of surface inflow angle in both axisymmetric and asymmetric framework.

Note that Zhang et al. [54] investigated the asymmetric structure of the hurricane boundary layer in relation to the environmental vertical wind shear. They focused on investigating how boundary layer thermodynamic structure is tied to the upper-level convection in an energy-cycling paradigm. They also studied the asymmetry of the boundary layer height scales relative to the environmental shear direction and found that the kinematic boundary layer height scales are larger in the downshear direction than in the upshear direction, while the thermodynamic boundary layer height scale is slightly larger left of the shear than right of the shear. The boundary layer being deeper in the downshear quadrants may be tied to the asymmetric distribution of convection that is usually initiated in the downshear-right quadrant and propagates to downshear-left quadrant [70]. The boundary layer is deeper in the quadrants where stronger convection occurs. Following this argument, storm-motion induced asymmetry of the boundary layer height scales may also be linked to asymmetric distribution of convective activity relative to storm motion, although storm motion and environmental wind shear are very different parameters.

The asymmetric structure above the boundary layer (i.e., convection) in tropical cyclones with respect to storm motion have been extensively studied using radar observations in both case studies and composite analysis studies. For instance, Jorgensen et al. [71] utilized flight-level data and found maximum upward mass transport in slow-moving storms occurred to the right of motion, with equal amounts occurring in the front and rear of the inner core. Our composite analysis supports strong convection occurs in the right-front quadrant, as both the maximum tangential wind speed and maximum inflow strength occur in this quadrant. Marks et al. [72] found that the maximum upward vertical velocities were in the left-front quadrant of Hurricane Norbert (1984). Marks [73] found that the maximum precipitation of Hurricane Allen (1980) in the eyewall region was in the left-front quadrant, while the maximum precipitation of Hurricane Elena (1985) was in the right-front quadrant. Through examining the patterns of reflectivity in the eyewall region of Hurricane Olivia (1994), Reasor et al. [74] found that the maximum radar reflectivity was located in the left quadrants relative to the motion direction, which was consistent with the structure in Hurricane Gloria (1985) documented by Franklin et al. [75]. Our composite analysis shows that the strongest inflow in the boundary layer is located in the front quadrants, which support large vertical motion and strong

convection due to mass continuity. Overall, it is hypothesized that the asymmetry of boundary layer inflow and jet strength are tied to the asymmetry of convection and precipitation.

Furthermore, the composite analysis result showing the front-back asymmetry of the boundary layer heights supports the theoretical argument that the boundary-layer convergence is larger in the front quadrants relative to the storm motion direction [45,47]. The surface wind asymmetry related to storm-motion effect induces the asymmetry of surface drag forcing, which in turn affects the distribution of boundary layer convergence and convection. Larger wind speed observed in the front quadrants induces larger turbulent mixing in the boundary layer, which supports larger kinematic boundary layer height according to dynamic scaling [42,45].

Of note, the combined effects of storm motion and environmental wind shear on the asymmetry of the boundary layer structure remains to be understood due to the limitation of data sampling size at the current stage. Future studies will combine the dropsonde and Doppler radar data to investigate the linkage between the boundary layer and convective processes and their asymmetry relative to both the storm motion and environmental wind shear. The extent of asymmetry of the boundary layer structure to the storm motion speed will also be evaluated when more observational data are available than those used in the present study.

**Author Contributions:** Conceptualization, J.A.Z., X.W.; methodology, Y.R., J.A.Z.; software, X.W.; validation, Y.R., J.A.Z., S.R.G. and X.W.; formal analysis, Y.R., X.W.; resources, J.A.Z., X.W.; data curation, Y.R.; writing—original draft preparation, Y.R., J.A.Z., S.R.G., X.W.; writing—review and editing, J.A.Z., S.R.G., X.W.; visualization, X.W.; supervision, J.A.Z.; project administration, Y.R.; funding acquisition, J.A.Z., S.R.G., X.W.

**Funding:** This study was partially supported by NSF Grant AGS1822128, NOAA Grant NA14NWS4680030, NASA Grant NNX14AM69G, and NASA Weather and Atmospheric Dynamics program (Grant NNH16ZDA001N-WEATHER) directed by Ramesh Kakar.

**Acknowledgments:** The authors would like to thank all the scientists and crew members who have been involved in the hurricane field program and operational reconnaissance hurricane missions to help collect the dropsonde data used in this study. This paper is a follow-up work of the first author's Master thesis and she would like to acknowledge Prof. Xiaolei Zou and Prof. Ming Cai for their helpful discussions. We also wish to thank the reviewers for their comments that led to improvement of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
