Kelvin-Helmholtz billows and their effects on mean state during gravity wave propagation

The Kelvin-Helmholtz (KH) billows which appear in the process of gravity wave (GW) propagation are simulated directly by using a compressible nonlinear twodimensional gravity wave model. The differences between our model and others include: the background field has no special initial configuration and there is no initial triggering mechanism needed in the mesosphere and lower thermosphere (MLT) region to excite the KH billows. However, the initial triggering mechanism is performed in the lower atmosphere through GW, which then propagate into the MLT region and form billows. The braid structures and overturning of KH billows, caused by nonlinear interactions between GWs and mean flow, can be resolved precisely by the model. These results support the findings in airglow studies that GWs propagating from below into the MLT region are important sources of KH billows. The onset of small scale waves and the wave energy transfer induce the shallower vertical wave number power spectral densities (PSD). However, most of the slopes are steeper than the expected k−3 z power law, which indicates that GWs with 10 km vertical wavelength are still a dominant mode. The results also show that the evolution of mean wind vary substantially between the different processes of GWs propagation. Before the KH billows evolve, the mean wind is accelerated greatly by GWs. By contrast, as the KH billows evolve and mix with mean flow, the mean wind and its peak value decrease.


Introduction
Numerical simulations and observational results indicate that Kelvin-Helmholtz (KH) billows are an important phenomenon in the atmosphere and have significant mixing effects on the mean atmospheric structures and species (e.g.Fritts et al., 2003;Hecht, 2004;Hecht et al., 2005).Up to now, the numerical simulations of KH billows have been performed using a background wind field exhibiting wind shear and excited by a regular or stochastic perturbation in the billows region.For example, the shear horizontal velocity and density background fields were used in the models developed by Patnaik et al. (1976) and Caulfield and Peltier (2000) (hereinafter designated PC).In addition, the models developed by Klaassen andPeltier (1985a, b, 1991) (hereinafter designated KP) and Hill et al. (1999) and Werne and Fritts (1999) (hereinafter designated HW) were all performed using the shear horizontal wind and temperature background fields, and the KH billows were excited by a small amplitude stochastic perturbation.Similarly, the shear horizontal field models were used by Fritts et al. (1996) and Palmer et al. (1996) (hereinafter designated FP), where KH billows were excited by a regular vertical velocity perturbation.The shear horizontal wind method was also employed by Fritts et al. (2003) to simulate KH instabilities.
As mentioned above, the initial background field has a special configuration, and an initial triggering mechanism to excite KH billows.In these situations, it is convenient to study the evolution of KH billows and their interactions with mean flow.However, the KH billows in the mesosphere and lower thermosphere (MLT) region may be caused by various triggering mechanism in various background fields.Therefore, there are three aspects of this problem which have to be addressed.The first question involves the triggering mechanism required to excite KH billows.The results from airglow and rocket studies have suggested that GWs play important roles in atmospheric dynamics and are capable of perturbing Published by Copernicus Publications on behalf of the European Geosciences Union.
the basic atmospheric state in the MLT region (e.g., Taylor and Hapgood, 1990;Li et al., 2005;Hecht, 2004;Hecht et al., 2005;Larsen et al., 2005;Yokoyama et al., 2005).The second problem relates to the shear background field.In fact, numerical simulations have shown that the nonlinear interactions between GWs and mean flow can also cause very large shear wind in the MLT regions even though the background atmosphere is initially at rest (e.g.Liu et al., 1999Liu et al., , 2006)).The third aspect deals with the interactions between the triggering mechanism and the shear background field.The simulation studies by Horinouchi (2004) showed that breaking of GWs generated by cumulus convection may create billows and induce conspicuous airglow ripples in the MLT region.Accordingly, GWs provide an important triggering mechanism to excite KH billows as they have the potential to induce strong wind shear.We speculate that GWs are capable of generating KH billows as they propagate in an initial windless background atmosphere.If this hypothesis is true, then non-linear interactions between GWs and initial zero mean flow are also an important physical mechanism to generate KH billows in the MLT region.Moreover, the characteristics of KH billows (such as the evolution, power spectrum and their effects on the mean wind) induced by this mechanism also need to be studied.
In this study we aim to test the feasibility of whether GWs formed in the lower atmosphere could cause KH billows in the MLT region, when GWs propagate into the MLT region in the background field without any initial configuration.A two-dimensional numerical model is used for the investigation (Liu et al., 2006(Liu et al., , 2008)).Although the two-dimensional model is unable to describe the secondary instability that occurred in a three-dimensional model (Fritts et al., 1996), it should be noted that, Fritts et al. (1996) have also showed that KH billows simulated by two-dimensional and spanwise averaged three-dimensional model are nearly identical before the maximum amplitude is reached.Therefore, it is reasonable to simulate the KH billows evolution before the maximum amplitude is reached by using a two-dimensional model.Consequently, we primarily focus our study on KH billows before its maximum amplitude has been reached.
Our paper is structured as follows.In Sect.2, we provide a brief description of the two-dimensional time-dependent non-linear GW model.The initial configurations of background atmosphere and triggering mechanism will also be described in this section.The evolution of KH billows is presented in Sect.3. Also, in this section, we focus our study on the evolution of KH billows that appear in the time interval after GWs become unstable and before they break.Sections 4 and 5 present the evolution of wave spectrum and mean wind as the KH billows evolve, respectively.The conclusions of this paper are given in Sect.6.

