Numerical Simulation of a Lee Wave Case over Three-Dimensional Mountainous Terrain under Strong Wind Condition

This study of a lee wave event over three-dimensional (3D) mountainous terrain in Lantau Island, Hong Kong, using a simulation combining mesoscale model and computational fluid dynamics (CFD) model has shown that (1) 3D steep mountainous terrain can trigger small scale lee waves under strong wind condition, and the horizontal extent of the wave structure is in a dimension of few kilometers and corresponds to the dimension of the horizontal cross-section of the mountain; (2) the life cycle of the lee wave is short, and the wave structures will continuously form roughly in the same location, then gradually move downstream, and dissipate over time; (3) the lee wave triggered by the mountainous terrain in this case can be categorized into “nonsymmetric vortex shedding” or “turbulent wake,” as defined before based on water tank experiments; (4) the magnitude of the wave is related to strength of wind shear.This study also shows that a simulation combining mesoscale model and CFD can capture complex wave structure in the boundary layer over realistic 3D steep terrain, and have a potential value for operational jobs on air traffic warning, wind energy utilization, and atmospheric environmental assessment.


Introduction
Lee wave is one of the key indicators of how mountains affect air flows [1].It also has significant impact on aviation safety, wind energy utilization, and air pollution.Over the years, scholars have studied lee waves through different approaches, including theoretical analysis [2][3][4][5][6][7], field observations [8][9][10][11], and numerical simulation [12][13][14].In the light of the previous studies, the formation of lee waves is closely related to atmospheric stratification and buoyancy force [15].When mountain-crossing winds occur in stable stratification of the atmosphere, leeward flows will form waves and even lee wave rotors if stable stratification is strong.
Though there were some studies discussing or summarizing the lee wave phenomena related to three-dimensional (3D) terrain [12,13,15,16], most of the previous studies mainly focused on two-dimensional (2D) or semi-2D ridge.
Compared with the 2D ridge, the 3D terrain is far more complex and the research conclusions for 2D issues may not be valid for 3D issues.Actually, the study on the lee wave phenomena related to 3D terrain is insufficient, and the characteristics and the mechanism of the lee wave triggered by 3D terrain are not fully understood till now.Thus, it is necessary to do more in-depth researches on the lee waves related to 3D terrain to improve the understanding in this field.
From analysis of Doppler radar data, Chan [17] noticed significant lee waves in Lantau Island, Hong Kong, under a typical strong wind condition, and he took it as "vortex/wave shedding." Lantau Island is a typical 3D steep mountainous area, as elevation change reaches more than 700 m within a 2 km range.Compared with simple and idealized mountains [12] or semi-2D ridges [14] discussed in previous numerical studies, the terrain of Lantau Island is far more complex and realistic, which makes the lee wave case observed in this area highly representative and worth further in-depth studies, especially these related to characteristics of lee waves such as spatial structure, impact range, and life cycle related to 3D terrain.
Based on the authors' previous studies, this paper leverages the advantageous application of computational fluid dynamics (CFD) in microscale numerical simulation to perform a detailed simulation on the "vortex/wave shedding" case in Lantau Island raised by Chan [17] and to analyze the process of lee wave formation triggered by 3D mountainous terrain under strong wind condition.

