High‐resolution simulation of a severe case of low‐level windshear at the Hong Kong International Airport: Turbulence intensity and sensitivity to turbulence parameterization scheme

Large eddy simulations are performed for a severe windshear case of terrain‐induced airflow disturbances at the Hong Kong International Airport. The objectives of the study include: (i) to investigate the performance of large eddy simulation in capturing the turbulence intensity to be encountered by the aircraft; and (ii) to find out the impact of the choice of turbulence parameterization scheme on the simulation results. The structure functions of the headwind profiles are calculated, and the simulated velocity and the LIDAR data are found to have similar behavior. The cube root of eddy dissipation rate is also calculated, which is taken to be the metric of turbulence intensity in aviation application. All EDR calculated has similar behavior. They have good correlation, but the best‐fit‐straight‐line of the scatter plot is found to have a slope far away from unity. The reason for this needs to be studied further. Finally, a simulation with the use of Smagorinsky scheme is performed. It is found to have better performance in capturing pilot windshear reports in terms of the relative operating characteristic curve. The use of Smagorinsky scheme alone may be an alternative in the choice of turbulence parameterization scheme in the alerting of windshear using large eddy dissipation, but the computational cost is too expensive. It is important to note that the present results are based on one case only and more cases would need to be studied in order to establish the robustness of the conclusions. More experiments would also need to be conducted, for example, in which nest(s) to use the Smagorinsky scheme.