Numerical model and initial conditions
The model solves a two-dimensional, fully nonlinear, compressible, non-hydrostatic equation of motion set in a nonisothermal atmosphere.The solution's method and the stability of the model were described by Liu et al. (2006).The main result by Liu et al. (2006) was that energy density and dispersion relation of the model coincide well with linear theory when the model satisfies the conditions of small amplitude and isothermal background.For finite amplitude GWs, the evolution of potential temperature contours, simulated by our model, are similar to those obtained by Walterscheid and Schubert (1990) and Liu et al. (1999).We should note that there are two improvements to the model.First, the model used in this paper includes eddy and molecular diffusion in the momentum and thermodynamic equations, respectively.The diffusion coefficients are similar to those used in the study by Xu et al. (2003).Consequently, the model provides a more realistic atmospheric viscosity.Second, the increased spatial and temporal resolutions are further improvements of the model.As a result, the model is now capable of resolving fine structures as a result of non-linear and unstable GW dynamics.
The initial perturbation is a GW packet which is introduced at the height of 30 km.The horizontal velocity perturbation has the following form (Zhang and Yi, 2001;Xu et al., 2003;Liu et al., 2008), where the wave numbers in the horizontal and vertical directions are k x = 2π l x and k z = −2π l z , respectively.The horizontal and vertical wavelength are l x = 100 km and l z = 10 km, respectively.The amplitude of the initial horizontal velocity perturbation is ū = 4.0ms −1 , the scale height is denoted as H , which is defined as H =RT(z)/g and depends on the initial background temperature profile.The center of the wave energy packet is located at z 0 =30 km, providing a large enough vertical region for GWs to freely propagate upward and interact with the mean flow.The half width of the wave packet is λ z =15 km.The other initial perturbations, such as w (x,z,t = 0), p (x,z,t = 0), ρ (x,z,t = 0) and T (x,z,t = 0), are derived from u (x,z,t = 0) using the polarization equations for a linear GW.The initial background temperature profile is calculated from MSIS-90 (Hedin, 1991) at 45 • N and day number 82 (vernal equinox).The initial background wind in the horizontal and vertical direction is 0 ms −1 .An arithmetic sum of the initial perturbations and the initial background quantities is regarded as the initial state of the model.The computational domain extends from the ground to 130 km in the vertical direction, and the horizontal range is from 0 km to 100 km and equal to one horizontal wavelength.
The vertical and horizontal grid spacings are set to 0.1 km and 0.5 km, respectively.In accordance with the Courant-Friedrichs-Lewy (CFL) conditions for the explicit third-order total variation diminishing (TVD) type Runge-Kutta method, the time step is set to 0.25 s for integration in the horizontal direction.While the symmetrical time splitting method requires that one horizontal time step matches with two half vertical time steps to increase the resolution in the vertical direction, therefore, each of the half vertical time steps is 0.125 s (Ferziger and Peric, 1996).Consequently, we may anticipate that the fine structures of GWs can be well resolved by the model with relatively higher resolution.We also hypothesize that KH billows that appear in the process of non-linear unstable GW propagation will also be resolved by the model.