Brief Descriptions of Simulation Method
The simulation period for this study is chosen to be 18:30-19:30 local time on September 29, 2011.During this period, Typhoon Nesat was transiting Hong Kong, steady southeast wind covered the entirety of Lantau Island, and surface wind speed was keeping steady at above 10 m/s.
The simulation uses a mesoscale model, namely, Regional Atmospheric Modeling System (RAMS) [18], which is offline coupled with FLUENT, a CFD software [19].Its flowchart is shown in Figure 1.The simulation first uses 20 km resolution data exported from the Hong Kong Observatory's Operational Regional Spectral Model (ORSM) to initiate a four nested grids RAMS simulation.It then exports 800 m resolution RAMS results as boundary conditions to initiate FLUENT calculations.This procedure is similar to those used in the authors' previous studies [20][21][22] and has proved to be reliable to illustrate flow field characteristics under complex underlying surfaces.
RAMS simulation has the same model setting as that of Chan [17,23].There are totally four grids for RAMS simulation, and grid 1 has a horizontal resolution of 4 km with 77 × 67 grids.The first vertical layer has a height of 40 m with the stretching ratio of 1.15.Grid 2 has a horizontal resolution of 800 m.The number of points and vertical levels is the same as grid 1. Mellor-Yamada scheme is used in grid 1 and Deardorff [24] scheme is used in grid 2.More details on the RAMS simulation can be found in [17,23].For comparison purpose, RAMS simulations in grid 3, namely, a horizontal resolution of 200 m, would be used.Deardorff [24] turbulence parameterization scheme is adopted in this grid for "large eddy simulation" mode.Cumulus parameterization is switched off.The first grid has a height of 20 m and a stretching ratio of 1.15 is maintained.
In FLUENT simulation, the atmosphere is assumed as incompressible, and the main governing equations are as follows.
Momentum equation: Continuity equation:  where   is the average velocity of air flow,  is the pressure,       is the Reynolds stress,  is the density of atmosphere,  is the molecular viscosity, and   is the body force.
The FLUENT simulation domain has a size of 22 km × 19 km, including Lantau Island, Hong Kong, and its surrounding areas.A computer aided design (CAD) model is built using terrain elevation data, as shown in Figures 2(a) and 2(b).Simulation domain is discretized by hexahedral meshing.The cell numbers for  and  directions are 140 and 120, respectively, which means that horizontal resolution is about 158 m.In vertical direction, simulation domain reached about 3000 m elevation and is divided into 22 layers.The size of each vertical grid cell sketches in a ratio 1 : 1.1 from the one below, which means that the vertical spaces of the bottom layer cells at sea level are about 45 m.
The FLUENT simulation is performed in two steps.The first step is a steady-state simulation in which time derivative terms and variables in dynamic equations are removed.In the first step of simulation, turbulent closure model uses the realizable k-e scheme in the Reynolds Averaged Navier-Stokes (RANS) framework.Once numerical solution reached convergence in steady state, we initiate the second simulation stage, a non-steady-state simulation.In the second stage simulation, variables depend on time and LES turbulent solution with Smagrinsky-Lilly scheme for sub-grid option are used in the simulation.This two-step simulation procedure guarantees that LES simulation yields higher quality and more physically meaningful initial value.
Considering there are steady southeast winds throughout the entire simulation period, the east, the north, and the south boundaries are set to velocity inlet type.Wind speed data at the three boundaries are taken from the second grid of the RAMS simulation output at September 29, 2011, 18:00, and are input into FLUENT computational model via the boundary profile (BP) module.For the west boundary, boundary conditions are set to outflow, which means that spatial derivatives at the boundary are zero.The top boundary is also set as the velocity inlet boundary, and the wind speed and direction are directly set as the RAMS output values.For all velocity inlet boundaries, the wind speed and direction are set as unchanged during the simulation.The boundary conditions setting used in this study can provide strong forcing for the whole simulation domain and will help to ensure the dominant flow direction in the domain in accordance with the observed one.The bottom of the simulation region is set as no-slip boundary condition, with the roughness heights of the sea surface, and the airport in Lantau Island and the Lantau mountains are set as 0.0005 m, 0.005 m, and 0.5 m, respectively.The domain setting and the boundary conditions setting used in this study are similar to those used in the authors' previous studies [22] and have proven to be reliable to capture the characteristics of the flow field over mountainous area.

Radial Velocity.
In order to help learn the background atmospheric condition of the vortex/wave event discussed in this study, the radiosonde sounding data in Hong Kong at 12 UTC, September 29, 2011, are shown in Figure 2(c).
From Figure 2(c), it can be seen that during the transiting of Typhoon Nesat, the wind direction generally keeps as southeast and only slightly varies with height.The wind speed generally increases with height within the atmospheric boundary layer, and the average wind speed of the layer within the height of 2000 m is around 20.0 m/s.The temperature decreases with altitude, whilst there are ground inversion and two elevated inversion layers at approximately 650 m and 1200 m above sea level.Detailed report on observation of lee wave during Typhoon Nesat can be referenced to Chan [17].The major observations in Chan [17] are the shedding of wave/vortex from the mountains of Lantau Island from the radial velocity data from the weather radar.The radar observations are reproduced in Figure 3, with the shedded waves/vortices enclosed in ellipses.The warm colour (red, orange, yellow, and brown) means that is the radial velocity that away from the radar, and the cool colour (green, blue, and purple) means that the radial velocity is heading towards the radar.As a result, the general wind direction is east-southeasterly.The green features highlighted are embedded in an area with generally away-from-the-radar flow and are considered to be waves/vortices.Figure 3 showing 6 instances of Dropper radar near Lantau Island, with Figure 3(a) shows the observation at 18:56 local time and around 1-minute interval between each subsequent figure.As shown in Figure 3, there is a wave/vortex structure distinctly different with surroundings (indicated by green patch in red circle in figure) formed in mountain, labeled as "A, " and then it shedded off and gradually moved downstream until dissipation.The entire procedure of this event lasts about 6 minutes.
Starting from an initial field obtained from a steady RANS simulation in the first step, the FLUENT LES simulation needs some time to adjust to get an equilibrium field.For this case, the adjusting time is about 25 minutes, and from the 27th minute of the simulation, the vortex/wave shedding event begins to appear.Thus, the simulation results starting from the 27th minute are used to analyze the characteristics of the simulated lee wave process.The radial velocity maps drawn by FLUENT simulation results are illustrated in Figure 4.It can be seen that FLUENT successfully captures the lee waves observed by radar.There is a wave structure different with surrounding wind direction generated in northwest Lantau Island and continued to move downstream, and the entire process also lasts about 6 minutes.Figure 4 also shows that the wave structure extends about 3 km downwind from mountains, which is a typical microscale lee wave phenomenon.
Comparing Figures 3 and 4, FLUENT can accurately reproduce the lee wave observations, and the simulated evolution and impact areas of the wave structure are consistent with the observed ones.
Apart from the wave/vortex "A, " which shedded off from Nei Lak Shan, there are also waves/vortices "C" and "B" that shedded off from Cheung Shan and Tai Tung Shan, respectively, as shown in Figure 3. Their sheddings are also reproduced successfully by FLUENT (see Figure 4).