| INTRODUCTION
The Hong Kong International Airport (HKIA) is situated in an area of complex terrain. To the south of the airport is the mountainous Lantau Island with peaks reaching 1000 m above mean sea level and valleys of 300 m height in between. As the dominant wind direction is from the east through southeast to southwest in Hong Kong, the airflow over the mountainous island may result in airflow disturbances in the airport area that could be hazardous to the arriving/departing aircraft. Detection of windshear has been conducted by using a number of equipment. However, prediction of windshear remains to be a challenging problem.
The simulation of terrain-induced windshear at HKIA has been demonstrated using 200 m horizontal resolution of Weather Research and Forecasting (WRF) model (Chan & Hon, 2016). It is able to capture the major feature of the terrain effect on the airflow, but still of insufficient resolution to resolve the tiny vortices/waves in the turbulent airflow. In Chan et al. (2021), large eddy simulation has been performed using Regional Atmospheric Modeling System (RAMS) with a spatial resolution of 40 m. It is found to successfully capture the tiny vortices and waves in a particular severe windshear case. However, the turbulence intensity of the airflow has not been considered.
In Chan et al. (2021), the Deardorff turbulence parameterization scheme (Deardorff, 1980) was chosen to perform large eddy simulation (LES) of terrain-disrupted airflow at the Hong Kong International Airport for a severe case of low level windshear using RAMS version 6.3. The model has been demonstrated to have skills in capturing the pilot windshear reports. There remain a number of questions with the study, namely: i. Statistical distribution of the simulated velocity has been shown in the previous study. But the more conventional metric of the turbulent structure based on the simulated velocity and the actual Doppler Light Detection and Ranging (LIDAR) velocity, such as the structure function, is not shown. The present study would compare with the structure function calculated from the actual LIDAR velocity and the simulated LIDAR velocity; ii. The common metric for turbulence intensity for aviation application is the cube root of eddy dissipation rate (EDR), EDR 1/3 , which is the dissipation rate of the turbulent kinetic energy, and this has not been discussed in the previous study. The present study would calculate the EDR from the actual LIDAR velocity, the simulated LIDAR velocity, and the energy dissipation term from the turbulent kinetic energy governing equation, and compare all three of them; iii. The previous study also has not addressed the question of the sensitivity of the modeling results to the choice of turbulence parameterization scheme of RAMS 6.3.
The present study attempts so address the above questions.
The case under study is a terriain-induced windshear event at the Hong Kong International Airport on March 5, 2015. This is the case with the largest number of pilot windshear reports for a nonthunderstorm and a nontropical cyclone day. The windshear event had been simulated before using a numerical weather prediction model with a spatial resolution of 200 m. However, probably due to limitation of the spatial resolution, this model manages to reproduce the return flow downstream of the terrain but is not capable of simulating the small-scale (in the order of several hundred meters) vortices in the terrain-disrupted airflow. This topic has been discussed in detail in Chan et al. (2021). The present article aims at studying this event using large eddy simulation method, with emphasis on the turbulence intensity aspect.
There are various ways of simulating small-scale atmospheric flow, including direct numerical simulation (DNS, which is the direct implementation of fluctuated values into the Navier-Stokes equation without any turbulence model), large eddy simulation (LES, which is an average turbulence model in which filtered Navier-Stokes equations are used for large-scale eddies with an appropriate model for solving small-scale eddies), and Reynolds averaged numerical simulation (RANS, which is a mathematical model based on average values of variables for both steady-state and dynamic flows). In the present study, LES with sub-grid scale parameterization based on a couple of schemes available in RAMS 6.3 would be tried out. Moreover, structure function of the actual and simulated LIDAR radial velocity would be calculated. The structure function is an important ingredient for calculating turbulent kinetic energy dissipation rate, which is the turbulence intensity internationally adopted for alerting clearair turbulence and other turbulence flow for the aircraft.
The state of the art of the meteorological development of large eddy simulation has been summarized, for example, in Lilly (2000). The application of Deardorff (1980) scheme is still considered to produce reasonable result. This point has been confirmed further, for example, in the article by Gibbs and Fedorovich (2016). The Smagorinsky (1963) scheme available in RAMS 6.3 would also be considered. Though with a rather long history, these schemes are still considered to give good results in the mesoscale to microscale numerical weather prediction models available nowadays.
LES study for meteorological applications have been considered in a number of studies in the literature. Ito et al. (2020) and Shimoyama et al. (2013) have considered applications of LES for a couple of airports in Japan. The former tried to use LES to explain an incident of the aircraft at the airport. Though sophisticated meteorological data (e.g., Doppler LIDAR) have been used for comparison with the simulation results, the simulated airflow has not been compared with the "sky truth" of windshear or turbulence, namely, the automatically generated windshear alerts based on LIDAR data and, most importantly, pilot windshear reports. Moreover, the terrain under consideration in these two papers is not particularly complex and rugged as compared with the complex terrain at the Hong F I G U R E 1 (a) Shows a sample of the headwind profile and the meaning of the various quantities in Equation (1). (b) Shows the structure function from headwind profiles, from the LIDAR and model simulation Kong International Airport. Complex terrain with wind energy application has been considered in the LES study of, for example, Uchida (2018), but only time series data at a specific point have been used for verification of LES results. Krishnamurthy et al. (2010) has considered the LES simulation of EDR for a complex terrain for T-REX experiment in the US. However, only a vertical profile of EDR has been considered for verification, instead of the spatial distribution of EDR (e.g., in the form of an EDR map in plan-position indicator scans of the Doppler LIDAR) over a larger area. Verification of LES has also been conducted for a long period of time at a relatively flat area in the Netherlands (Schalkwijk et al., 2015).
Compared with all the previous studies, the present article aims at making a small step forward in the aviation application of LES with the following new features: (i) comparison with an EDR map from the Doppler LIDAR at an extended period of time, (ii) comparison with pilot windshear reports in order to verify the LES results with actual observations, and (iii) comparison of two sub-grid scale turbulence parameterization schemes on an actual operational environment.
The statistics of the wind field in terms of the structure function are considered in Section 2. This sets the starting point for the calculation of EDR. The comparison of EDR is shown in Section 3. Finally, the choice of different parameterization schemes is compared in terms of the structure function and EDR as in Section 4. The various sections are related to the calculation of turbulence intensity and the choice of the turbulence parameterization scheme and are included in a coherent manner.