The evolution of KH billows
In the study of fluid or atmospheric dynamics, KH billows are usually described by the potential temperature or vorticity field (e.g.Patnaik et al., 1976;Klaassen and Peltier, 1985;Palmer et al., 1996;Fritts et al., 1996).The evolution of KH billows will also be described by the potential temperature in this paper for the following two reasons.First, the potential temperature is a conservative quantity following the motion of an air parcel, and second, its evolution is similar to that of the sodium mixing ratio (Xu et al., 2006).
Figure 1 gives the spatial distribution of the potential temperature field (contour), convectively unstable (N 2 <0, shadow), and dynamically unstable (0<Ri<0.25, the region between the shadow and the thick line) regions during GW propagation.Where the buoyancy frequency N and the Richardson number Ri are defined as, and C p is the specific heat at constant pressure.From Fig. 1, we can find that the dynamically unstable region is approximate to the convectively unstable region in two-dimensional simulations (e.g.Fritts, 1984).The potential temperature and the spatial distribution of the Richardson number indicate that a GW becomes unstable before the KH billows appear.
The initial instability of the GW occurred at 7300 s (which is regarded as the initial unstable time) and at the height of about 93 km.Subsequently, the overturning of the potential temperature increases greatly due to the non-linear interactions between the GW and the mean flow.
Figure 2 presents the evolution of the total potential temperature during the process of non-linear unstable GW propagation.The simulation results show that the potential temperature is in an overturning state at 8950 s, however, the overturning is consistent even before 8950 s.Since the potential temperature at the time interval between the initial unstable time of 7300 s and 8950 s (at this time new bending appears, denoted as "a" in Fig. 2) are similar to other simulation results (e.g.Walterscheid and Schubert, 1990;Liu et al., 1999Liu et al., , 2006;;Zhang et al., 2001;Xu et al., 2003), they will not be shown here.However, when time equals 8950 s, as shown in the frame drawn in Fig. 2, new bending appears below the tip of the overturn (denoted as "a" in Fig. 2).Since the aim of this paper is to study the evolution of KH billows that appear in the process of GW propagation, the potential temperature that describe the new bending (at 8950 s) and the subsequent evolution of KH billows are presented in Fig. 2.
It should be noted that, compared with the time scale of GW propagation, the time scale of the formation and decay of KH billows is very short and about several minutes (Hecht et al., 2005).In addition, the spatial scale of KH billow, which is less than 10 km in vertical direction (Hecht et al., 2005), is also smaller than the vertical domain where the GW propagates.Therefore, it is difficult to describe the evolution of KH billows in a direct numerical simulation if the model does not have sufficient spatial and temporal resolutions.Because of the improved spatial and temporal resolutions, our model can resolve the KH billows that appear before the GW breakdown into turbulence.It should also be noted that the initial background atmosphere and the initial triggering mechanism required to excite the KH billows have no initial configuration.In fact, the initial background atmosphere and the GW perturbation used in our model are similar to those used in Zhang et al. (2001) and Xu et al. (2003) and Liu et al. (2006), and the initial GW perturbation is in the lower atmosphere and far from the MLT region where the KH billows appear.Therefore our model is different from KH billow models with shear background field and external triggering factors which are located in the region where the billows will appear.Now, we analyze the evolution and fine structures of KH billows with the help of Fig. 2  isolines below the tips of the overturn (denoted as "a" in Fig. 2) also continue to overturn upward.Meanwhile, the upper parts of the overturning isolines (noted as "b" in Fig. 2) become gradually thinner (as shown in Fig. 2 at 9000 s).At 9050 s, the isolines below the tips of the overturn (noted as "c" in Fig. 2) continue to overturn upward as they had at 9000 s.At 9100 s, the upper part of the overturning isoline (noted as "d" in Fig. 2) is engulfed into the inner part of the other isolines (noted as "e" in Fig. 2).At 9150 s, the situations are similar to that of 9100 s.In addition, the tip of one isoline has been engulfed by the others, and the braids between two KH billows appear.
At 9200 s, KH billows and the braids can be seen more clearly than those at 9150 s (as shown in the frame drawn in Fig. 2 at 9200 s).Compared with the braids at 9150 s, the width of the braids becomes thinner.Meanwhile, the KH billows begin to overturn, and this is one of the main characteristics of the formation of KH billows.Subsequently, the KH billows continue to overturn, as shown in Fig. 2 at 9250 s and 9300 s.The fact that the width of the braids further decreases is the other characteristics of the evolution of KH billows.
After 9300 s, smaller scale structures appear in KH billows.Subsequently, these smaller scale structures and their evolutions can not be resolved precisely due to the limited resolution of the model.However, as for the potential temperature at 9350 s and 9400 s, the model can resolve them qualitatively.The large-scale structures seen from the potential temperature are the same as the evolution of KH billows.Afterward, the potential temperature is very chaotic and develops to turbulence after 9400 s.
In order to present the evolution of KH billows more clearly, in Fig. 3, we present the detailed structures of KH billows at 9250 s, 9300 s, 9350 s and 9400 s.The horizontal range is confined to 60 km-100 km, and the height range is confined to 94 km-108 km.At 9300 s, the braid structures can be seen more clearly (as shown in the frame drawn in Fig. 3 at this time).Subsequently, the braids have been ruptured due to the rotation of the KH billow core (at 9350 s).The rupture of the braids is also a main characteristic of KH billows evolution.Then, the potential temperature field is in turbulence, as seen in the contours at 9400 s.
The horizontal scale of KH billows at 9300 s is about 15 km (from x=73 km to about x=88 km, judged by the KH billows in the middle of the figure at this time).Similarly, the vertical scale of this billow is about 8 km (from 96 km to about 104 km).The core of this billow is at 100 km.We should note that there is no similar model to compare the KH billows with our model; however, we may qualitatively compare our results with those obtained by Hill et al. (1999).When the ratio between the streamwise size of KH billow L 0 and the shear depth of mean wind h equals 4π, the most rapidly growing mode can be excited efficiently.Judging from the mean wind profiles at 9200 s, 9250 s and 9300 s, we find that the shear depth is about 1 km (see the layer of positive wind shear abound 100 km altitude in Fig. 5c).Therefore, the streamwise size of a KH billow is about 13 km and qualitatively supports our simulation results.In addition, the Fig. 1b in Caulfield and Peltier (2000) indicates that the ratio between the horizontal and vertical scale of KH billows is about 2. This also supports our results.
In fact, the onset of turbulence starts from the instability when wave energy starts to cascade to smaller scales.From onset to well-developed it will take about 1 buoyancy period, which is about 300 s for the parameters used in our model.Our simulation result is well consistent with this (from 9100 s to 9400 s).Airglow images show that the timescale of KH billows from the beginning of formation to approaching maximum amplitude is about 5 min (Hecht et al., 2005).Our simulation results support this observational result.Unfortunately, the subsequent evolution of KH billows can not be resolved by our model with the current resolution.It should be noted that we have tested the stability of the model with respect to the representation of the KH billows and their temporal evolution by improving the resolution still further.The result indicates that the KH billows can evolve over a longer time after the GW overturning and the fine structures induced by breaking GWs can be resolved more clearly.Hecht et al. (2005) also used four models to compare their observations and concluded that the timescale depends on the initial conditions of the mean flow chosen for the simulations.The timescale obtained from the FP model is about 4.5 min and qualitatively consistent with our results.The comparisons between our results and other model results and observations at mesopause heights indicate that the KH billows that we obtained in our simulation are realistic.The analysis in this section indicates that, in the real atmosphere, the KH billows can be brought out and evolve for one buoyancy period.The KH billows appear after GW overturning and before the excitation of turbulence induced by breaking GWs.These results also support the findings in airglow studies that the GWs in the lower atmosphere are important sources in generating KH billows when they propagate to the MLT region (e.g., Taylor and Hapgood, 1990;Li et al., 2005;Hecht, 2004;Hecht et al., 2005).Although the spatial and temporal scale of KH billows evolutions are smaller than that of the whole process of GW propagation, the KH billows not only have significant triggering effects on turbulence, but also have significant mixing effects on the mean wind.The variation of mean wind caused by the mixing effect of KH billows will be discussed in detail in Sect. 5.

