The importance of friction in mountain wave drag amplification by Scorer parameter resonance

A mechanism for amplification of mountain waves, and their associated drag, by parametric resonance is investigated using linear theory and numerical simulations. This mechanism, which is active when the Scorer parameter oscillates with height, was recently classified by previous authors as intrinsically nonlinear. Here it is shown that, if friction is included in the simplest possible form as a Rayleigh damping, and the solution to the Taylor–Goldstein equation is expanded in a power series of the amplitude of the Scorer parameter oscillation, linear theory can replicate the resonant amplification produced by numerical simulations with some accuracy. The drag is significantly altered by resonance in the vicinity of n/l0 = 2, where l0 is the unperturbed value of the Scorer parameter and n is the wave number of its oscillation. Depending on the phase of this oscillation, the drag may be substantially amplified or attenuated relative to its non‐resonant value, displaying either single maxima or minima, or double extrema near n/l0 = 2. Both non‐hydrostatic effects and friction tend to reduce the magnitude of the drag extrema. However, in exactly inviscid conditions, the single drag maximum and minimum are suppressed. As in the atmosphere friction is often small but non‐zero outside the boundary layer, modelling of the drag amplification mechanism addressed here should be quite sensitive to the type of turbulence closure employed in numerical models, or to computational dissipation in nominally inviscid simulations. Copyright © 2012 Royal Meteorological Society