| STRUCTURE FUNCTION
A re-cap to structure function is given here. Spatial statistics of the wind field b v R, θ,k ð Þ, as measured in the form of the radial velocity of the Doppler LIDAR at radial distance R, azimuthal angle θ and time step k, are assumed to be described by the structure function b D R 1 , R 2 ð Þ. The structure function measures correlation between radial velocities at positions R 1 and R 2 along the direction of the headwind. The headwind is taken to be the direction along the flight path/glide path/runway. The structure function in this article is based only on the headwind profile measured by the LIDAR, that is, the radial velocity along the flight path/glide path/runway direction, with the headwind being the wind in the runway orientation to be encountered by the aircraft for lift. The structure function is not meant for windshear alerting/ detection, but just shows the statistics of the wind field. The structure function is calculated up to a distance of 3 nautical miles (about 6 km) from the runway end, as this is the region for alerting windshear of the aircraft.
The headwind profile is constructed by considering the headwind component along a particular glide path of a runway corridor. For each headwind profile, the fluctuating velocity component b v 0 at every step in space and time is first obtained by removal of the spatial mean v, which can be calculated by fitting a straight line to the profile of the radial velocity along the flight path/runway: An example of the headwind profile, the spatial mean and the fluctuating component is shown in Figure 1a.
The structure function at a particular separation S can then be written as: The summation over time (index k) is performed over a 1-h period and the summation over R is to go through all the data points along the flight path/glide path. N is the total number of valid data parts at position R in the 1-h period. For a given separation S, those pairs of positions with the separation S are considered. An example of the meaning of R and S is shown in Figure 1a. The structure function is expressed as a function of the separation. The same calculation is performed for both the LIDAR headwind profile and model-simulated headwind profile for the most commonly used runway corridor, namely, 07LA, that is, landing at the north runway of the Hong Kong International Airport (HKIA) from the west. Some results of the structure function are shown in Figure 1b, namely, the hours 03-04 UTC, 04 to 05 UTC, and 07 to 08 UTC. It could be seen that the LIDAR-based structure function and model-simulated structure function are very similar. As a result, apart from the statistical distributions of velocity fluctuations in Chan et al. (2021), it is again shown here the turbulent velocity components of the actual measurements from the LIDAR and the simulation results of RAMS have similar statistical behavior. This adds confidence to the use of computer simulation to derive the air turbulence quantity, such as the cube root of eddy dissipation rates, which is the same as the turbulent kinetic energy (TKE) dissipation rate.
The structure function is calculated along the glide path of the aircraft only, namely, up to 3 nautical miles (about 6 km). With larger separation, the amount of paired data points for calculation of the structure function would decrease and the so-calculated structure function may appear to be fluctuating. It is anticipated that only data up to a separation of 2 km would be reliable, showing the general increase of the structure function with separation. Normally only data points in this region are used to fit with the theoretical curve in the calculation of eddy dissipation rate. This point should be borne in mind in the interpretation of the results in Figure 1.
It is noted that structure function is a measure of the fluctuation of the wind speed. This may be related to windshear and turbulence to certain extent. However, according to the international civil aviation, windshear refers to headwind change and turbulence intensity is based on the cube root of the eddy dissipation rate, which would be discussed in the next section and is calculated based on the structure function, instead of using the structure function itself.

| EDR 1 / 3
There are three sources of EDR 1/3 , namely: i. Based on the actual LIDAR measurements, following the approach of Chan (2010)-the LIDAR radial velocity data are used to calculate the average over a sector and the deviation of the individual measurement from the average is used to calculate the structure function and the calculated structure function is fitted with the theoretical curve to determine the EDR; ii. For the RAMS-simulated velocity, the LIDAR-based radial component is extracted, and then input into the algorithm of Chan (2010), that is, treating the F I G U R E 2 EDR maps from actual LIDAR data (a), model-calculated LIDAR velocity (b) and model direct output (c). The comparison among the three diagrams show that the part encircled in the figures has the best comparison (between the actual [a] and the simulated results in [b] and [c]). The model tends to over-estimate the EDR at the region to the west of this circle RAMS-simulated radial velocity as the "true" LIDAR-based radial velocity, and the following the method as described in (i) above to calculate the EDR; and iii. Based on the TKE dissipation term of Deardorff scheme of RAMS, in a way similar to Chan (2009)the direct model output from the TKE dissipation rate term in the TKE equation is used and regarded as EDR.
As in Chan et al. (2021), the model in the article is initialized at 00 UTC. It would take about 1 h for the model to be fully developed. The study of EDR, structure function, windshear, etc. is considered after the first hour of running of the model.