The evolution of GW spectrum
The shape of the GW spectrum is an important parameter in GW parameterization schemes and plays a significant role in the study of the global structure of the atmosphere in the MLT region.Therefore, in this section, we will discuss the evolution of the GW spectrum.
The vertical grid spacing used in this paper is 0.1 km.Since a wave can be described reasonably by 5 grid points per wavelength in numerical simulations, the corresponding minimum wavelength is 0.5 km and the maximum wave number is 4π×10 −3 (cyc/m).Because the vertical wavelength of the initial perturbation is 10 km, which corresponds to a vertical wave number of 2π×10 −4 (cyc/m), we have to calculate the spectral slopes in the wave number range as 2π×10 −4 -4π×10 −3 (cyc/m).
When we analyze the vertical wave number power spectral densities (PSD) during the stages of KH billow evolution, the vertical range is confined to 90 km-110 km, because the KH billows evolve in this height range.Therefore, the variations of wave spectrum caused by the evolution of KH billows may be described more precisely.However, as we can see from Fig. 2, the evolutions of the individual KH billows are not identical to each other at this height range.Therefore, the wave spectrum embodies the KH billow spectrum only qualitatively.In order to ensure that all individual spectral densities contribute equally to the shape of the mean PSD, a normalized individual PSD of each vertical profile was averaged as the mean PSD, similar to the ones used in Allen and Vincent (1995) and Zhang and Yi (2008).The mean vertical wave number PSD is presented in Fig. 4, Fig. 4a and Fig. 4b are the mean PSD of the horizontal velocity and temperature perturbations, respectively.The PSD and their slopes at different times are identified by different colors.
Our model can resolve the wave number range 2π×10 −4 -4π×10 −3 (cyc/m) precisely.However, the smaller scales with wave numbers larger than 4π×10 −3 (cyc/m) may contaminate the waves with wave numbers in the range 2π×10 −4 -4π×10 −3 (cyc/m) due to the non-linear wavewave interactions.In fact, a real turbulence spectrum would form only after the turbulence is well-developed.Therefore, the PSD and its slope at 9500 s only represents the turbulence spectrum qualitatively.
Although the PSD shown in Fig. 4 exhibits temporal evolution, we find that, the wave number corresponding to the peak of the PSD is 2π×10 −4 (cyc/m), and the corresponding wavelength is 10 km.This indicates that, although the small scale waves are excited during this process, the wave with 10 km vertical wavelength is still prominent.
The spectral slope of the horizontal velocity perturbations (herein after, when we mention the spectral slope, we refer to its magnitude) is increasing until 9300 s.Afterward, the slope begins to decrease.For the temperature perturbations, the slope is increasing until 9200 s and slightly earlier than that of the horizontal velocity perturbations.Compared with the KH billows shown in Fig. 2, the KH billows are welldeveloped before 9300 s.This indicates that the GW is unsaturated and its amplitude is increasing.By contrast, after the KH billows are well-developed, the small scale waves appear and their intensities increase due to the non-linear wavewave and wave-flow interactions.According to the dispersion relation of a linear GW l z = |c − u| N, where u is the mean wind and c is the phase speed of a GW (e.g., Fritts and Alexander, 2003), the vertical wavelength of a GW will decrease in the regions where the buoyancy frequency increases and/or the intrinsic phase speed decreases.Caused by the wave-flow interactions, the mean wind will be accelerated by the GW in the direction of the GW phase speed during the stage of non-linear GW propagation.It should be noted that we have also calculated the shear and curvature term in the dispersion relation derived by Nappo (2002, P29, Formula 2.34).The results showed that the shear term is an order of magnitude smaller than the curvature term, and the curvature term is an order of magnitude smaller than the buoyancy term.Therefore, the impacts of the shear and curvature terms on vertical wavelength are not very significant and can be neglected.Moreover, wave-wave interactions cause the wave energy transfer from larger scale waves to smaller scale waves.As the KH billows evolve, the small scale structures also appear in the billows, which intensify the small scale waves.As a result, the PSD increases at higher wave numbers and decreases at lower wave numbers, and the spectral slope decreases gradually after the KH billows are well-developed.
The time-averaged mean spectral slope for the horizontal velocity perturbations during the period of KH billow evolution is 3.51, which is larger than the value of 3.07 obtained by Zhang and Yi (2008).While the time-averaged mean slope for the temperature perturbations is 3.34, and it approximates to the value of 3.33 obtained by Zhang and Yi (2008).The results in Dewan and Good (1986) and Smith et al. (1987) indicated that, the saturated GW PSD obeys the k −3 z power law.However, the slopes simulated by our model are larger than 3, and the spectrum is dominated by the GW with 10 km vertical wavelength.Limited by the resolution used in the current model, the subsequent evolution of PSD can not be easily resolved.The PSD presented in Fig. 4 shows the non-linear wave-flow and wave-wave interaction, we tentatively associate most of these slopes with the unsaturated GW PSD.In addition, the spectral slope is consistent with the theoretical result obtained by Liu (2007) that the PSD slopes of horizontal velocity and temperature perturbations are in the range between −2 and −4.
The differences in spectral slopes may be attributed to the following two reasons.First, there are rarely purely monochromic GWs in the real atmosphere.However, in this paper we focus our study on the KH billows which are excited by a quasi-monochromic GW in the lower atmosphere.Therefore the slopes represent the spectral evolution of a quasi-monochromic GW.Second, the KH billows evolved only in a local horizontal range 50 km-100 km and not in the whole horizontal range.As a result, the mean PSD for the whole horizontal range 0 km-100 km may not present a well-developed KH billow spectrum.The slope at 9500 s is shallower than that at 9400 s, this indicates that the slopes may be even more shallow in the subsequent evolution of KH billows.However, this question can only be answered if we further improve the resolution of the model and extend the model to three dimensions.However, currently our model is limited by the resolution and can not resolve the KH billows after the maximum amplitude is reached.