Introduction
Gravity wave drag is one of the key phenomena that must be parametrized in global weather prediction and climate models. A multitude of processes affect this topographically generated force, making it exceed leadingorder estimates from linear theory, in what are called 'highdrag states'. It is important to understand drag enhancement mechanisms, since they provide a dominant contribution to the globally integrated drag, having a significant impact on the deceleration of the atmospheric circulation.
In a recent paper, Wells and Vosper (2010) (hereafter referred to as WV10) assessed the accuracy of linear theory for diagnosing mountain wave drag in stratified flow over 2D ridges. The focus of their study was on situations where the actual drag might significantly exceed the drag predicted by linear theory, even for extremely low mountain heights. One such situation, identified by them, is when the spectrum of the Scorer parameter profile contains harmonics that can lead to a resonant amplification of the mountain waves. In order to isolate this effect, they considered an idealized case where a small sinusoidal perturbation is added to a constant Scorer parameter profile, with a wavelength half that of the dominant vertical wavelength of the internal gravity waves. In that case, they showed that the drag may exceed its linear value by a factor of 2 or more, depending on the phase of the Scorer parameter oscillation. Since the linear model of WV10 failed to predict this behaviour, they attributed the drag amplification to nonlinear wave-wave interactions, as investigated originally by Phillips (1968) in an oceanographic context, and more recently by Nance and Durran (1998) and Lee et al. (2006). However, although in WV10's study the nonlinear drag was normalized by the corresponding linear value, the behaviour of the linear drag was never explicitly shown for this case.
The characteristics of the drag amplification outlined above immediately suggest that it results from parametric resonance, since the role played by the Scorer parameter in the Taylor-Goldstein equation is akin to that of the coefficient multiplying the position in an equation describing a simple harmonic oscillator. This resonance is therefore of a different kind from those investigated, for example, by Miranda and Valente (1997), Wang and Lin (1999), Leutbecher (2001), and Teixeira et al. ( , 2008, which resulted from discontinuities in the mean atmospheric parameters, or their derivatives. It is also different from the resonance investigated by Grubišić and Stiperski (2009) and Stiperski and Grubišić (2011), which results from the horizontal distribution of the topographic forcing of lee waves (see also Grisogono et al., 1993;Vosper, 1996).
One aim of the present paper was to address an idealized situation akin to that considered by WV10 using linear theory, thus showing that the kind of resonance they investigated is possible under its assumptions, and can lead to very substantial drag enhancement. A second aim was to show how the inclusion of friction in the model is crucial to obtain a drag behaviour qualitatively similar to that produced in the numerical simulations of WV10, or roughly similar ones. Friction fulfils two roles: firstly, it limits the drag magnitude in resonant conditions (something that is presumably effected by nonlinear processes in more realistic circumstances); secondly, since friction prevents the singular behaviour of the drag that occurs in inviscid conditions, it also widens the drag extrema at resonance, rendering them detectable in a representation similar to that of Figure 9 of WV10.
We therefore speculate that inviscid linear theory, such as employed by WV10 in part of their calculations, should be unable to properly represent this type of resonance. For the same reasons, it is likely that the magnitude of the drag enhancement in resonant conditions is quite sensitive to dissipative processes in numerical models, be they due to the type of adopted turbulence closure, or to the more or less diffusive character of the discretization scheme employed. These conjectures are tested in the present study by comparing results from linear theory, where a small sinusoidal variation of the Scorer parameter is treated using a perturbation approach, and fully nonlinear numerical simulations of the same situation.
This paper is organized as follows: section 2 contains a description of the analytical theoretical model and of the numerical simulations. In section 3, the variation of the drag with friction and non-hydrostatic effects is explored, and the pressure and velocity fields are analysed. Results from the theoretical model and from the numerical simulations are compared. Finally, section 4 summarizes the main findings of this study.

Analytical drag model
Stationary, non-rotating flow over a 2D mountain ridge is considered. The Boussinesq approximation is assumed, and the equations of motion are linearized with respect to a reference state, but the flow is assumed to be nonhydrostatic. Friction, which turns out to be a crucial effect, is included in the simplest possible form as a Rayleigh damping. The adopted equations of motion (cf. Lin, 2007, sections 5.1 and 13.2.2) are Here x is the horizontal coordinate perpendicular to the ridge, z is the height, U(z) is the incoming wind velocity (aligned with x), N(z) is the Brunt-Väisälä frequency of the incoming flow, ρ 0 is a reference density (assumed to be constant) and λ is a (constant) Rayleigh damping coefficient. u, w, p and b are, respectively, the horizontal velocity, vertical velocity, pressure and buoyancy perturbations associated with the mountain waves. Note that, for simplicity, friction terms are only included in the momentum equations and not in the heat balance equation. It can be shown that, for the small values of λ to be considered, this has a negligible effect on the model behaviour apart from a rescaling of λ by a factor of 2 (if the Rayleigh coefficient for heat was assumed to be the same).
Equations (1)-(4) may be differentiated and combined, yielding a single equation for w. All perturbed variables may be expressed as 1D Fourier integrals. For example, w can be written (see, for example, Lin, 2007, Appendix 5 whereŵ is the corresponding Fourier transform, k is the horizontal wave number and i = √ −1. From Eqs (1)-(4) and (5), it can be shown thatŵ satisfieŝ where l 2 = N 2 /U 2 − U /U is the square of the Scorer parameter, and the primes denote differentiation with respect to z. In order to reproduce conditions akin to those considered in Figure 9 of WV10, the Scorer parameter squared is assumed here to take the form where l 2 0 is a constant, ε is a small dimensionless parameter, n is the vertical wave number of the perturbation imposed on l 2 0 and φ is the corresponding phase. Equation (7) defines a Scorer parameter that oscillates with height with a relatively small amplitude (see Figure 1).
For the particular case λ = 0, Eq. (6) along with Eq. (7) is a Mathieu equation, which describes parametric resonance, and its solutions must be expressed in terms of Mathieu functions (Abramowitz and Stegun, 1972). However, when ε is assumed to be small (it takes the value 0.1 in WV10), it is possible to solve Eq. (6) approximately in terms of elementary functions using a perturbation approach, by expanding its solution in a power series of ε: In fact, for the present purposes, and as will be seen next, it is sufficient to consider this series only up to first order in ε.
If Eqs (7) and (8) are introduced into Eq. (6), the following two equations result at zeroth and first order in ε: In this formulation, it becomes especially clear that the zeroth-order solutionŵ 0 , in conjunction with the perturbation to the Scorer parameter, acts as a source term for the first-order solutionŵ 1 . Thus this equation set already contains the possibility of resonance.
The solution to Eq. (9) is (see Lin, 2007, section 5.2.1) where U 0 = U(z = 0),ĥ is the Fourier transform of the surface elevation and m is the vertical wave number of the internal gravity waves, defined, using Eqs (9) and (11), as In the above passage it was implicitly assumed that U is constant, otherwise Eqs (11) and (12) would not be valid in the generic case λ = 0. However, if λ = 0, this assumption is not necessary, so U will continue to be treated as a function of z in the following, for maximum generality. The coefficient multiplying the exponential in Eq. (11) takes into account the boundary condition at the surfaceŵ 0 (z = 0) = iU 0 kĥ and the sign of m is determined by the upper radiation boundary condition. If m is decomposed into its real and imaginary parts m R and m I m = m R + im I (which are non-zero simultaneously as long as λ = 0), then it can be shown that the definitions of m R and m I that are physically consistent (see Appendix A) are which result from Eqs (12) and (13). In Eq. (15) m I > 0, as it must be for the wave perturbation to decay with height according to Eqs (11) and (13). On the other hand, in Eq. (14) m R has the same sign as Uk, which corresponds to upward wave energy propagation (see, for example, Miranda, 2005, 2006). Inserting Eq. (11) into Eq. (10), the solution to the latter equation (see Appendix B) iŝ where it has again been assumed that U is constant when λ = 0. This form emphasizes upward and downward propagating components of the first-order wave perturbation (corresponding to exponential terms with positive and negative exponents that are a function of z, respectively). The solution (Eq. (16)) satisfies the lower boundary condition w 1 (z = 0) = 0 and the upper boundary condition that the wave energy decays or propagates upward as z → +∞. It can be seen from Eq. (16) that the upward and downward propagating wave components are of similar magnitude at z = 0, whereas the downward propagating component vanishes as z → +∞. The integrals in Eq. (16) may be calculated analytically, yieldinĝ In this equation it can already be seen that parametric resonance will occur for n = ±2m R , since 4m 2 − n 2 appears in a denominator on the right-hand side. When λ = 0, m 2 is real and there is the possibility that this denominator becomes zero. When λ = 0, however, Eq. (15) shows that m I is never zero, even if the waves are vertically propagating. Thusŵ 1 will not diverge to infinity in those circumstances (which would invalidate the power series solution), but will be strongly amplified for relatively small λ.
The focus in the present study is on the calculation of the surface drag associated with the mountain waves. The drag per unit spanwise width of the ridge is given by (Teixeira and Miranda, 2004) where h is the surface elevation, the asterisk denotes complex conjugate, and Parseval's theorem has been used. In order to calculate the drag, it is therefore necessary to obtain the pressure perturbation at the surface. Using Eqs (1) and (4), the Fourier transform of the pressure perturbationp, which is related to p in the same way asŵ is related to w in Eq. (5), is given byp An inviscid (λ = 0) but 3D version of Eq. (19) was presented, for example, by Teixeira and Miranda (2006) as their Eq. (8). From Eqs (8) and (19), it is clear thatp may also be expressed as a power series of ε, aŝ If Eq. (11) is differentiated and evaluated at z = 0, and used in Eq. (19), it can be shown (see Appendix C) that where U 0 = U (z = 0), and Eqs (13) and (20) were also used. Equation (17) may also be differentiated, evaluated at z = 0, and inserted into Eq. (19) (see Appendix C), yieldinĝ where, again, Eqs (13) and (20) have been used.
In the case of an orography that is symmetric in x,ĥ is real and so, according to Eq. (18), only the imaginary parts ofp 0 andp 1 contribute to the drag. The drag is obviously also expressed as a power series in ε: Here the total drag, which is only calculated up to first order, is normalized by the drag in the absence of resonance, D 0 . This is accomplished by inserting Eqs (21) and (22) into Eq. (20) and using the latter equation in Eq. (18). The final result is where k = ka,ĥ =ĥ/(h 0 a), m R = m R /l 0 , m I = m I /l 0 and n = n/l 0 are dimensionless parameters and functions defined in terms of the corresponding dimensional quantities, specified previously. a and h 0 are, respectively, a representative width and height of the orography.
Although the results of the model would not be essentially changed by adopting a different form for the surface elevation, following WV10 it will be assumed here that the orography is a bell-shaped ridge: The normalized drag given by Eq. (24) is a function of five dimensionless parameters n/l 0 , φ, ε, l 0 a and λa/U. Since, as was seen above, resonance occurs when n = ±2m R , and |m R | ≈ l 0 when the flow is approximately inviscid and hydrostatic (i.e. when l 0 a is large and λa/U is small-see Eq. (14)), resonance will occur in the vicinity of n/l 0 = 2. For that reason, the drag will be represented next as a function of n/l 0 , for particular values of the other flow parameters.

Numerical simulations
Numerical simulations were performed using the FLEX numerical model, which is a 2D nonlinear and nonhydrostatic microscale to mesoscale model using generalized curvilinear coordinates (see Arga´ın et al., 2009). For all simulations, and as in WV10, the rotation of the Earth was neglected and flow over a bell-shaped mountain (Eq. (25)) with a = 10 km, h 0 = 10 m was considered. The unperturbed Brunt-Väisälä frequency and the mean wind speed were N 0 = 0.01 s −1 and U = 20 m s −1 . Therefore, the flow was strongly linear (N 0 h 0 /U = 5 × 10 −3 ) (at least in non-resonant conditions) and reasonably hydrostatic (N 0 a/U = 5). When testing non-hydrostatic effects, a = 20 km and a = 4 km, were also employed, so that N 0 a/U = 10 (strongly hydrostatic flow) or N 0 a/U = 2 (strongly nonhydrostatic flow). The adopted Scorer parameter profile was constructed according to Eq. (7), with l 0 = N 0 /U and its vertical oscillation was imposed by adding a perturbation to N 2 0 , of relative amplitude ε = 0.1. This approach differs from that of WV10, where the perturbation is imposed on U instead.
The model was run in a domain of length 24a in the horizontal (≈ 240 km for a = 10 km) by 3.5λ z in the vertical (λ z = 2π U/N 0 is the vertical wavelength of the gravity waves). This corresponds to ≈ 44 km for the values of N 0 and U quoted above. A grid of 160 × 525 points, with a resolution of 0.15a in the horizontal (1.5 km for a = 10 km) and λ z /150 in the vertical (≈ 84 m for the values of N 0 and U quoted above) was used. Lateral sponges with width 4a and a sponge at the top of the domain with depth 1.5λ z were also employed.
The model was run with a time step of 2 s for a period of up to 500a/U (corresponding to 69 h for a = 10 km and U = 20 m s −1 ), until the drag stabilized to a constant value. In many cases, the time necessary for this to happen was not larger than 120a/U, or 17 h, but in resonant conditions it became considerably larger (cf. WV10).
Most of the simulations were carried out in inviscid mode, i.e. using a free-slip boundary condition at the surface and no turbulence closure. The simulations of section 3.5 used a no-slip boundary condition, and either a Smagorinskytype turbulence closure (Lilly, 1962) or a K-ε turbulence closure. In these cases, Monin-Obukhov profiles were used to initialize the model up to the top of the boundary layer height (which is consistent with the absence of rotation). A roughness length of z 0 = 0.05 m, a Monin-Obukhov length of L MO = 255 m and a friction velocity of u * = 0.46 m s −1 were employed. The latter two values were determined using the method proposed by Arga´ın et al. (2009). In these simulations, the drag took considerably less time to stabilize than in inviscid conditions (at most 150a/U or 21 h in all cases).

The importance of friction
Since friction appears to be of crucial importance in the type of flow being considered, its effect on the behaviour of the analytical model presented above will be analysed first, and compared with nominally inviscid numerical results. Figure 2 shows the normalized drag, as given by Eq. (24), for ε = 0.1 and l 0 a = 5 (as in WV10), for three values of the friction coefficient λa/U (lines). Also shown are results for the same conditions, but from inviscid simulations of the FLEX numerical model (symbols). Although, in all cases displayed, the normalized drag only departs substantially from 1 in the vicinity of n/l 0 = 2, which is consistent with the existence of parametric resonance, the behaviour is strongly dependent on φ. For φ = 0 (Figure 2(a)), there is a drag maximum to the left of n/l 0 = 2, followed by a minimum to the right of n/l 0 = 2. The magnitude of these extrema increases as λa/U decreases. In Figure 2(b), where φ = π/2, there is a single drag minimum, of larger magnitude than the extrema displayed in Figure 2(a), approximately centred on n/l 0 = 2 (slightly to the left). The magnitude of this minimum increases as λa/U decreases, but the minimum vanishes altogether for λ = 0. For φ = π (Figure 2(c)) and φ = 3π/2 (Figure 2(d)) the behaviour of the drag predicted by Eq. (24) is exactly symmetric with respect to the line D/D 0 = 1 relative to that displayed in Figure 2(a) and (b), respectively. Hence for φ = π there is a drag minimum followed by a maximum, whereas for φ = 3π/2 there is a single maximum. This last maximum, which is perhaps the most relevant result from the viewpoint of drag parametrization, has the same peculiar behaviour as the minimum. Namely, while it becomes higher as λa/U decreases, it vanishes for λ = 0. It should be noted that the drag minima predicted by Eq. (24) at φ = 0, φ = π/2 or φ = π would become negative for sufficiently low λa/U, an obviously unphysical result which stems from the limitations of the perturbation approach employed in the analytical model.
The results from the FLEX numerical model are broadly in agreement with the analytical model for λa/U = 2 × 10 −2 . They cover a much more limited number of values of n/l 0 , and for that reason are not able to sample as effectively the drag variation near the peaks. For example, since the drag minimum in Figure 2(b) and the drag maximum in Figure 2(d) occur slightly below n/l 0 = 2, they are not perfectly captured in the numerical simulations. Unlike in the analytical model, the magnitude of the drag minimum in Figure 2(b) is smaller than the magnitude of the drag maximum in Figure 2(d), essentially because the drag cannot become negative in the numerical simulations. Nevertheless, apart from an underestimation of the drag by Eq. (24) for φ = 0 at n/l 0 = 1.75, the agreement between both models is quite good, even quantitatively. This means that a simple Rayleigh damping is able to capture the essential aspects of frictional effects in the present flow. It also means that, while being nominally inviscid, the simulations of the FLEX model obviously contain a small amount of numerical diffusion. Preliminary tests, where the FLEX model was run at a higher resolution, suggest that the drag extrema existing near n/l 0 = 2 become more pronounced, which is consistent with a decrease of computational dissipation, but further tests are required to confirm the robustness of this behaviour.
The drag maximum in the inviscid numerical results of WV10, presented in their Figure 9, is slightly larger than that in Figure 2(d), but these results are not directly comparable with ours, since WV10 imposed a perturbation on U rather than on N to produce the vertical oscillations of the Scorer parameter (see section 4).
The present results suggest that friction, which always exists in reality, is essential for the existence of the drag maximum for φ = 3π/2 and of the drag minimum for φ = π/2. Curiously, that is not the case for the double extrema produced for φ = 0 or φ = π , which were overlooked by WV10, but are essentially inviscid in the present framework (see Figure 2(a) and (c)). The possibility that a modulation of the drag with n/l 0 may be produced when λ = 0 for φ = π/2 or φ = 3π/2 if the wave solution is considered up to higher order in ε cannot be ruled out, but when ε = 0.1 this should amount to a relatively small correction, since the following power series coefficient is ε 2 = 0.01.
The mathematical reason for the behaviour of the drag in the analytical model can be sought in Eq. (24). When φ = π/2 or φ = 3π/2, cos φ = 0 and so only the second term in the numerator of the fraction inside the integral is non-zero. This term contains m I and is multiplied by m R outside the fraction, so it only gives a non-zero contribution to the integral when m R and m I are simultaneously nonzero. Now, in inviscid conditions, either m R or m I must be zero, i.e. the vertical wave number of the mountain waves must either be pure real or pure imaginary. In those circumstances, the correction to the drag due to resonance vanishes.
A clearer explanation for this behaviour can be achieved if one notes that for the case φ = 3π/2 (corresponding to a single drag maximum) the integral in the numerator of Eq. (24) can be written where Eqs (14) and (15) have been used. Obviously when λa/U → 0 this integral is zero if n/l 0 = 2, but it can be shown that it diverges if n/l 0 = 2. This becomes even more evident if the hydrostatic limit is considered. In that case, and it follows immediately that the integral in Eq. (27), apart from approaching zero (making D/D 0 = 1) when λa/U → 0 and n/l 0 = 2, tends to infinity proportionally to (λa/U) −1 if n/l 0 = 2. In fact, it can be shown from Eqs (24) and (27) that, in the latter limit: Thus the correction to the drag due to parametric resonance behaves in this case like a Dirac delta function as λa/U → 0, although the perturbation approach used to obtain this result becomes invalid in that limit. Obviously, when frictional effects are excluded from the outset, this behaviour is not uncovered, since the limit λa/U → 0 is taken before the limit n/l 0 → 2.
The linear model used by WV10 was based on an equation analogous to Eq. (6), but with λ = 0 (their Eq. (3)). The above arguments could explain the inability of these authors to obtain a drag maximum such as shown in Figure 2(d) (or a minimum such as that in Figure 2(b)) in their linear results. However, caution is necessary, since their model, although linear, was solved numerically, and the discretization scheme should introduce some spurious numerical diffusion. Additionally, as noted above, WV10 perturbed the Scorer parameter in a different way from that adopted here, which makes comparisons difficult.
Even an inviscid linear model should be able, however, to produce the double drag extrema displayed in Figure 2(a) and (c), since when φ = 0 or φ = π , sin φ = 0 and so only the first term in Eq. (24) in the numerator of the fraction inside the integral is non-zero, but this term does not vanish if m I = 0.
The disappearance of the single drag extrema for φ = π/2 or φ = 3π/2 is an example of a situation where λ → 0 (vanishing friction) does not coincide with λ = 0 (zero friction). Other situations of this kind exist in fluid dynamics, in connection with boundary layer theory, e.g. D'Alembert's paradox. Since the real atmosphere always possesses some friction, the subtle difference between these two limits is not very relevant from an experimental perspective. In the context of numerical modelling, however, it illuminates the fact that the behaviour of the drag in these resonant flows must be quite sensitive to computational diffusion in inviscid simulations, and certainly requires an adequate turbulence closure to be accurately represented (see section 3.5 below).

The pressure perturbation
In order to better understand the behaviour of the drag, it is worth analysing the pressure perturbation at the surface, which is ultimately responsible for it. The cases of greatest interest are those where a single drag maximum or a single drag minimum exists, because they presumably correspond to extreme flow configurations. Figure 3 shows the normalized pressure perturbation at the surface as a function of streamwise distance across the ridge from the analytical model (Figure 3(a)) and from the FLEX numerical model (Figure 3(b)). Cases corresponding to a drag maximum (n/l 0 = 2 and φ = 3π/2) and to a drag minimum (n/l 0 = 2 and φ = π/2) are shown, with other flow parameters similar to those used by WV10 (l 0 a = 5 and ε = 0.1). In the analytical model, λa/U = 2 × 10 −2 is assumed. A case without resonance (ε = 0) is also presented for reference. This last case is characterized (both in Figure 3(a) and (b)) by a pressure distribution that is anti-symmetric with respect to the ridge axis, with a positive maximum upstream and a negative minimum downstream of the ridge, as is typical of approximately hydrostatic flow over this kind of orography (see, for example, Queney, 1948). Although the flow is not perfectly hydrostatic here, non-hydrostatic effects on the surface pressure are almost imperceptible, especially in Figure 3(a). The situation of drag enhancement is characterized by a substantial amplification of the pressure perturbation, which nevertheless maintains a roughly anti-symmetric form. The situation of drag attenuation, on the other hand, is characterized by a decrease of the positive pressure maximum upstream of the ridge and a translation of the negative minimum upstream towards the ridge top (without appreciable reduction in magnitude).
In all of these aspects, there is considerable agreement between the analytical model (Figure 3(a)) and the numerical simulations (Figure 3(b)). There are, however, some slight discrepancies, particularly in the resonant cases. For example, in the high-drag case the positive pressure perturbation extends more upstream of the ridge in the numerical simulations than in the analytical model, and the negative pressure perturbation extends less downstream and is centred further upstream. Conversely, in the analytical model, the pressure perturbation is positive over the ridge top and is almost zero (as in the non-resonant case) in the numerical results. Additionally, in the low-drag case in the analytical model the pressure perturbation downstream of the ridge oscillates and has a positive maximum downstream of the minimum, which does not exist in the corresponding numerical result. This feature, which appears to be a nonhydrostatic effect, might reflect the unreliability of the analytical model for reproducing the details of the pressure distribution in resonant conditions, due to the substantial amplitude of the flow perturbations, despite the smallness of the forcing. Indeed, it is this maximum of the pressure perturbation that, for a sufficiently high ε, is responsible for making the drag become negative, so it is likely a spurious feature. Nevertheless, and as would be expected considering the drag behaviour, the essential effects of resonance on the pressure field seem to be captured by the analytical model for all three cases.
Given the peculiar behaviour of the drag produced by the analytical model in the inviscid limit, it might appear relevant to present plots of the pressure perturbation for this case in the same resonant conditions as illustrated in Figure 3. However, in the inviscid limit the first-order term of the Fourier transform of the pressure perturbation has a singularity at n/l 0 = 2, which makes the corresponding inverse Fourier transform diverge. This is not surprising since the drag amplification, if it existed in this limit, would be infinite. However, the diverging pressure field (not shown) is perfectly symmetric with respect to the orography, producing no additional drag, consistent with Figure 2(b) and 2(d). For these reasons, this inviscid pressure perturbation is not presented here.

The flow field
A better understanding of the behaviour of the pressure perturbation can be achieved by analysing the velocity field associated with it. From Eq. (19) it is clear that the pressure perturbation is determined by the structure of the vertical velocity perturbation. In particular, it can be shown that the term in this equation that contributes to the pressure that produces drag is that proportional to the vertical derivative ofŵ.
In Figure 4, the normalized vertical velocity perturbation w/(Uh 0 /a) is shown in a vertical cross-section perpendicular to the ridge, as calculated from the analytical model for ε = 0.1, l 0 a = 5 and λa/U = 2 × 10 −2 . The three situations displayed in Figure 4(a), (b) and (c) are the same as those considered in Figure 3. Namely, in Figure 4(a) a non-resonant case with ε = 0 is considered. Figure 4(b) displays a case of drag enhancement with ε = 0.1, n/l 0 = 2 and φ = 3π/2, while Figure 4(c) displays a case of drag attenuation with ε = 0.1, n/l 0 = 2 and φ = π/2. In all panels, the structure typical of propagating mountain waves can be seen, with elongated maxima and minima over the mountain, tilted upstream. Non-hydrostatic effects are visible, with some downstream propagation of the wave pattern and some attenuation as one moves upwards. In Figure 4(a) (the non-resonant case) the negative lobe of the vertical velocity sitting directly above the ridge has a minimum value below −0.9, while in Figure 4(b) (the high-drag state) that minimum is lower than −1.3 and in Figure 4(c) (the low-drag state) the minimum is merely below −0.6. Obviously, these differences in magnitude of the vertical velocity explain the differences in the pressure perturbation described in the previous section. There are some additional differences. In Figure 4(b) the pressure perturbation extends somewhat more in the upstream and downstream directions than in Figure 4(a). In particular, the non-hydrostatic wake of the mountain waves is considerably enhanced. In Figure 4(c), on the other hand, the maxima and minima of the vertical velocity have a two-lobe structure. Figure 5 shows results similar to those of Figure 4, but produced by the FLEX numerical model. The minima of the vertical velocity perturbation above the ridge are just below −0.8 in Figure 5(a), below −1.2 in Figure 5(b) and below −0.6 in Figure 5(c). Concerning this aspect, there is thus quite good agreement with Figure 4. However, the configuration of the flow shows some differences. For example, in Figure 5(b) the vertical velocity perturbation extends less upstream of the ridge than in Figure 4(b), and the way in which it extends downstream is rather  conditions are considered, are perhaps the most similar, which is expectable, since they should correspond to the least nonlinear flow. Unlike in Figure 4, in all panels of Figure 5 the flow perturbation tends to decay markedly above l 0 z/π = 4. This happens because the sponge applied in the numerical simulations at the top of the domain extends approximately down to that height. There is also more moderate attenuation below l 0 z/π = 4, which is clearly larger than in the analytical results, and a tendency for the flow pattern to become more horizontal. This could result from numerical dissipation associated with the discretization scheme. Either effect does not seem to have a large impact on the flow near the surface, as would be desirable for the present purposes.
It would be interesting to visualize the vertical velocity field given by the analytical model in invicid and resonant conditions. Unfortunately, for the same reasons as invoked for the pressure perturbation, the corresponding field diverges, and so is not presented here. The first-order term of this field, which is responsible for this divergence, does not display any tilting in its vertical structure (not shown), being therefore unable to increase the drag.

Non-hydrostatic effects
WV10 pointed out that, when the value of U is decreased, the drag maximum displayed in their Figure 9 for n/l 0 ≈ 2 increased in magnitude. The present analytical model suggests that this behaviour may be explained by the variation of l 0 a, which increases as U decreases or as N 0 increases. This motivates an exploration of the drag dependence on this parameter.
In Figure 6 the normalized drag predicted by Eq. (24) is presented as a function of n/l 0 for ε = 0.1 and λa/U = 2 × 10 −2 , for the same values of φ as before and three values of l 0 a: 2, 5 and 10. It can be seen that the drag maxima and minima become more pronounced and narrower as l 0 a increases. For example, in Figure 6(d) the drag maximum has a magnitude of ≈ 2.6 when l 0 a = 10. This is still somewhat lower than the value of 3.5 predicted by Eq. (28) for perfectly hydrostatic conditions. By symmetry, this implies that in Figure 6(b) the drag minimum is negative for l 0 a = 10, an unphysical result which stems for the limitations of the analytical model, as pointed out before. For l 0 a = 2, a case where the flow is highly non-hydrostatic, the minimum in Figure 6(b) and the maximum in Figure 6(d) are decreased, translated to slightly lower values of n/l 0 , and become wider on their left slopes. This should be caused by wave dispersion, which tends to 'unfocus' the resonance process. Since resonance relies on trapping by vertical reflections of the wave energy, dispersion attenuates it by allowing that energy also to propagate downstream instead of becoming only concentrated over the mountain. Figure 7 shows similar results, but from the FLEX numerical model. Clearly, the main features in the drag variation with n/l 0 are in agreement with those of Figure 6; nevertheless, because the drag is only plotted as a function of n/l 0 at intervals of 0.25, some details are necessarily lost. For example, the increase in the magnitude of the drag extrema as l 0 a increases is less pronounced. This happens because, despite the fact that the drag minimum in Figure 7(b) and the drag maximum in Figure 7(d) become centred closer to n/l 0 (for which there is a computed value of D/D 0 ), their width becomes smaller. Therefore, the magnitude of the drag extrema seems to be progressively more underestimated as l 0 a becomes larger. A similar phenomenon can be seen for the drag minimum in Figure 7(a) and 7(c) for l 0 a = 10. This implies that the drag extrema become more difficult to detect as the flow becomes more hydrostatic. There is an additional subtle difference between Figures 6 and 7, whose cause is not clear. The drag to the left of the minimum in Figure 7(b) and to the left of the maximum in Figure 7(d) departs slightly more from 1 for l 0 a = 10 than for l 0 a = 5. Despite these differences, the behaviour of the analytical

Numerical simulations with friction
It was seen in the preceding sections that friction is a crucial effect for the type of resonance being addressed in this study, in particular for producing the single drag maximum (φ = 3π/2) or minimum (φ = π/2) at n/l 0 = 2. Since, among these two, the most relevant is undoubtedly the drag maximum, attention will be focused next on this case, with a preliminary analysis of the effect of physical friction (as opposed to numerical friction) on its behaviour. The following numerical results do not aim at more than illustrating how the drag variation is modified when a turbulence closure is adopted in the FLEX model, instead of running it in inviscid mode. A more comprehensive exploration of these effects is left for future studies. Figure 8 reproduces the results from the FLEX numerical model for φ = 3π/2, ε = 0.1 and l 0 a = 5, also presented in Figure 2(d) (squares and dashed line). The circles and solid line correspond to the drag calculated with the same parameters, with the difference that a simple Smagorinskytype turbulence closure, based on Lilly (1962), was turned on. For the triangles and dotted line, a K − ε turbulence closure was adopted instead. Clearly, by comparison with the inviscid drag, the presumably more realistic drag produced with the turbulence closures is somewhat reduced for all values of n/l 0 . For n/l 0 = 2 this reduction is larger using the Lilly closure than for the K − ε closure. The level of the nonresonant drag in the former case (≈ 0.7) is lower than that of the latter (≈ 0.8), which is comparable to that found in the viscous simulations of WV10. However, for n/l 0 = 2, the drag using both turbulence closures is comparable (≈ 1.2) and smaller than both the inviscid result and that shown in Figure 9 of WV10.
If the drag is normalized, not by its inviscid nonresonant value but by its value at the origin of n/l 0 (which corresponds to the non-resonant viscous limit) (Figure 8(b)), it can be noticed that there are almost no differences between the inviscid result and that using the Lilly turbulence closure. This means that the drag is attenuated proportionally by friction in the latter case in all circumstances. Using the K − ε turbulence closure, on the contrary, the drag amplification is somewhat reduced compared with the inviscid simulations. This behaviour is also different from that displayed in Figure 9 of WV10, where a simple turbulence closure appears to have been used. Presumably, the K − ε turbulence closure, unlike the Lilly closure, becomes more active in resonant conditions, which makes some sense, since the flow is then more likely to become turbulent. The drag behaviour using the K − ε turbulence closure could be mimicked using the present analytical model by increasing the value of λa/U, but, of course, the selected value would only be suitable for this particular case. These results further emphasize the sensitivity of the drag to the representation of frictional effects, a finding which parallels those of previous authors for various types of resonant or high-drag orographic flows, for examplé Olafsson and Bougeault (1997), Peng and Thompson (2003) and Stiperski and Grubišić (2011). Clearly, in order to achieve a realistic representation of the drag in nature, particularly for the type of resonant flow being investigated here, much attention needs to be devoted to the formulation of turbulence closures in numerical models.

Concluding remarks
A mechanism of parametric resonance leading to the amplification of mountain waves, and their associated surface drag, was investigated, inspired by the recent study of WV10. This resonance relies on the existence of a vertically oscillating Scorer parameter, although this oscillation may be of relatively small amplitude. WV10 suggested that this mechanism is intrinsically nonlinear, being related to the wave-triad interaction originally addressed by Phillips (1968) in an oceanographic context. A linear model using a perturbation approach, where the solutions are expanded in powers of a small parameter proportional to the amplitude of the Scorer parameter oscillation, was developed, showing that the wave and drag amplification mechanism under consideration is essentially linear, although it is obviously altered by nonlinearity when the drag is enhanced or attenuated by a large factor.
The results presented in the previous sections have substantial implications for numerical modelling of mountain waves. Since friction has a crucial impact on the magnitude, and even on the existence or not of drag maxima caused by the kind of parametric resonance being studied, a good representation of this process in numerical models appears to be essential to obtain accurate mountain wave drag estimates. The behaviour of the drag should be quite sensitive to the particular turbulence closure employed, as is suggested by the preliminary viscous results presented in the preceding section. Inviscid calculations, on the other hand, are probably of limited quantitative value, since their resolution, or the degree of numerical diffusion associated with their particular discretization scheme, may alter substantially the way in which the resonance process is represented.
When comparing the results of the present study to those of WV10, some intriguing aspects were found. For example, on their page 436 these authors remark that the drag maxima produced by their numerical model corresponds to drag minima in their linear model. No such discrepancy in behaviour was ever found in the present calculations.
Instead of using Eq. (7) to define the Scorer parameter, WV10 added an expression of the form A sin(nz − φ wv ) to the mean velocity U. At first glance, a simple relation between the perturbation imposed on the Scorer parameter . Normalized drag as function of n/l 0 for ε = 0.1, l 0 a = 5 and φ = 3π/2 from the FLEX numerical model, in inviscid conditions (dashed lines and squares), using Lilly's (1962) turbulence closure (solid lines and circles) and using a K − ε turbulence closure (dotted lines and triangles). (a) Same normalization as used previously. (b) Drag normalized by its value at n/l 0 = 0. by WV10 and that employed here could be described by when A/U is small. Then, the small parameter used in the present study should be defined in terms of the quantities employed by WV10 as ε = 2A/U. However, the fact that in the present analytical model the oscillation in the Scorer parameter profile was implicitly imposed on N 2 rather than on U may lead to important differences in the results (Vosper, private communication). A possible explanation for this behaviour is that some of the integrals calculated in the solution procedure involve λa/U (see section 2.1), which would become a function of z instead of a constant if U oscillates in the vertical. This would considerably complicate the analytical treatment. Several refinements and additions to the idealized situation considered here would be possible, and of substantial interest for mountain wave modelling. For example, 3D orography, which is obviously more realistic than a 2D ridge, could be adopted. This modification would be expected to weaken the resonance process, since it leads to a higher degree of wave dispersion (then not only associated with non-hydrostatic effects, but also with the variable orientation of the mountain waves).
A further step towards making the present model problem more realistic could be the prescription of Scorer parameter profiles with a more complicated form (i.e. containing more harmonics in the vertical). One of the motivations presented by WV10 for considering one single harmonic, as is done in the present study, results from their analysis of the spectrum of a more realistic Scorer parameter profile. Consideration of this effect would also, in principle, tend to weaken the resonance process under study, by spreading it over a wider range of n/l 0 than presently.
Finally, higher mountains, for which N 0 h 0 /U is closer to unity, and hence where the flow becomes nonlinear even in the absence of resonance, could be considered. Clearly, this more realistic situation could only be investigated using a set of numerical simulations. This, as well as the other developments alluded to above, are left as suggestions for future work.

Appendix A
The vertical wave number of the waves Using Eq. (13), Eq. (12) may be decomposed into its real and imaginary parts, yielding the following two equations: In order to calculate m I , Eq. (A3) must be solved for m 2 I and then the positive root for m I must be selected. Combining Eqs (A1) and (A2) in an alternative way, a relationship that gives m R in terms of m I is obtained: It can be seen from Eq. (A4) that if m I > 0 then m R has the same sign as Uk, which automatically satisfies the radiation boundary condition. The m I found as a solution to Eq. (A3) and the corresponding m R calculated from Eq. (A4) are presented in Eqs (14) and (15).