F I G U R E 3
Scatter plot (density plot) of the EDR calculated from the three methods. The equation of the best fit straight line to the data point and the corresponding correlation coefficient are shown at the top of each plot. The correlation coefficient is described as r^2, and this "r" does not mean the radial of the LIDAR. The fitting passes the significance test of 5% An example of the EDR map for 3-degree plan position indicator (PPI) scan of the LIDAR for the three method is shown in Figure 2. It could be seen that all of them are showing higher turbulence for the east to southeasterly flow closer to the mountains. They appear to have quite similar spatial distributions. The various panels of Figure 2 are separated by a number of seconds only, namely, a separation of 6 s between Figure 2a,b, and a separation of 1 s between Figure 2b,c.
In order to compare the EDR values obtained by the various methods more quantitatively, scatter plots in the form of density plots has been prepared for the whole F I G U R E 4 The model-simulated 3 PPI LIDAR velocity imagery for the use of Smagorinsky scheme in all five nests 14-h simulation period of RAMS for this severe windshear case. Figures 3a,b show the scatter plot of the model simulated EDR, namely, direct EDR output in Figure 3a and the calculated EDR based on the modeloutput radial velocity in Figure 3b, against the actual LIDAR measurements. They show reasonable degree of correlation, though the slope is not unity and the yintercept is 0.1, instead of 0. Based on the correlation coefficient, the former has higher degree of correlation (r 2 = 0.185) than the latter (r 2 = 0.156). Note that this r refers to the correlation coefficient, and should not be confused with the radial R in Equation (1) for the radial direction of the LIDAR.
The two model based EDR calculation methods are compared in Figure 3c. The correlation coefficient is even higher (0.377). However, it is not quite natural to note F I G U R E 5 (a) Shows a sample of the simulated surface wind field with the use of Smagorinsky scheme in all five nests. (b) Is the corresponding surface wind field with the use of Smagorinsky scheme in the first two nests, and Deardorff scheme in the next three nests that the slope is far from unity. Probably this is just a single case to verify and the results and the comparison includes data for 14 h. With longer simulation time, the simulated wind pattern becomes more different from the actual LIDAR observations. More cases need to be studied in order to see the correlation of the EDRs obtained from the two methods.

| ANOTHER TURBULENCE PARAMETERIZATION SCHEME
In the previous simulation, Smagorinsky Scheme Smagorinsky (1963) is used in the first two nests, and Deardorff scheme is used in the next three nests. To study the sensitivity of the modeling results to the choice of the turbulence parameterization scheme, a simulation is performed with Smagorinsky scheme being used in all five nests. The simulation is very slow and takes about 7 to 8 times of the physical time required for the previous simulation study. Again, the severe windshear case is studied for a 14-h period from 00 to 14 UTC, March 5, 2015. Deardorff (1980) scheme is applied to both vertical and horizontal mixing, so that the turbulence so simulated is isotropic and the diffusion coefficients are the same in all directions. The prognostic TKE equation is solved. The dissipation term in the TKE equation, viz. the EDR (ε), is given by: where l is a subgrid-scale mixing length which depends on the atmospheric stability (see Deardorff (1980) for details), E the TKE and C D = 0.19 + 0.51 l/(Δx Δy Δz) 1/3 (Δx is the grid size in the x-direction, etc.). The Smagorinsky mixing-length model has been designed to simulate the energy transfer from resolved to unresolved scales across an inertial subrange of locally isotropic three-dimensional turbulence. The resolved motions are separated from residual motions by implicit filtering of the governing equations in space. The residual stress-tensor is defined by a linear eddy viscosity model. An implementation of this scheme has been described in F I G U R E 6 Similar to Figure 1b but with the use of Smagorinsky scheme for five nests detail in https://www.ocf.berkeley.edu/$langhans/papers/ Langhans_2012_COSMO_Smag.pdf. Figures 4a,b show the simulation results (in the form of the simulated LIDAR PPI scan velocity imagery at 3 elevation angle) at two times, namely, around 04 UTC and 08 UTC. Instantaneous shots are given, and they could be compared with similar figures in Chan et al. (2021). The general airflow patterns are similar, but the timing and size of the small-scale vortices as discussed in Chan et al. (2021) could be different in this simulation. In particular, at 04 UTC, the present simulation does not seem to give the small scale vortex. There is only towards-the-LIDAR flow emerging from the mountain gap to the southwest of the airport (an elongated area of green). This emerging flow is quite persistent and also shows up at the latter time of 08 UTC.
To better visualize the difference between the simulations with the use of Smagorinsky scheme in all five nests and with the use of Smagorinsky scheme in the first two nests, and Deardorff scheme in the next three nests, the surface wind fields from the two simulations are shown in Figure 5. It could be seen that the former is more "turbulent", and there are many small-scale wave/vortex features in the airflow. On the other hand, the latter is relatively smooth and appears to be less turbulent. Given the limited number of surface wind observations in the airport region, it is not certain which simulation result is more "realistic" in comparison with the actual observations.
To study the statistical behavior of the simulated wind, the structure function of the headwind profile is calculated again. The results are shown in Figure 6, for the times 02 to 03 UTC, 05 to 06 UTC, and 07 to 08 UTC. The structure functions of the simulated flow and the actual LIDAR velocities are very similar. Again this shows that the simulated airflow has similar statistical behavior as the actual observations.
To see which choice of parameterization scheme would perform better, a relative operating characteristics (ROC) curve based on 07LA pilot windshear reports is calculated for the two simulations. The results are shown in Figure 7. It could be seen that the present case (Smagorinsky scheme F I G U R E 7 ROC curve for the two sets of simulations. The alerting thresholds are shown near the data points at all nest) is found to perform better than the previous calculation, with an optimum alerting threshold of around 12 knots. Again this is based on a single simulation only.
More results would be necessary to determine which model setup would perform better in capturing the pilot windshear reports.
F I G U R E 8 (a) Shows a sample of the headwind profile with the use of Smagorinsky scheme in all five nests for a case with pilot windshear report. The windshear ramps are highlighted in red. (b) Is the corresponding headwind profile with the use of Smagorinsky scheme in the first two nests, and Deardorff scheme in the next three nests In order to study the performance of the windshear capturing, sample heqdwind profiles from the two simulations, namely, the use of Smagorinsky scheme in all five nests and the use of Smagorinsky scheme in the first two nests, and Deardorff scheme in the next three nests are compared for a case of pilot windshear report in Figure 8. As shown previously in the surface wind field, the former simulation tends to give more "turbulent" wind field, and thus, for a given windshear ramp threshold, say, 12 knots, the former is able to capture the wind field with a fixed percentage of time of windshear alerts. On the other hand, the latter simulation tends to be less turbulent more "smooth" in the airflow, and thus the windshear ramp has a magnitude of 11 knots only, which is not able to capture the pilot windshear report with a threshold of 12 knots given a fixed percentage of time of windshear alerts. This example shows the difference in the two simulation models.