The evolution of mean state
Similar to the previous simulation results, the mean wind is accelerated by the stable and unstable GWs before the KH billows appear (e.g., Liu et al., 1999Liu et al., , 2006)).In this section, we focus our study on the evolution of the mean wind during the process of KH billow evolution.In Fig. 5, we present the vertical profiles of mean wind at the initial unstable time and during the process of KH billow evolution.According to the evolution of KH billows shown in Fig. 2 and the mean winds shown in Fig. 5, the impact of a GW or KH billows on the mean winds has three stages, with different controlling mechanisms.The first stage happens from the time when a GW becomes unstable (7300 s) to the time when KH billows appear initially (9050 s) (as shown in Fig. 5a and  b).During this stage, the mean wind is accelerated to the peak value about 75 ms −1 (at the height of about 99 km).This stage may be regarded as the fast acceleration process due to the GW momentum deposition (e.g., Liu et al., 1999).
The second stage happens when KH billows appear (9050 s) until KH billows are well-developed (we use "appear" to describe the early stage of KH billows) (before 9200 s).In order to see the mean wind in the height range where the peak value appears, the height range in Fig. 5b and c is confined to 93 km-102 km.During this stage, the peak value at the height of 99 km increases slowly; and there is a new peak at the height of about 96 km (as shown in Fig. 5b).Comparing with the height range where the KH billows appear, we find that the evolution of mean wind during this stage may be attributed to the following two reasons.First, the core of a KH billow is at the height of about 100 km.Therefore, the KH billows mix the mean wind in this range.Second, the GW momentum deposition accelerates the mean wind in the region where the GW has not evolved into KH billows.The tradeoff between mixing and acceleration results in the slow acceleration of mean wind.
The new peak at the height about 96 km may be attributed to the following two reasons.First, the half width of the wave packet is 15 km.Despite the fact that the upper part of the wave packet evolved into KH billows, the mean wind below 96 km is still accelerated by the lower part of unstable GWs and results in the increase in peak value.Second, the edges of the KH billows occur at a height just below 96 km.Therefore, the mixing effect of KH billows on the mean wind is weak.The strong accelerative effect of GWs and the weak mixing of KH billows result in the increase of the peak values at the height of 96 km.
During the third stage (after 9200 s), mixing due to the KH billows becomes the dominant mechanism.As shown in Fig. 5c, induced by the mixing effects of KH billows on the mean wind, the peak value at the height of about 99 km decreases gradually.In addition, the peak value at the height of about 96 km increases before 9300 s.This is because this height is at the edge of the KH billows and the mixing effect on mean wind is weaker.However, as the height range of KH billows increases, the peak value decreases gradually at the height of about 96 km.This results in the decrease of mean wind in the height range of 94 km-102 km.Therefore, this stage may be regard as the deceleration stage of mean wind.
Compared with the previous simulation results, the mean wind evolution in the third stage is difficult to be resolved with the lower resolution model.This is because the fluid field in this stage is approaching turbulence and can not be resolved well by the lower resolution model.By using a model with a relatively higher resolution, we can simulate the nonlinear unstable GW evolution further.Consequently, we find that the mean wind evolution is different during the different stages of GW propagation.Our result indicates that the mean wind accelerates prior to the formation of KH billows and decreases afterwards.