Horizontal Wave Structure.
Figure 5 depicts the simulated flow field at 220 m above sea level in northwest Lantau Island, and this elevation is about where radar observed lee waves.Figure 5 clearly shows that the green patch in Figures 3  and 4 is a wave structure, and FLUENT shows the process in which the wave structure detached from the mountain and moved downstream.Horizontal scale of each wave structure is about 3 km (in the east-west orientation) and about 1.5 km (in the north-south orientation) and corresponds to the horizontal cross-section of the mountain where it formed at 220 m elevation.The entire process of lee wave moving downstream lasts about 6 minutes for a distance about 3-4 km.

Vertical Wind Speed and Wave Structure.
To analyze wave structure in vertical direction, we chose a cross section in simulation domain to analyze the changes in vertical wind speed and the wave structure.This cross section runs in southeast-northwest direction roughly consistent with dominant background wind direction and is through three control points within simulation domain, p1 ( = 0 m,  = 15800 m,  = 0 m), p2 ( = 0 m,  = 15800 m,  = 3000 m), and p3 ( = 16900,  = 0 m,  = 3000 m).
Figure 6 shows the vertical wind speed component () in the above-mentioned cross section and can find the following characteristics: (1)  alternates between positive and negative values along the background flow direction on the vertical cross section, indicating that vertical wind speed is clearly fluctuating in a typical wave-like fashion; (2) from the distance between negative  value regions (blue area), it can be determined that the wavelength of the lee wave is about 3 km; (3) the movement of negative  value regions over time shows that wave structures continue to move downstream and dissipate.Deduced from the changes from Figure 6(a) to 6(e), wave moves about 2.5 km downstream within 4 minutes; (4) Figure 6 also shows that wave not only moves downstream but also regenerates in original location and the whole process continuously repeats.
For better visualization of the shedded wave, Figure 6 is replotted by showing the vectors of the flow field as projected on the cross sectional plane, as shown in Figure 7.The shedded wave is highlighted as a black arrow.It could be readily seen the shedding of the wave from the mountain from this wind field plot.
3.4.Pathline.Sections 3.2 and 3.3 describe the structure of shedded wave/vortex by using horizontal and vertical cross sections across the structure.In order to investigate the 3D structure of the wave/vortex, the "pathlines" of the flow (namely, following the trajectories of the air particles) are also analyzed, though the figure of the pathline is not provided in this paper.
It could be seen that, for structures "A" and "C" in Figure 3, the flows appear as oblique waves in the 3D space; that is, the wave axis does not lie purely horizontally or vertically, but appears as an oblique line in the 3D space.As a result, when horizontal or vertical cross sections are cut through this structure, a wave shows up in the cross sectional planes.
On the other hand, the structure downstream of Tai Tung Shan, namely, feature "B" in Figure 3, appears to be more complicated, at least as given in the model simulation.There appears to be a rotor downstream of the mountain.The rotor has an axis generally parallel to the orientation of the mountain ranges of Lantau Island, that is, roughly northeast to southwest orientation.areas of higher vorticity values, as coloured in green/yellow in Figure 8.The patterns are rather similar to the results of idealized model simulations of vortex shedding in, for instance, Schär and Smith [12].The vorticity plot provides another view of the vortex/wave shedding from the three major mountains of Lantau Island.