| CONCLUSION
The present study is a study of a severe windshear case to study the turbulence of the model results. For either choice of turbulence parameterization scheme, the structure function based on the headwind profile is similar to that from the actual LIDAR observations. This adds weights to the credibility of the simulated airflow in terms of the statistical behavior of the wind velocity.
The use of Deardorff scheme permits the consideration of EDR calculation. Three sets of EDR data are obtained, and they are found to have reasonable correlation. However, the slopes of the scatter plots are far from unity, and it appears to be rather surprising. More simulation results would be required to address this point.
Also ROC curves are calculated for the two sets of simulations. The use of Smagorinsky scheme all the way appears to have better performance in capturing the pilot windshear reports. Again, this result is obtained on a single simulation only and is still subject to larger-scale study. Even though the area under ROC curve for Deardorff (1980) scheme is smaller than that based on Smagorinsky scheme, it does not undermine the application value of the former scheme in the successful simulation of the terrain-disrupted airflow in the present case, for example, in Chan et al. (2021) and this article. More windshear cases based on these two schemes would be considered in the future studies, and we would see if the comparison results in the present article would be applicable to a larger sample of windshear cases.
Moreover, the dynamical core of RAMS 6.3 and that of Weather Forecasting and Research (WRF) model in the operational 200 m simulation may be considered to have similar merits and it would be hard to draw conclusion based on a single case that RAMS 6.3 is superior in terms of dynamics, thermodynamics and model physics. Again, a larger sample of terrain-induced windshear cases would be considered in the future study to see if there is significant difference in the model behavior, particularly when both are used for LES modeling with a horizontal spatial scale in the order of several tens of meters.
Further experiments would also be conducted for the various combination of the turbulence parameterization schemes in the nested runs, for example, with the use of Smagorinsky scheme in the first three and four model nestings. This would test the sensitivity of the modeling results to the particular choice of parameterization schemes.
The authors are now performing calculation for a month long period in order to address some questions raised in the above discussions. The results would be reported later. ORCID Pak Wai Chan https://orcid.org/0000-0003-2289-0609