Conclusions
The numerical model used in this paper is able to simulate the non-linear unstable GW propagation in an atmosphere with a realistic background temperature and diffusions.The increased resolution of the model permits us to simulate the KH billows which appear in the process of non-linear GW propagation.The acceleration effect of GWs and the deceleration effect of KH billows on the mean wind were also studied.Moreover, the evolutions of the PSD of the horizontal velocity and temperature perturbations were also analyzed.
The evolution of KH billows is simulated using a twodimensional model.The differences between our model and others include: the background atmosphere is initially at rest, and there is no external triggering mechanism required in the MLT region to excite the KH billows.In fact, the initial triggering factor (GW) is forced in the lower atmosphere and far from the MLT region where the KH billows will appear.Our simulation shows that the KH billows in the MLT region excited by GWs in the lower atmosphere appear in the time interval between the time of a GW overturning and it's breakdown into turbulence.During this time interval, the KH billows were excited directly by the non-linear interactions between GWs and mean flow.The braids and overturn of KH billows can be resolved precisely by the model.
The studies on the vertical wave number PSD of the horizontal velocity and temperature perturbations indicate that, before the KH billows are well-developed, the increasing amplitude of the unsaturated GW results in an increasing slope.However, after the KH billows are well-developed, the small scale waves and the energy transfer from larger scale waves to the smaller scale waves lead to a decreasing slope.
The mean wind was accelerated greatly by the stable and unstable GWs before the KH billows appear.However, the impact of the KH billow evolution on the mean wind experiences significant variations.During the initial stage of the KH billows evolution, the mean wind was accelerated slowly due to the GW momentum deposition.Afterward, the mixing of KH billows causes the mean wind and its peak value decreases.
It should be noted that, although the results in this paper are simulated with relatively high resolution, the twodimensional model can not precisely resolve the fine structure of KH billows and the secondary instability.In addition, the further evolutions of KH billows can not be resolved well, and thus a further improvement of the resolution is required.The model also needs to be extended to include threedimensions.

Fig. 1 .
Fig. 1.Contour plot of the total potential temperature and Richardson number when a GW becomes unstable.

Fig. 2 .
Fig. 2.Contour plots of the total potential temperature during GW propagation.These figures also present the evolution of KH billows.

Fig. 3 .
Fig. 3. Contour plots of the total potential temperature as the KH billows evolve.

Fig. 4 .
Fig. 4. Normalized PSD of the horizontal velocity (a) and temperature perturbations (b) versus vertical wave number.The PSD and their slopes at different times are shifted downward by a factor of 10 −1 and identified by different colors.

Fig. 5 .
Fig. 5. Vertical profiles of mean wind at several times during the process of KH billows evolution.