Vorticity. As an alternative view of the flow field in the model simulation, vorticity plot is given in
From the simulation results including pathlines and vorticity, it could be seen that vortex/wave shedding occurs downstream of isolated mountains of Lautau Island.Each wave/vortex is associated with each mountain.Compared to the whole simulation domain, each mountain occupies a relatively small area only.

Comparison with RAMS Simulation
Since there is RAMS simulation at horizontal resolution of 200 m, which is comparable with that of FLUENT simulation, there is a question on the comparison of the results from these two model simulations.In particular, there is a question on whether the vortex/wave shedding features as observed in actual data could show up more clearly in FLUENT than in RAMS.
The RAMS simulation results at a height of 220 m are shown in Figure 9.As in Figure 4, the simulated wind fields are resolved along the measurement radials of the weather radar in order to reproduce the radar-observed features for ease of comparison.As shown in Figure 9 (highlighted in a black ellipse), vortex/wave shedding is also observed downstream of Nei Lak Shan in RAMS simulations at 200 m resolution.However, compared to the FLUENT simulation, the magnitude of the wave in RAMS simulation is smaller.In general, the flow field downstream of the mountains of Lantau Island appears to be smoother.No significant wave/vortex shedding could be reproduced downstream of Cheung Shan and Tai Tung Shan in RAMS.Though the 50 m resolution RAMS simulation can reproduce the wave shedding process, it is too time-consuming and the wave structure is not as clear as that in FLUENT with a much coarser resolution.In terms of the computational time, the use of grid 4 (50 m resolution) would require around 10 times more about the time required for running the RAMS simulation up to grid 3 (200 m resolution) only.The use of FLUENT is computationally more efficient and it would be sufficient to simulate the vortex/ wave shedding with a resolution of about 200 m by properly resolving the complex terrain near the airport.From the comparison result, it appears that FLUENT is better in capturing the vortex/wave shedding features.On the other hand, the waves do not show up equally clearly in RAMS simulation.It shows that FLUENT has the potential of forecasting the occurrence of vortex/wave shedding in real time applications, for example, in providing an earlier alert about the occurrence of low-level windshear/turbulence arising from terrain-disrupted airflow over an airport.

Sensitivity of Wave Shedding to Wind Shear
In order to investigate the effect of vertical wind shear on the occurrence of wave/vortex shedding downstream of the mountains, two more numerical simulations are carried out with the modification of the vertical wind speed profile.The profiles in use are shown in Figure 10.They include the original profile (control, CTL), no wind shear (NS), and a stronger shear below 500 m above sea level (SS).For NS case, the average wind speed of CTL over the range of height under consideration is used.
The pathlines of the model simulation results for NS and SS shows that there are still waves/vortices and even rotors downstream of the mountains (figures not shown in this paper).As such, the generation and shedding of the waves/ vortices and the occurrence of rotors do not appear to be dependent on the vertical wind shear.It is believed that they are more related to the steep slopes associated with the mountains of Lantau Island under consideration.
However, the various simulations do show some differences.Figure 11 illustrates the vorticity plot on a horizontal plane.As could be expected, the SS case shows larger vorticity at a couple of kilometres downstream of the mountain.This is followed by CTL (maximum vorticity of about 0.8 s −1 ), and NS shows the smallest vorticity (vorticity generally in the order of 0.6 to 0.8 s −1 after a couple of kilometres downstream of the mountains, and the area of 0.8 s −1 vorticity is less extensive).Moreover, from the vertical cross section of the vertical velocity (figures not shown), SS gives the largest value of vertical velocity (in the order of 4 m/s or above), followed by CTL and NS (in the order of 2 to 3 m/s for the maximum value of the vertical velocity in the vertical plane).A larger vertical shear in the background flow is found to increase the vorticity and the upward/downward flow, which may be more hazardous to the aircraft.

