Frequency variations of gravity waves interacting with a time-varying tide

Using a nonlinear, 2-D time-dependent numerical model, we simulate the propagation of gravity waves (GWs) in a time-varying tide. Our simulations show that when a GW packet propagates in a time-varying tidal-wind environment, not only its intrinsic frequency but also its ground-based frequency would change significantly. The tidal horizontalwind acceleration dominates the GW frequency variation. Positive (negative) accelerations induce frequency increases (decreases) with time. More interestingly, tidal-wind acceleration near the critical layers always causes the GW frequency to increase, which may partially explain the observations that high-frequency GW components are more dominant in the middle and upper atmosphere than in the lower atmosphere. The combination of the increased ground-based frequency of propagating GWs in a time-varying tidal-wind field and the transient nature of the critical layer induced by a time-varying tidal zonal wind creates favorable conditions for GWs to penetrate their originally expected critical layers. Consequently, GWs have an impact on the background atmosphere at much higher altitudes than expected, which indicates that the dynamical effects of tidal–GW interactions are more complicated than usually taken into account by GW parameterizations in global models.


Introduction
Atmospheric solar tides and gravity waves (GWs) are ubiquitous in the middle and upper atmosphere, while they also play a significant role in atmospheric dynamics Xu et al., 2009;Li et al., 2010;Huang et al., 2012). Tides are believed to strongly modulate the propagation conditions experienced by upward-propagating GWs (Beldon and Mitchell, 2010). In the context of energy conservation, tidal and GW amplitudes grow with increasing altitude because of the exponentially decreasing atmospheric density with height. They can eventually become large enough (several tens of m s −1 in terms of horizontal-wind amplitudes) to induce mutual nonlinear interactions in the mesosphere and the lower thermosphere. Therefore, interactions between GWs and tides have become an issue of great importance in middle-and upper-atmosphere studies.
Limited by the spatiotemporal resolution and coverage of different observation techniques, further understanding of the tidal-GW interaction mechanism and the associated alterations of the background and the tidal-and gravity-wave fields depends sensitively on sustained development of numerical studies. Most previous numerical studies of tidal-GW interactions can be sorted into two categories. The first explores these interactions on a global scale (Mayr et al., 1998(Mayr et al., , 2001Meyer, 1999;Manson et al., 2002) and focuses on the effects of the interactions on the variability of tidal and background structures rather than on the GW parameters. In this type of study, GW effects are realized using GW drag parameterization models because of the sub-grid scale of GW perturbations, which cannot be resolved explicitly in global models. The other type deals with GW phenomena on local scales (Liu et al., 2008) and focuses mainly on the effects of the shear-wind field induced by tides on GW propagation. In most such studies, the tide is assumed to be timeinvariant. Liu et al. (2008) numerically studied the propagation of GWs in a time-invariant tidal wind and found that the tidal wind-induced critical layer could almost completely absorb GWs propagating from below and prevent them from reaching higher altitudes. Nevertheless, it is noteworthy that some other studies (Eckermann and Marks, 1996;Senf and Achatz, 2011) have emphasized the importance of the tidal wind's time variability in the interactions. Eckermann and Marks (1996) suggested that the tidal wind's time variability (i.e., its acceleration) could increase the intrinsic frequency of the interacting GW. Senf and Achatz (2011) showed that critical layers disappear when the background wind is a timedependent wave.
For the purpose of investigating tidal-GW interactions, we simulated the propagation of GW in a time-varying tide based on a nonlinear, 2-D time-dependent numerical model. The model is based on the Navier-Stokes equations and extends from the ground up to the lower thermosphere. By introducing a time-variable diurnal tide as the background environment in which GWs propagate, we concentrate on the impact of time-variable diurnal tidal winds on GW parameters and propagation, and in particular on the GWs' groundbased frequency and energy propagation near the critical layer. The remainder of this paper is organized as follows. In Sect. 2, our numerical models, including the relevant equations governing the dynamics, the background state, specification of the tide, and the initial GW perturbation, as well as the numerical scheme, are introduced in detail. Section 3 presents the nonlinear simulation results, including the evolution of GW packets in the presence of a tidal background, the temporal and altitudinal variations of the GWs' groundbased frequency, as well as the heights of the critical layers. Finally, the conclusions are provided in Sect. 4.

Governing equations
GWs in the mesosphere and lower thermosphere are usually characterized by medium to high frequencies and medium to small scales. Their propagation can thus be regarded as a local dynamic process. Hence, we start from a set of hydrodynamic equations in a 2-D atmosphere that is compressible and isothermal; rotation effects are ignored. In our model, the sum of the time-independent zonally averaged fields and the time-variable tidal fields, including the wind, temperature, and density fields, constitutes the entire background environment experienced by upward-propagating GWs. The time-independent fields are assumed to be purely altitudedependent. In addition, to investigate the effects of tidal winds on GW propagation and isolate the possible effects of the time-independent background winds, we set the timeindependent background winds to zero. The tidal fields are time-and altitude-dependent, while their horizontal variation is disregarded in our simulations, because the horizontalpropagation domain of the GWs is less than 4000 km. GW fields are treated as perturbations on such a background field. Therefore, the governing equations can be written as where x and z are the horizontal and vertical (positive corresponding to the upward-increasing direction) coordinates, respectively; u and w are the horizontal and vertical perturbation velocities, respectively; T and ρ are the perturbation temperature and density, respectively; and T 0 and ρ 0 are the time-independent temperature and density, respectively. For simplicity, the time-independent background is assumed to be a horizontally homogeneous, isothermal atmosphere in hydrostatic equilibrium (T 0 = 288 K). R = 286.9821 J kg −1 K −1 is the gas constant; c v = 718 J kg −1 K −1 is the heat capacity at constant volume; c p = c v + R is the heat capacity at constant pressure; and γ = c p /c v .

Specification of initial tidal disturbance
u T , w T , ρ T , and T T in Eq. (1) are the zonal and vertical wind, density, and temperature components of the tide, respectively, which are of the form f a (z, t) cos[n t + f p (z, t)], where = 7.292 × 10 −5 rad s −1 is Earth's rotation rate, n = 1 represents the diurnal component, and f a (z, t) and f p (z, t) are the tidal amplitude and phase, respectively. Since strong tides likely induce strong tidal-GW interactions, we take the diurnal tidal results from the Global Scale Wave Model 2000 (GSWM-00; Hagan et al., 2001) at 30 • N in October as our initial tidal field. This diurnal tide is more significant than those at other latitudes and in other months. Since GSWM-00 only provides tidal results at altitudes of 0-124 km, the tidal-wind and temperature amplitudes are artificially attenuated to zero using an exponentially declining function above 124 km. If tides are regarded as large-scale perturbations of the basic state, tidal forcing should be included in the model. However, in our simulations we focus only on tidal-GW interactions and eliminate interference from the nonlinear evolution of the diurnal tide itself. Therefore, we regard the timevariable tide as part of the atmospheric basic state by artificially introducing residual terms (i.e., R u , R w , R ρ , and R T ), which keep Eq. (1) balanced in the absence of perturbations of this basic state; i.e., they ensure that the instantaneous tidal values can evolve freely with time. The temporal evolution of the tide's zonal-wind profile is illustrated in Fig. 1.

Specification of initial GW perturbation
An upgoing Gaussian GW packet is specified as the initial perturbation. Its horizontal perturbation velocity has the following form: where u c is the initial wave amplitude, which is set to 1 m s −1 , x c = 480 km and z c = 48 km are the initial geometric (x, z) center coordinates of the wave packet, and k x and k z are the horizontal and vertical components of the wavenumber vector, respectively. According to observations of GWs at mesospheric altitudes (Lue et al., 2013), the value of the horizontal wavelength, λ x , is set to 120 km, and an initial vertical wavelength λ z = 6 km, is adopted. The dominant GW ground-based frequency is then 1.02×10 −3 rad s −1 at z c and the full width at half maximum of the power spectrum for the ground-based frequency (not shown) is calculated to be 0.23 × 10 −3 rad s −1 . σ x = 120 km and σ z = 6 km are the half widths of the wave packet in the horizontal and vertical directions, respectively. The other initial perturbation quantities (i.e., w (x, z, 0), ρ (x, z, 0), and T (x, z, 0)) are derived from the polarization equations of GWs.

Computational method
We applied different methods to calculate the spatial differentiation in the two directions, including the Fourier collocation method in the horizontal direction and the finitedifference method in the vertical direction. A detailed description of the computational method was presented by Huang et al. (2006Huang et al. ( , 2007. Considering the spatial scales of the initial GW, x and z are set to 10 and 0.4 km, respectively. In this simulation, the horizontal and vertical calculation ranges are 0 ≤ x ≤ 4000 km and 0 ≤ z ≤ 200 km, respectively. Because the Fourier collocation method is applied in the horizontal direction, the lateral boundaries are periodic. As for the upper and lower boundaries, projected characteristic line boundaries are employed. A detailed specification of the latter can be found in Zhang and Yi (1999). Numerical experiments show that these boundary conditions work well.

Wave propagation
To observe GW propagation in a tidal environment, we show the temporal evolution of the horizontal perturbation velocity in Fig. 2. The horizontal perturbation velocity is composed of the GW's horizontal-wind perturbation, the altered background and horizontal tidal winds caused by the interaction among the GWs, tide, and background flow. Since we focus on the effects of tidal-GW interactions on GW propagation before gravity wave breaking, we only present 15 h of propagation process. Generally speaking, the GW packet propagates upward and eastward from the initial center position (x = 480 km, z = 48 km) and significantly changes its packet shape.
The initial GW packet covers several tens of kilometers in height. It can be divided into three parts with respect to the positive or negative vertical shear of the zonal tidal wind; i.e., S T = ∂u T ∂z , which is shown as the dotted curve in Fig. 1. The top and middle parts are located in regions of stronger shear; the former is in the positive-shear region while the latter is in the negative-shear region. It is well-known that, based on ray-tracing theory for internal GWs (Jones, 1969), the vertical wavelength/propagation velocity of upward-propagating internal GWs would become shorter/lower in the presence of a horizontal background wind with positive vertical shear and longer/higher in one with negative vertical shear. Therefore, at 1 h the vertical wavelength of the top part decreases while that of the middle part increases, leading to a phase tilt in the middle part that is steeper than that in the top part, and the entire packet seems to shrink. At 3 h, the top part has propagated beyond the height where peak A (see Fig. 1) is located and propagates into the negative-shear region. Its vertical propagation is accelerated until it enters the higher positive-shear region. Because of the alteration of the GW's vertical propagation velocity in different shear regions, the packet seems to distort significantly and turns into two separate segments (see the result at 3 h in Fig. 2).
After 3 h, the top segment of the GW packet seems to encounter the first critical layer near peak B and does not propagate upward for several hours, while the bottom segment continues going upward. Subsequently, the two separate segments merge and become a single entity again at 6 h. According to linear GW theory, the critical layer height is the altitude where the local horizontal phase velocity of the GW equals the horizontal background wind. Based on this criterion, we know that between 3 and 9 h, the top part of the GW packet encounters the critical layer at an altitude of approximately 91 km.
In previous simulations by Liu et al. (2008) of nonlinear interactions between GWs and the diurnal tide, the tidalwind and temperature components were time-invariant. If the tidal zonal wind were to remain unchanged with time (refer to Fig. 1, time = 0 h), the GW packet would encounter the critical layer below peak B (where the corresponding maximum of the tidal horizontal wind is 54.7 m s −1 , 3 times larger than the initial central phase speed of the GW packet, 17.3 m s −1 ) and could not propagate upward any longer. That is, the packet could not penetrate peak B. However, in our simulation, part of the packet eventually penetrates peak B and reaches higher regions. This is different from the results of the simulations of Liu et al. (2008). Our simulation results imply that the time variability of the tidal field should be taken into account when exploring the effects of the complex nonlinear interactions between GWs and tides on GW propagation.
After 11 h, part of the packet penetrates the first critical layer near peak B, rapidly passes through the negative-shear region, and encounters the second critical layer near peak C (see Fig. 1). Two packet centers form at that stage: one below peak C and another below peak B. After 13 h, part of the packet even penetrates the second critical layer near peak C and reaches heights of more than 120 km.
To summarize, while propagating in a time-varying tidal field, GWs would encounter several critical layers at different altitudes. After encountering the critical layers, the amplitude of the GWs would decrease with time because of absorption by the critical layer. The time variability of the tidal zonal wind is beneficial for the GWs to penetrate the critical layers and propagate upward to higher altitudes.
A similar conclusion was reached by Walterscheid (1981), who simulated the interactions of GWs and semidiurnal tides and found that although waves are strongly absorbed near the critical level, considerable transmission of wave energy occurs as well. Our work also supports previous results (Zhong et al., 1996;Broutman et al., 1997;Eckermann, 1997;Walterscheid, 2000;Sartelet, 2003a, b) in that the finite lifetime of some critical levels allows waves to survive the critical layer and reach higher altitudes, or waves can escape critical level removal when propagating in a time-dependent background.

Ground-based frequency of the GWs
From inspection of Fig. 2 we find that when a GW packet propagates in a tidal-wind field, it seems to be eventually divided into three segments: one is located below the first critical layer, the second between the first and second critical layers, and the third ends up above the second critical layer. We now explore the temporal and height variability of the ground-based frequency of the GW packets around the centers of these three segments. Taking the bottom segment as an example, we calculated the ground-based frequency of the GW, as follows. First, we locate the center position of the GW horizontal perturbation velocity (specified by the velocity's maximum amplitude) below the first critical layer from 3 to 12 h at each integer hour. We obtain 9 center positions for these 9 integer hours; i.e., (x k km, z k km), with k = 3, 4, . . . , 12. To clearly demonstrate the temporal evolution and height variation of the GW frequency around the segment center, for the integer hour k, we record a 6 h temporal series of the GW's horizontal perturbation velocity at temporal intervals of 15 min at each height node from (x k km, z k −6 km) to (x k km, z k + 6 km); i.e.,from 6 km below the center height up to 6 km above the center height. The height range is only a display scope, which does not assume a constant vertical wavelength at each integer hour. The 6 h time window starts at 3 h before and ends at 3 h after this integer hour; i.e., from k − 3 to k + 3 h. Next, we perform Lomb-Scargle periodogram analysis (Scargle, 1982) on each time series. The periodogram's maximum is regarded as the ground-based frequency of the GW at a certain hour and height. Figure 3 shows the height and time variability of the calculated ground-based GW frequency around the segment center below the first critical layer. We find that the groundbased frequency varies significantly with height and time. At 3 h, the ground-based frequency increases more significantly with height above the center than at 0 h (i.e., around the initial value of 1.02 × 10 −3 rad s −1 ), while it remains almost unchanged below the center. From 5 to 6 h, the groundbased frequency decreases with time above 66 km (the center at 6 h), and then becomes almost constant over 60-72 km (the entire altitude range shown at 6 h). From 7 h (when the segment center suddenly moves up to a height of 84 km) to 11 h, the ground-based frequency clearly increases with height over the full height range shown. The calculated maximum frequency can reach 2.04×10 −3 rad s −1 , which is twice the initial frequency. These maximum frequencies occur at the greatest heights shown in the panels in Fig. 3, near the transient first critical layer of the GWs determined by the time-variable tidal-wind field. This obvious increase of the GW ground-based frequency indicates that these GWs are potential to penetrate the critical layer, which should not happen for these GWs if they propagate in an invariant wind field. As for the temporal evolution of the ground-based frequency, it increases (decreases) clearly with time approximately above (below) the segment center. The resulting frequency is higher (lower) than the initial ground-based frequency. We also note that the frequency at the lowest height shown is very low. For example, the frequency at 11 h at a height of 75.2 km (0.65×10 −3 rad s −1 ) is notably lower than the initial frequency.
Based on these results, we conclude that the groundbased frequency of a GW packet propagating in a timevariable background wind can vary significantly with both height and time. This result cannot be explained on the basis of linear GW theory in a time-invariant backgroundfield frame, which predicts that the ground-based frequency of GWs does not change. Zhang et al. (2000) and Zhang and Yi (2002) have numerically shown that nonlinearity and height-dependent molecular viscosity decrease the GW frequency, but they have not found a similar frequency increase as in the present study. More interestingly, the timevariable frequency increases at some heights but decreases at other heights. For the case explored here, the frequencies above and below the segment center show different temporal variations. A relevant investigation by Broutman and Young (1986) stated that the time dependence of the background flow removes the constraint that the absolute frequency of a small-scale wave is constant along the ray, and the time dependence of the oscillating inertial current also eliminates critical layers, which would otherwise be important for short-wave dynamics. Walterscheid (2000) also found that the low-frequency approximation overpredicts variations in the intrinsic frequency in a time-dependent background wave field because the variation of the observed frequency is exactly out of phase with the variation of the Doppler term. However, it is the time dependence of the intrinsic frequency rather than that of the observed frequency that was presented in his work. Senf and Achatz (2011) pointed out that when the time dependence of thermal tides is included in the description of GW propagation, GW frequencies are modulated in time. This leads to penetration of their transient critical level, which was also found in our simulation. Some interesting questions arise from our simulation results: what causes the variation of the GW frequency in the absence of molecular viscosity? And which factors dominate the increase or decrease of the frequency?
On the basis of ray-tracing techniques, Jones (1969) deduced the temporal derivative of the ground-based frequency for internal GWs, Here, u and w are the background horizontal and vertical wind components, respectively. Clearly, Eq.
(2) indicates that the temporal variations of the background winds change the GW frequency. In the case of our simulation, since the initial background flow is set to zero, the background wind of the upward-propagating GWs may consist of two parts. These include the tidal field due to the low frequency with respect to the frequency of the initial GW (period < 2 h) and the background winds arising from the nonlinear interactions between the background and the waves (including tidal and gravity waves; hereafter "nonlinear winds"); i.e., u N and w N , respectively. The nonlinear background winds are regarded as the wind perturbation components with spatial scales larger than the GW spatial scales. Because the GW vertical wavelength changes greatly during the propagation while the horizontal one almost remains unchanged, we calculated the nonlinear background winds at k h by averaging the horizontal and vertical perturbation velocities in the horizontal direction at x k . The average scale should be larger than the GW dominant horizontal wavelength, but a toolarge value leads to an underestimation. So we choose it to be 4 horizontal wavelengths. We found that the horizontal and vertical background wind speeds are no greater than 2.5 and 0.002 m s −1 , respectively. This is at least one order of magnitude smaller than the amplitudes of the tidal horizontal and vertical winds, respectively, at the same height. Moreover, we also calculated the nonlinear horizontal-wind acceleration ( ∂u N ∂t ) and found that its absolute value was smaller than 1.7 × 10 −4 m s −2 . This is much lower than the absolute value of the tidal horizontal-wind acceleration ( ∂u T ∂t ). The same situation applies to the nonlinear vertical-wind acceleration. Thus, in the following analysis, the effects from nonlinear winds can be ignored, so that we have u ≈ u T and w ≈ w T . In addition, we compared k x ∂u T ∂t and k z ∂w T ∂t , and found that the former term is at least 14 times larger than the latter. Therefore, Eq. (2) can be rewritten as Equation (3) predicts that in the present case, the tidal horizontal-wind acceleration ∂u T ∂t dominates the GW frequency variation. Positive (negative) accelerations cause the frequency to increase (decrease) with time.
To quantitatively investigate the impact of the tidal horizontal-wind acceleration, we also provide the mean values of the tidal horizontal-wind acceleration averaged over 6 h time windows (see the dotted curves in Fig. 3). The 6 h time window (starting at 3 h before and ending at 3 h after each integer hour) is chosen for consistency with our calculation of the GW frequency. Based on Fig. 3, we note that at 3 h the tidal horizontal-wind acceleration above the center is positive, and the ground-based frequency in that region simultaneously increases with time. The tidal horizontal-wind acceleration below the center is very weak, so that the associated ground-based frequency remains almost unchanged. At 5 and 6 h, the tidal horizontal-wind accelerations are negative above but positive below the centers. Correspondingly, the frequency above the centers decreases with time while that below the centers increases with time. From 7 to 11 h, the tidal horizontal-wind acceleration is almost positive above the wave centers and negative below them. Therefore, the ground-based frequency increases (decreases) with time above (below) the centers. These quantitative comparisons of the tidal horizontal-wind acceleration and the gravity-wave frequency confirm that the time-variable tidal wind has a significant impact on the GW frequency and that the positive (negative) tidal horizontal-wind acceleration increases (decreases) the GW's ground-based frequency.
Generally, in Fig. 3 we can see that for the present case, after propagation for a significant length of time, the GW frequency near the first critical layer (above the centers) becomes much higher than the initial frequency, while at the lowest height shown (below the centers), the frequency becomes much lower than the initial frequency. This significant increase of the gravity-wave frequency near the first critical layer will lead to GWs penetrating the critical layer, which was shown in Fig. 1. Another possible cause for the penetration of GWs is that since the critical layer is specified based on an instantaneous tidal-wind field, in combination with tidal evolution this critical layer may move up or down. In other words, the critical layer is essentially transient. We will now provide an estimate of critical layer height. The first critical layer should be at a height where the tidal horizontal wind equals the critical layer threshold value determined by the local horizontal phase velocity of the GWs. Considering the tidal-wind acceleration-induced GW frequency variation, at 9 h the ground-based frequency of the GWs at 89.2 km is approximately 2.04 × 10 −3 rad s −1 , corresponding to a horizontal phase velocity of 34.6 m s −1 . This threshold value is much larger than the initial critical layer threshold value (17.3 m s −1 ). This estimate implies that, after 9 h of propagation, the transient first critical layer moves to 89.2 km, where the tidal horizontal wind is exactly 34.6 m s −1 . On the other hand, as shown in Fig. 1, at 9 h the tidal amplitude of the transient peak B (at 89.6 km) is 34.8 m s −1 . This is slightly larger than the GWs' horizontal phase velocity, implying that GWs can hardly penetrate peak B. These propagation scenarios agree well with the results shown in Fig. 2. Meanwhile, at 10 h the ground-based frequency of a GW at 87.2 km is also approximately 2.04 × 10 −3 rad s −1 , identical to that at 89.2 km at 9 h. The corresponding horizontal phase velocity of 34.6 m s −1 is larger than the tidal amplitude (33.2 m s −1 ) of the transient peak B (at 88.8 km) at this time. Thus, the transient first critical layer disappears and the GWs can penetrate peak B, which was also clearly shown in Fig. 2. Therefore, except for the tidal wind-induced transience of the critical layer, the increase of the ground-based frequency of GWs induced by the tidal wind also contributes significantly to the GWs penetrating the first critical layer.
To further explore the variability of the ground-based frequency of the GWs in a time-variable tidal field, we provide the height variability of the calculated GWs' ground-based frequency, as well as the tidal horizontal-wind acceleration between the first and second critical layers, and above the second critical layer; i.e., around the amplitude center positions of GWs that have penetrated the first and second critical layers. Since the GWs penetrated the first critical layer at about 10 h and the second critical layer at about 12 h, and we only showed 15 h of GW propagation, here we present the height variation of the GW frequency at 12 h (Fig. 4), which we calculated from the 6 h (9-15 h) time series using the same method as above. In the left-hand panel, similar to the discussion above, the tidal horizontal-wind acceleration is almost positive above the segment center and negative below it, thus leading to an increase of the GW frequency  near the second critical layer and eventually to penetration of the second critical layer (at approximately 110 km). It can be seen that the GW frequency is 0.95 × 10 −3 rad s −1 at 98.4 km while it is 2.04 × 10 −3 rad s −1 at 110 km. Clearly, the ground-based frequency at 98.4 km is much lower than that at 88 km (see Fig. 3), which illustrates that the groundbased frequency of the GW segment having penetrated the first critical layer decreases significantly due to local, large negative tidal horizontal-wind acceleration. As shown in the right-hand panel in Fig. 4, above the second critical layer the tidal horizontal-wind acceleration is very weak. Therefore, the frequency of the GW segment that has penetrated the second critical layer is 2.04 × 10 −3 rad s −1 (refer to the left-hand panel) and remains almost constant above the second critical layer. It should be noted that although the tidal wind above 124 km is unrealistic, since it was designed to attenuate exponentially to zero, our result clearly demonstrates that a weak tidal-wind field can hardly affect the GW frequency. Therefore, our simulation confirms that the timevariable tidal wind is the crucial factor leading to the variation of the GW frequency.
By performing a control study only with tidal temperature perturbation, i.e., T T calculated from GSWM-00 while u T , w T , and ρ T are set to zeros, we found that the time-varying tidal temperature perturbation has very little influence on the GW ground-based frequency.

Additional cases
To validate the analysis presented in the previous sections and clarify the influence of the initial phase of the tidal field on GW propagation, we carried out three additional case studies with the same parameters, except that the tidal field was delayed by 6, 12, and 18 h with respect to the first case. For convenience, we call our first exploration "case A" and these three additional cases B, C, and D. We found that in cases A, B, and C the GW is divided into three parts by two critical layers, but only into two parts by a single critical layer in case D. Figure 5 shows the height variations of the GWs' ground-based frequency (solid curve) and the tidal horizontal-wind acceleration (dotted curve) at 12 h between the first and second critical layers for cases B and C, and below the first critical layer for case D. Based on Fig. 5 we further confirm that, after propagation for a significant length of time, the GW frequency near the second or first critical layers becomes much higher than the initial frequency, which contributes significantly to GWs penetrating the critical layers. However, because of the different initial phases of the tidal field, the highest GW frequencies over the full height range shown are very different from each other; i.e., 2.04 × 10 −3 rad s −1 (Fig. 5, left), 2.69 × 10 −3 rad s −1 , 2.55 × 10 −3 rad s −1 , and 1.75 × 10 −3 rad s −1 for cases A, B, C, and D, respectively. In summary, the GW frequency is without a doubt manipulated by tidal-wind acceleration. More interestingly, although both increases and decreases of the GW frequency caused by time-variable tidal winds are possible, tidal-wind acceleration near the critical layers always causes the GW frequency to increase, irrespective of the initial phase of the tidal field. This leads to critical-layer penetration of GWs. Thus, when a GW propagates in a timevarying tidal-wind field, the tidal-wind acceleration is always favorable for the GW to penetrate its critical layers. The initial phase of the tidal field can only change the incremental quantity of the GW frequency near the critical layers.
A question arises from these simulations: although timevariable tides can exert both increasing and decreasing effects on the ground-based frequency of GWs propagating in a tidal-wind field, why does the tidal wind near the critical layer always cause a frequency increase of GWs and create favorable conditions for the GWs to penetrate the critical layers? The horizontal propagation direction of the GWs is naturally in the direction of the tidal wind around the critical layer, and the magnitude of the tidal wind increases with height just below the critical layer; i.e., at the lower edge of the critical layer. Combined with tidal evolution, the tidal zonal-wind phase front moves downward, leading to a positive acceleration of the tidal zonal wind at the lower edge of the critical layer along the GWs' horizontal propagation direction. This eventually causes the ground-based frequency of GWs propagating from below to increase.
To clearly demonstrate the time-variable tidal windinduced transience of the critical layer, and GW penetration from the perspective of wave energy, we provide the height distribution of the GW energy. First, we calculate the spatial distribution of the average GW energy density E d , which is calculated from the wave energy density E d by applying a low-pass filter with horizontal and vertical cut-off wavelengths of 240 and 12 km: E d = 1 2 (u 2 g + w 2 g ) + 1 2 g 2 T 2 g N 2 T 2 . In this equation, N is the buoyancy frequency, T the background temperature, and u g /w g /T g is the GW perturbation  component, which is filtered from the total perturbation component u /w /T by a high-pass filter with horizontal and vertical cut-off wavelengths of 180 and 9 km, respectively. From inspection of Fig. 2, we know that GWs deposit part of their energy at the first and second critical layers, and carry the residual energy to altitudes above the second critical layer.
The GW energy is divided into three segments at different height ranges. We calculated the height distribution of the GW energy density E h by integrating E d over the horizontal domain. Since a time-variable tidal field can lead to formation of the transient critical layers, the initial phase of the tidal field plays a dominant role in determining the heights of the transient critical layers. Figure 6 shows E h for cases A, B, C, and D. To clearly demonstrate that the GW energy could reach heights above 120 km, we present E h at 15 h for cases A, C, and D and at 16.5 h for case B (because of the late arrival time for case B). The energy peaks above 120 km show that some of the GW energy could penetrate their critical layers and be transported to heights above 120 km.
As for the heights of the critical layers, they can be clearly identified by the altitudes exhibiting a significant decrease of the energy density. However, in view of the unrealistic tidal wind above 124 km, we focus only on the critical layers below 120 km. Here, we used the initial (not the actual) critical layer heights determined by the horizontal phase speeds of the GWs and the horizontal tidal winds at the start time for comparison. They are located at 73. 6, 88.8, 82.4, and 76.8 km, respectively, for cases A, B, C, and D. Because of the time-variable tidal horizontal wind, as well as the variation of the GW frequency caused by the time-variable tidal wind, the heights of the actual critical layers, determined by the instantaneous value of the horizontal background wind and the local value of the horizontal phase velocity of the GWs, are not fixed; i.e., they are time-variable or transient. Just as mentioned above, there are two evident critical layers for cases A, B, and D but only one for case B after more than 10 h of propagation. At 15 h, the first/second critical layer is located at 82.8/102.4 km for case A and at 69.6/90.0 km for case C, and the only critical layer is located at 83.6 km for case D. For case B, at 16.5 h, the first/second critical layer is located at 79.2/92.8 km. Those heights are all very different from the initial critical layer heights. Liu et al. (2008) pointed out that the GW penetration height increases with vertical wavelength. However, our simulation results show that even for GWs with definitive parameters, the critical layer is changeable and transient, and GWs could reach higher altitudes above the initial critical layer. Moreover, for a specific GW, the height of the critical layer sensitively depends on the initial tidal phase.
Our results agree well with the conclusion of Zhong et al. (1995) that the temporal variation of the tidal component of the wind changes the observed frequency, sometimes substantially, and that different start times relative to the tidal phase cause the temporary critical levels to appear at different heights.

Dissipation effects
In these simulations, we did not take dissipation into account in order to clearly show that the wave energy leaks from critical layers and is transmitted to higher altitudes. Here, we add two more realistic cases (E and F). In case E, only eddy diffusion is introduced, while in case F eddy and molecular diffusion are both considered. The formulations assumed and the coefficients adopted for the eddy and molecular diffusion of heat and momentum are the same as those used by Hagan et al. (1993). Now, the governing equations are where µ 0 is the dynamic molecular viscosity coefficient; ν eddy is the kinematic eddy viscosity coefficient; K 0 is the molecular thermal conductivity coefficient; and K eddy is the eddy thermal conductivity coefficient. Obviously, eddy diffusion dominates below the mesopause, while molecular diffusion becomes important at a height of approximately 100 km. For comparison, in Fig. 7 we present the propagation images for cases A (no dissipation), E (only eddy diffusion), and F (eddy and molecular diffusion). We see that, in the presence of diffusion, GW energy absorption would be noticeable and the GW packet less sharply concentrated near critical layers. As regards GW propagation, molecular diffusion has little impact below 100 km but great impact above that altitude. The weaker the dissipation is, the more wave energy is transferred to higher altitudes. Note that in the absence of dissipation, there also exists wave energy leakage of critical layers at the altitudes above 120 km. However, this could not be displayed because the maximum absolute value is more than 5 times larger than those in the presence of dissipation. Our results are consistent with those obtained by Walterscheid (1981). In conclusion, although dissipation can significantly affect GW propagation in a time-varying tide, it cannot eliminate the transient nature of the critical layers and prevent GWs from penetrating the critical layers and reaching higher altitudes. We also explored the height and time variations of the GW ground-based frequency below the critical layer for cases E and F. As they are almost the same, only the results for case F are presented in Fig. 8. Although the variation of the ground-based frequency is somewhat different from that for case A, the main conclusion regarding the frequency variations of GWs interacting with a time-varying tide remains unchanged; i.e., the ground-based GW frequency changes significantly and the frequency near the critical layer always increases. This is favorable for GWs penetrating their critical layer.

Conclusions and comments
Aiming at investigating nonlinear interactions between GWs and tides, we simulated the propagation of GWs in timevariable tidal winds using a fully nonlinear model. In this numerical study, we focused mainly on the impact of tidal-GW interactions on GW propagation, in particular on the GW ground-based frequency and the critical layers. Senf and Achatz (2011) previously showed that the temporal change of the background wind due to tides enables a GW to avoid a critical level by changing its ground-based frequency and phase speed. This investigation was based on a ray tracing model; i.e., a WKB method. The current paper uses an explicitly non-WKB method to study this effect.
To specify simple, well-controlled conditions, several simplifications were implemented: (1) we force the model to assume the background wind and temperature structure of a tide, hence full nonlinear interactions with the waves are suppressed; (2) in contrast to ray tracing technique, a wave packet in the nonlinear model has a finite vertical extent and, accordingly, the frequency of the GW varies over this range. In addition, the limited vertical resolution of the model will affect the results for shorter vertical wavelengths, in particular when the phase speed reaches a value close to the background wind.
Our simulations show that when a GW packet propagates in a tidal-wind environment, not only the intrinsic frequency but also the ground-based frequency of the GW packet changes significantly. The tidal horizontal-wind acceleration dominates the GW frequency variation. Positive (negative) accelerations induce the frequency to increase (decrease) with time. Compared with the tidal winds, the nonlinearity-induced background wind has only negligible effects on GW propagation characteristics.
Although both increases and decreases of the GW frequency caused by time-variable tidal winds are possible, tidal-wind acceleration near the critical layers always causes the GW frequency to increase, no matter the initial phase of the tidal field.
For GWs propagating in a tidal-wind field, the combination of their increased ground-based frequency and the transient nature of the critical layer induced by the time-variable tidal zonal wind create favorable conditions for them to penetrate their initial critical layer.
To verify the general applicability of our results, we calculated four additional cases (G, H, I, and J), characterized by almost the same parameters as case A, except that the horizontal or vertical wavelength was changed. For cases G, H, I, and J, the horizontal/vertical wavelengths are 60/6 km, 240/6 km, 120/3 km, and 120/9 km, respectively. The simulation results show that, although the horizontal or vertical wavelength changes over a wide range, the primary conclusions remain unchanged (figures not shown): the groundbased GW frequencies near the critical layers increase significantly; the increased frequency and the transient nature of the critical layers induced by a time-variable tidal zonal wind are combined to create favorable conditions for energy leakage of critical layers. That is to say, for the GWs in this range of scales, i.e., 60-240 km in the horizontal wavelength and 3-9 km in the vertical wavelength, the ground-based GW frequencies are greatly affected by a time-variable tidal wind and part of wave energy can penetrate the expected critical layer.
Previous studies (Liu et al., 2008) have suggested that the critical layer induced by a tidal wind could almost completely absorb the GW energy and prevent GWs from propagating to greater heights. However, this assessment was based on simulations in a frame of time-invariant tidal-wind field. We note that our wave-energy analysis demonstrates that part of the GW packet could penetrate the expected critical layer and carry energy to higher altitudes. That is to say, GWs would reach much higher altitudes than expected and accelerate the background atmosphere there. This indicates that dynamical effects induced by tidal-GW interactions are more complicated than commonly taken into account by GW parameterizations in global models. These complex effects should be included in GW drag parameterization schemes to calculate a more realistic dynamical structure of the middle and upper atmospheres.
On the other hand, GW perturbations in the mesosphere are believed to be excited mainly in the lower atmosphere. However, decades of observations have revealed that GW perturbations in the lower and mid-upper atmosphere are, respectively, dominated by low-frequency Zhang et al., 2010Zhang et al., , 2012 and high-frequency (Diettrich et al., 2005;Dou et al., 2010;Li et al., 2011) components. As discussed in this paper, the evidently higher ground-based frequency of GWs observed in the mesosphere may result from the significant time-variable tidal windinduced frequency increase of low-frequency GWs excited in the lower atmosphere.