Analysis of the Dynamics of the Vortex/Wave Shedding
In the present study, both observational data and numerical simulation suggest that lee waves can be formed in the lee of three-dimensional mountains under strong wind condition.So, a new task for the present study is to explain the dynamics of the phenomenon.A series of water tank experiments were performed by Lin et al. [25] to investigate the characteristics of the pattern for flows past a three-dimensional sphere.According to Lin et al. [25], the flow past sphere can be categorized into eight regimes, namely, steady two-dimensional attached vortices, unsteady two-dimensional attached vortices, two-dimensional vortex shedding, lee-wave instability, nonaxisymmetric attached vortex, symmetric vortex shedding, nonsymmetric vortex shedding, and turbulent wake.The eight regimes were located on an Fi-Re regime diagram, where Fi is the Froude number, and Re is the Reynolds number.
The study of Lin et al. [25] provides very valuable hints for explaining the dynamics of the lee wave/vortex shedding found in the present study.For the present case, based on the radiosonde data upstream of Lantau Island, the Brunt-Vaisala frequency, , is in the order of 0.01 rad/s and  0 is in the order of 20 m/s.As a result, the Froude number (ℎ/ 0 ) is about 0.376, where ℎ is the height of Nei Lak Shan (751 m) and  0 is the average background wind speed.According to Lin et al., it is in the regime of "nonsymmetric vortex shedding" or "turbulent wake" triggered by 3D terrain.
Vortex shedding in tropical cyclone situation at Hong Kong International Airport has been discussed in Shun et al. [26].Shun et al. [26] indicated that the vortex shedding may be related to "nonsymmetric vortex shedding" or "turbulent wake" with the Froude number in the order of 0.36, which is in accordance with the results in the present study.Nevertheless, the new result in this paper is that the vortex/wave shedding event is found to be successfully reproducible by numerical simulation, which provides many more data (such as vertical cross section and vorticity field) to help understand the vortex/wave shedding event.

Conclusions
This study of a typical lee wave case in Lantau Island, Hong Kong, using RAMS/FLUENT combined simulation concludes the following.
(1) As a typical 3D steep mountainous terrain, Lantau Island can trigger small scale lee waves under strong wind condition.Both radar observation and CFD numerical simulation captured clear evidences of lee waves, and some waves are in a formation of rotor.
(2) The lee wave triggered by the mountainous terrain can be categorized into "nonsymmetric vortex shedding" or "turbulent wake, " as defined before based on water tank experiments.The spatial dimension of the lee wave corresponds to the horizontal cross-section dimension of Lantau mountain, with a wavelength of about 3 km.
(3) In the present study, the life cycle of each lee wave is about 6 minutes.Wave structures will continuously form in roughly the same location, then gradually move downstream, and dissipate over time.
(4) The steep terrain is the major reason leading to the lee waves, though the magnitude of the wave is related to strength of wind shear.
Besides the above conclusions, this study once again shows that simulation combining mesoscale model and CFD can capture complex flow movement in boundary layer.CFD software, especially FLUENT, uses body-fitted mesh structure for domain discretization and uses finite volume method (FVM) in numerical calculation, which ensures CFD's capability to perform stable numerical calculations with nonconformal mesh.Therefore, CFD can obtain stable numerical solutions without the terrain being smoothed and ensures the accuracy of terrain.In addition, some special modules in CFD, such as user defined function (UDF, which has a grammar similar to C language in FLUENT) and boundary profile (BP), could help users code specific programs to describe boundary conditions and define some physical properties of fluid in simulations and thus provide an interface between the CFD and mesoscale models or observational data.The above special features make CFD have advantages over mesoscale models when describing wind field characteristics over steep terrain.The method of combining mesoscale model and CFD can ensure the accuracy of simulation by obtaining largescale circulation data in mesoscale simulation and detailed microscale topographical description in CFD simulation simultaneously.
Finally, it can be expected that the simulation method of combining mesoscale model and CFD would be applied in operational air traffic warning, wind energy utilization, and atmospheric environmental assessment in a near future.

Figure 2 :
Figure 2: Terrain within FLUENT simulation is shown in (a) and (b).(c) gives the background atmospheric condition for the simulation as obtained from the radiosonde sounding in Hong Kong at 12 UTC, September 29, 2011.The figures include the wind direction profile, wind speed profile, and temperature and humidity profiles.

Figure 4 :
Figure 4: FLUENT simulated radial wind speed at 220 m elevation ((a) represents the 27th minute of the simulation period with 1-minute intervals between each figure).

Figure 5 :Figure 6 :
Figure 5: FLUENT simulated wind field in m/s at 220 m elevation ((a) represents the 27th minute of the simulation period with 1-minute intervals between each figure.The location of the vertical cross section in Figures 6 and 7 is shown as a red dotted line in (a)).

Figure 7 :
Figure7: The wind vector projected on the vertical cross sectional planes.The wave shedding is highlighted in a black arrow (unit: m/s, time for each figure is same as that in Figure5).

Figure 8 :
Figure 8: Vorticity plot for the model simulation result at 31 minutes (magnitude of vorticity is in s −1 ).

Figure 9 :Figure 10 :
Figure 9: RAMS model simulation at a height of 220 m, with the flow field resolved in the direction of the measurement radials of the weather radar.A vortex/wave shedding feature is highlighted in a black ellipse.

Figure 11 :
Figure 11: Vorticity field at a height of 220 m above sea level at 31 minutes for NS (a) and SS (b).