Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere

The surface drag force produced by trapped lee waves and upward propagating waves in non-hydrostatic stratiﬁed ﬂow over a mountain ridge is explicitly calculated using linear theory for a two-layer atmosphere with piecewise-constant static stability and wind speed proﬁles. The behaviour of the drag normalized by its hydrostatic single-layer reference value is investigated as a function of the ratio of the Scorer parameters in the two layers l 2 / l 1 and of the corresponding dimensionless interface height l 1 H , for selected values of the dimensionless ridge width l 1 a and ratio of wind speeds in the two layers. When l 2 / l 1 → 1 , the propagating wave drag approaches 1 in approximately hydrostatic conditions, and the trapped lee wave drag vanishes. As l 2 / l 1 decreases, the propagating wave drag progressively displays an oscillatory behaviour with l 1 H , with maxima of increasing magnitude due to constructive interference of reﬂected waves in the lower layer. The trapped lee wave drag shows localized maxima associated with each resonant trapped lee wave mode, occurring for small l 2 / l 1 and slightly higher values of l 1 H than the propagating wave drag maxima. As l 1 a decreases, i.e. the ﬂow becomes more non-hydrostatic, the propagating wave drag decreases and the regions of non-zero trapped lee wave drag extend to higher l 2 / l 1 . These results are conﬁrmed by numerical simulations for l 2 / l 1 = 0 . 2 . In parameter ranges of meteorological relevance, the trapped lee wave drag may have a magnitude comparable to that of propagating wave drag, and be larger than the reference single-layer drag. This may have implications for drag parametrization in global climate and weather-prediction models. Copyright


Introduction
Much attention has been devoted in recent years to the problem of gravity wave drag produced by stratified flow over orography, a subgrid-scale force that must be parametrized in global weather-prediction and climate models. Following pioneering contributions, where leadingorder effects of the incoming flow and orography on the drag were evaluated analytically using linear theory for an atmosphere with constant wind and static stability (Smith, 1980;Phillips, 1984), various refinements to these conditions have been considered. For example, Smith (1986) and Grubišić and Smolarkiewicz (1997) addressed the effect of a linearly varying wind profile on the drag, and Grisogono (1994) assessed boundary-layer effects, in conjunction with a more general variation of the wind with height.  and Miranda (2004, 2006) applied a WKB approximation to compute corrections to the drag due to a generic, but relatively slow, variation of the wind with height. Shutts and Gadian (1999) and, more recently, Teixeira and Miranda (2009) considered effects of directional shear on the momentum flux associated with mountain waves, which is a quantity that has a direct impact on the deceleration of the large-scale atmospheric circulation. Discontinuities in the atmospheric profiles of the wind shear and static stability, and their impacts on the drag, have been examined by Wang and Lin (1999), Leutbecher (2001) and Teixeira et al. (2005Teixeira et al. ( , 2008, among others. These authors have shown that the drag can be substantially modulated by resonance associated with partial reflections of the gravity waves at sharp variations of the atmospheric parameters, or their derivatives. Most of these investigations treated the flow as purely hydrostatic, because non-hydrostatic effects are typically relatively weak in the atmosphere at the mesoscale; it is generally believed that the main contributions to the drag arise in nearly hydrostatic conditions, and the hydrostatic approximation facilitates the calculation of closed-form analytical expressions for the drag. However, when the atmospheric parameters vary in a realistic way, in particular when the wind speed increases sufficiently, or the static stability decreases sufficiently, the flow becomes significantly non-hydrostatic at some levels. The hydrostatic approximation implies that all gravity waves in the spectrum of disturbances forced by a given orography propagate vertically. In non-hydrostatic conditions, on the other hand, there are parts of the wave spectrum (high wavenumbers) which are evanescent, and thus do not contribute to the drag, while others (low wavenumbers) propagate vertically, contributing to the drag. If an evanescent layer lies above a wave-propagating layer, wave trapping can occur, with total vertical wave reflection of some wave components at particular heights, and consequent resonant drag amplification. This phenomenon co-exists with the partial wave reflections that may occur even in hydrostatic conditions (Leutbecher, 2001;Teixeira et al., 2005).
It is known that vertically trapped internal gravity waves, generally known as trapped lee waves, produce a drag force (Bretherton, 1969;Smith, 1976;Lott, 1998;Broad, 2002), but, unlike drag in a hydrostatic atmosphere, its behaviour has not been explored in detail (Wurtele et al., 1996). Extending the calculations of Bretherton (1969) for flow confined vertically by a rigid lid (see also Tutiś, 1992), Smith (1976) presented a formula for trapped lee wave drag in an unbounded atmosphere, based on linear theory, where this quantity is expressed as the ratio of two terms, one involving an integral of the Fourier transform of the vertical velocity perturbation and the other the vertical derivative of this Fourier transform. However, Smith (1976) did not provide any explicit rigorous calculation based on this expression, and limited himself to giving a rough estimate of the trapped lee wave drag based on measured atmospheric parameters. More recently, Vosper (1996), Grubišić and Stiperski (2009) and Stiperski and Grubišić (2011) studied the trapped lee wave resonance that occurs in flow over two successive mountains. Using numerical simulations, they showed, among many other results, the variation of the drag with the horizontal separation of the mountains. For sufficiently large separations, and smooth atmospheric profiles, which preclude partial wave reflections, this variation can only be attributed to trapped lee waves. The amplitude that may thus be inferred for the trapped lee wave drag from the results of Stiperski and Grubišić (2011) is comparable to that of the remaining drag, which must be associated with vertically propagating waves. This is surprising, given that trapped lee waves, being highly non-hydrostatic, would be expected to produce relatively little drag. This result motivates the present study, where the behaviour of the drag associated with trapped lee waves and with vertically propagating waves will be investigated as a function of input parameters for much simpler atmospheric profiles than considered by Grubišić and Stiperski (2009) or Stiperski and Grubišić (2011). Scorer (1949) was the first to provide a satisfactory theoretical explanation for the existence of trapped lee waves, assuming an atmosphere with two layers, where wave propagation is permitted in the lower layer, but not allowed in the upper layer for an important fraction of the wavenumbers. This author used analytical techniques to obtain the flow configuration associated with these waves relatively far from the orography that generates them, but did not calculate the drag. In the present study, a two-layer atmosphere with piecewise-constant parameters, such as adopted by Scorer (1949), will be assumed to evaluate the trapped lee wave drag and the drag associated with vertically propagating waves, and compare the magnitude of these two drag components. It will be seen that, in some circumstances, the trapped lee wave drag may be comparable, or even larger, than the drag associated with vertically propagating waves, and substantially larger than the drag for a hydrostatic atmosphere with a constant Scorer parameter equal to that existing in the lower layer. These results have implications for gravity wave drag parametrization (Lott, 1998).
This article is organized as follows. In section 2, the linear model that will be used to calculate the trapped lee wave drag and propagating wave drag is described. Section 3 presents illustrative calculations of the drag as a function of input parameters, for three different atmospheric profiles, and a brief comparison with numerical simulations. Finally, in section 4, some concluding remarks are presented.

Theoretical model
Since a systematic exploration of trapped lee wave drag has not been carried out previously, this problem will be reduced here to a form as simple as possible. Inviscid flow over orography will be considered, and calculations will be carried out essentially using linear theory, although a few results from numerical simulations will also be presented.
Consider inviscid, stationary, non-rotating, stratified, flow over a 2D mountain ridge aligned in the y direction. If the equations of motion with the Boussinesq approximation are linearized with respect to a reference mean state, differentiated and combined, and the dependent variables are expressed as Fourier integrals along x, it can be shown that the Fourier transform of the vertical velocity perturbation, w, satisfies (Lin, 2007) where is the Scorer parameter of the atmosphere. k is the horizontal wavenumber of the waves in the x direction, N 2 (z) > 0 is the static stability of the reference state, and U(z) is the incoming wind velocity (aligned in the x direction and thus perpendicular to the orography). A two-layer atmosphere similar to that prescribed by Scorer (1949) will be assumed henceforth, with the lower layer existing between z = 0 and z = H and the upper layer extending indefinitely above. In each of the layers, both the static stability and the wind velocity are assumed to be constant, so that the corresponding Scorer parameters are also constant. However, the Scorer parameter is assumed to be discontinuous at z = H. The incoming wind velocity, static stability and Scorer parameter in the lower layer are denoted by U 1 , N 2 1 and l 1 , respectively, while in the upper layer the same quantities are U 2 , N 2 2 and l 2 . Although the situation where U 1 = U 2 would produce Kelvin-Helmholtz instability, and that aspect is not taken into account in the model, this situation will nevertheless be contemplated here, as a crude representation of the effect of the wind variation with height. Since trapped lee waves are the main focus of the present study, l 2 < l 1 will always be considered, as in Scorer (1949), which is a necessary condition for the occurrence of these waves.
In such a model set-up, three possibilities exist for waves of a given horizontal wavenumber: these waves may propagate vertically in both layers, they may propagate only in the lower layer, or they may be evanescent in both layers. The wave solutions in the lower layer, when there is propagation or when the waves are evanescent, are respectively. The wave solutions in the upper layer, in the same circumstances, are where m 1 = l 2 1 − k 2 1/2 sgn(U 1 k), n 1 = k 2 − l 2 1 1/2 , m 2 = l 2 2 − k 2 1/2 sgn(U 2 k), n 2 = k 2 − l 2 2 1/2 , and i = √ −1. In (3) the two terms correspond, respectively, to waves whose energy propagates upward and downward, while in (4) they correspond to evanescent waves that decay or amplify exponentially with height, respectively. In (5) and (6), only the upward propagating and decaying terms are considered, since the upper layer is vertically unbounded. The coefficients a 1 , b 1 , c 1 , d 1 , a 2 and c 2 are functions of k.
The boundary conditions to which these solutions are subjected result from the requirement that the flow follows the terrain elevation at the surface, and continuity of the vertical streamline displacement and pressure at the interface separating the two layers. The Fourier transform of the vertical streamline displacement η is defined by (Smith, 1979) where the index q = 1, 2 denotes values in the lower or in the upper layer, respectively. The Fourier transform of the pressure perturbation, on the other hand, is given by Teixeira et al., 2012), where ρ 0 is a reference density (assumed to be constant). In the present case, the first term on the right-hand side of (8) vanishes because dU 1 /dz = dU 2 /dz = 0. Therefore, from (7) and (8), in the linear approximation, the boundary conditions are expressed as where h is the Fourier transform of the surface elevation. If these boundary conditions are applied to (3)-(4) and (5)-(6), the coefficients in these expressions may be determined. For calculating the drag, it is sufficient to know a 1 , b 1 , c 1 and d 1 . These coefficients, which determine the flow in the lower layer, are provided in Appendix A.
Following Wang and Lin (1999) and Teixeira et al. (2005), a wave reflection coefficient at z = H may be defined as the modulus of the ratio between the downward propagating and the upward propagating wave components in the lower layer, R = |b 1 /a 1 |. For trapped lee waves (l 2 2 < k 2 < l 2 1 ), it may be checked from the definitions of a 1 and b 1 in Appendix A that R = 1, as would be expected. For vertically propagating waves (k 2 < l 2 2 ), the reflection is only partial and R = U 2 1 l 2 1 − k 2 1/2 − U 2 2 l 2 2 − k 2 1/2 U 2 1 l 2 1 − k 2 1/2 + U 2 2 l 2 2 − k 2 1/2 , which satisfies R ≤ 1. Unfortunately, R depends in general on k, so its relation to the behaviour of propagating waves is not trivial, as will be seen. The gravity wave drag exerted on the mountain ridge per unit spanwise length is given by ) where the asterisk denotes complex conjugate. The second equality used the fact that both p 1 (the pressure perturbation in the lower layer) and h (the surface elevation) are real, and the properties of Fourier transforms. If (8) is evaluated at z = 0 and (3)-(4) are taken into account, it follows that when the waves are propagating in the lower layer, whereas when the waves are evanescent in this layer, then From these two equations and (A1)-(A6), it is possible to obtain explicit expressions for p(z = 0). These expressions may then be used in (13) to obtain the internal gravity wave drag. From the partition of the wave solutions in wavenumber ranges, according to (3)-(6), it turns out that the drag comprises three terms: The first term, D 1 , receives contributions from wavenumbers |k| < l 2 , corresponding to waves that propagate in both atmospheric layers. This will be called 'propagating wave drag'. The second term, D 2 , receives contributions from wavenumbers satisfying l 2 < |k| < l 1 , and corresponds to trapped lee wave drag. The third term, D 3 , receives contributions from wavenumbers |k| > l 1 , corresponding to evanescent waves. Using all the previously derived results, the first term may be written which was obtained by taking the imaginary part of the corresponding integral in (13). In general, the integral in (17) must be evaluated numerically. This drag term is the only one that exists in hydrostatic flow, and its behaviour for a variation of l only due to N 2 was studied, for example, by Leutbecher (2001). The third drag term can be shown to have the form: k| h| 2 n 1 U 2 2 n 2 cosh(n 1 H)+U 2 1 n 1 sinh(n 1 H) U 2 1 n 1 cosh(n 1 H)+U 2 2 n 2 sinh(n 1 H) dk .
Since the only contributions to this term come from the imaginary part of the integral, but the integrand is real, and without singularities on the real axis, it follows that The second drag term is the main focus of the present analysis. Following the same procedure as for the previous terms, it may be expressed as While, as in the previous case, the contributions to the drag only come from the imaginary part of the integral, now these contributions actually exist in some circumstances, due to the singularities of the integrand along the real axis. These singularities are given by which results from equating the denominator of the fraction in (20) to zero. Because of the periodicity and range of variation of the tan function, and taking into account the definition of m 1 , it may be shown that the number of singularities, i.e. the number of trapped lee wave modes N L , is related to other flow parameters by This implies, in particular, that there will only be at least one singularity when as originally stated by Scorer (1949). Note that none of these conditions explicitly depends on U 1 or U 2 . Provided (23) is satisfied, (20) is best evaluated through a change of variable, as described in Appendix B. Going back to the original variable k, (B5) obtained in Appendix B can be expressed as where the sum is carried out over all trapped lee wave modes and k j are the corresponding resonant wavenumbers, given by (B1) with φ = φ j , or alternatively directly from (21). Equation (24) constitutes perhaps the most important result of the present study and, in conjunction with (17), will be used next to compare the propagating wave drag and the trapped lee wave drag. When U 1 = U 2 = U, (24) simplifies considerably, reducing to Equation (25) could be confirmed by an alternative method, as a particular case of the formula derived by Smith (1976) for trapped lee wave drag, when the corresponding wave solutions are substituted in the integral and derivative contained in this formula, and the integral is calculated analytically. This provides a good consistency check for the present calculations.

Results
The behaviour of propagating wave drag and trapped lee wave drag will now be investigated for a bell-shaped ridge given by where a is the half-width of the ridge and h 0 is its maximum height. This kind of orography is only used for illustrative purposes, and the results to be presented should remain essentially valid for a different type of orography. D 1 and D 2 will be normalized by the drag for a hydrostatic atmosphere with a constant Scorer parameter equal to that existing in the lower layer and the same type of orography. This reference drag is given by (Smith, 1979) From (17), (24) and (28), the normalized propagating wave drag and trapped lee wave drag may be expressed, respectively, as and where, in (30), the summation is carried out over all trapped lee wave modes. D 1 /D 0 and D 2 /D 0 are clearly functions of l 1 H, l 2 H, a/H and U 1 /U 2 (k j H itself is a function of U 1 /U 2 , l 1 H and l 2 H; see (21)). Any alternative combination of these parameters that yields four independent input parameters can be considered. Here the key parameters will be taken as: l 1 H, l 1 a, l 2 /l 1 and U 1 /U 2 , since they allow an easier exploration of parameter space. l 1 a quantifies non-hydrostatic effects in the lower layer, and is adequate to distinguish the present results from previous ones calculated using the hydrostatic approximation (e.g. Teixeira et al., 2005). Three values of this parameter will be considered next, similar to those employed by Teixeira et al. (2012). l 1 H is a dimensionless height of the interface separating the two atmospheric layers, essentially similar to that defined by Miranda and Valente (1997) and Teixeira et al. (2005). l 2 /l 1 quantifies the jump in Scorer parameter at the interface, and in the present context may be varied between 0 and 1. Finally, for the jump of the wind velocity at the interface, U 1 /U 2 , three cases will be considered, to be described next in turn.
In the hydrostatic approximation, i.e. when l 1 a, l 2 a → ∞, (29) reduces to and obviously D 2 /D 0 = 0, because wave trapping becomes impossible. Unlike the previous non-hydrostatic results, (31) is independent of the form of the orography (27), in the same way as the hydrostatic results of , or Teixeira et al. (2005Teixeira et al. ( , 2008. According to (31), D 1 /D 0 has an oscillating behaviour with l 1 H, with a period of π in terms of this parameter. This is clearly caused by resonant reinforcement or weakening of the waves in the lower layer, depending on the phase of the downward propagating waves reflected at the interface z = H. D 1 /D 0 can only be 1 if either l 1 = l 2 and U 1 = U 2 (which is fairly intuitive), or when U 1 /U 2 = (l 2 /l 1 ) 1/2 , which is a non-trivial result.
In hydrostatic conditions, the reflection coefficient (12) reduces to which helps to explain the drag behaviour described above, because if U 1 /U 2 = (l 2 /l 1 ) 1/2 , then R = 0, i.e. there is no reflection of vertically propagating waves at z = H.

Discontinuity in N
If U 1 /U 2 = 1, U is continuous at the interface separating the two atmospheric layers, z = H, and the discontinuity of l is purely due to the variation of N, namely l 2 /l 1 = N 2 /N 1 . This was the case originally considered by Scorer (1949). The normalized drag will be presented first as a function of l 1 H and l 2 /l 1 for l 1 a = 10, l 1 a = 5 and l 1 a = 2. These three cases correspond to weak, moderate and strong non-hydrostatic effects. As in Leutbecher (2001) and Teixeira et al. (2005), l 1 H will be varied between 0 and 3π . Figure 1 presents the number of trapped lee wave modes, N L , as a function of l 1 H and l 2 /l 1 , calculated from (22). The behaviour of this quantity is independent both of U 1 /U 2 and of l 1 a. It can be seen that there are no trapped lee wave modes when l 2 /l 1 = 1, because wave trapping requires a decrease of l with height. So the number of trapped lee wave modes decreases as l 2 /l 1 varies between 0 and 1. On the other hand, the number of trapped lee wave modes increases with l 1 H, reaching a value of 3 at the highest values of l 1 H considered. There are no trapped lee wave modes for l 1 H/π < 1/2 because of condition (23). Figure 2 shows the (horizontal) wavelengths of the first three trapped lee wave modes (Figure 2(a, b, c)), defined as λ j = 2π/k j (with j = 1, 2, 3), normalized by the vertical wavelength of hydrostatic waves in the upper layer, 2π/l 2 ,  as a function of l 1 H and l 2 /l 1 . The normalized wavelengths depend on U 1 /U 2 , as can be inferred from (21), but not on l 1 a. Since trapped lee waves receive contributions from wavenumbers between l 2 and l 1 , λ j l 2 /(2π ) vary between l 2 /l 1 and 1, being increasing functions of l 2 /l 1 and decreasing functions of l 1 H. Figures 3, 4 and 5 present the drag associated with mountain waves (for l 1 a = 10, l 1 a = 5 and l 1 a = 2, respectively) normalized by the drag D 0 for hydrostatic flow in a single-layer atmosphere with Scorer parameter l 1 , given by (28). Figures 3(a), 4(a) and 5(a) show the drag associated with propagating mountain waves, Figures 3(b), 4(b) and 5(b) the drag associated with trapped lee waves, and Figures 3(c), 4(c) and 5(c) the total drag, which is a sum of the two. Unlike the quantities presented in Figures 1 and  2, the normalized drag depends on l 1 a in addition to U 1 /U 2 .
In Figures 3(a), and 4(a) it can be seen that the propagating wave drag D 1 /D 0 approaches values ≈ 1 when l 2 /l 1 → 1, as would be expected for reasonably hydrostatic flow, but in Figure 5(a), these values become noticeably lower because of non-hydrostatic effects. As l 2 /l 1 decreases towards zero, D 1 /D 0 increasingly shows a periodic oscillatory behaviour with l 1 H, with maxima at l 1 H/π = 0.5 + n, where n is an integer, and minima at l 1 H/π = n. The location of the maxima and minima may be explained by (31) with U 1 /U 2 = 1. This behaviour also resembles that shown in Figure 2 of Leutbecher (2001), which assumes the hydrostatic approximation, with the difference that the maxima and minima are exchanged, because Leutbecher considered a case where l 1 < l 2 . These maxima and minima are obviously produced by resonance associated with partial wave reflections at z = H, and make the normalized drag at low l 2 /l 1 take values between D 1 /D 0 ≈ 0 and D/D 0 ≈ 9.3 in Figure 3(a), and D/D 0 ≈ 4.6 in Figure 4(a) (occurring for slightly higher l 2 /l 1 ). In Figure 5(a), the drag maxima are considerably weakened (not exceeding ≈ 1.8), and become distinctly smaller as l 1 H increases. This behaviour, which is also visible in Figure 4(a) (but to a lesser extent), is due to decrease of the wave energy density in the lower atmospheric layer (when this energy propagates downstream, as happens in non-hydrostatic conditions) as the top of this layer is lifted. In Figures 4(a) and 5(a) there are gaps in the D 1 /D 0 field near the right edges of the maxima.
In Figures 3(b), 4(b) and 5(b), it can be seen that the behaviour of the trapped lee wave drag D 2 /D 0 is substantially different. The drag is zero except in localized regions of parameter space, which, of course, are contained in the wider regions where at least one trapped lee wave mode is allowed to exist. In Figure 3(b), there is only appreciable drag for l 2 /l 1 < 0.5 and inside intervals of l 1 H with their lower limits coinciding with the locations of the D 1 /D 0 maxima in Figure 3(a). The width of these intervals increases with l 1 H. These intervals are much less localized in parameter space in Figure 4(b) (for example, the non-zero drag regions surrounding the second and third maxima overlap), and occupy a yet much larger fraction of parameter space in Figure 5(b), with the three displayed maxima having merged together. In Figure 3(b), the magnitude of the first maximum is larger than those of D 1 /D 0 (D 2 /D 0 ≈ 20.6), but the maxima are attained at l 2 /l 1 = 0, requiring an extreme contrast in Scorer parameter between the two layers. In Figure 4(b) the maximum of the trapped lee wave drag is D 2 /D 0 ≈ 9.9, but values of O(1) are attained at higher l 2 /l 1 than in Figure 3(b), which is more meteorologically relevant. In Figure 5(b), D 2 /D 0 reaches a value of ≈ 3.2, which is still larger than the maximum of the propagating wave drag. The trapped lee wave drag takes values of O(1) for l 2 /l 1 as high as 0.4, which is not hard to satisfy in approximations to a realistic atmosphere.
In Figures 3(c) and 4(c), the total drag D/D 0 considerably resembles the propagating wave drag D 1 /D 0 , for l 2 /l 1 ≈ 1 but is dominated by the trapped lee wave drag for l 2 /l 1 ≈ 0. In Figure 5    The results shown in Figures 3-5 clearly stress the importance of trapped lee wave drag, showing that it may become an important fraction of the total drag, and be comparable to the hydrostatic one-layer drag, especially in considerably non-hydrostatic conditions. These aspects are examined in more detail in Figure 6 for l 2 /l 1 = 0.2, a value chosen for purely illustrative purposes. This may seem a low value, but a large decrease of the Scorer parameter with height (although not necessarily discontinuous, and not necessarily due to the variation of N) is possible in real atmospheres (Grubišić and Stiperski, 2009;Stiperski and Grubišić, 2011). Figure 6(a) shows the number of trapped lee wave modes, N L , as a function of l 1 H/π . There is one trapped lee wave mode between l 1 H/π ≈ 0.5 and l 1 H/π ≈ 1.5, two trapped lee wave modes between l 1 H/π ≈ 1.5 and l 1 H/π ≈ 2.5, and three trapped lee wave modes above l 1 H/π ≈ 2.5. As pointed out for Figure 1, this result is independent both of U 1 /U 2 and of l 1 a. In Figure 6(b), the normalized trapped lee wave wavelengths, λ j l 2 /(2π ), are shown as a function of l 1 H/π . It can be seen that these quantities always depart from a value of 1 when a new trapped lee wave mode is established, but as l 1 H/π increases they asymptotically tend to l 2 /l 1 , which in the present case is 0.2. Figure 6(c) shows the propagating wave drag as a function of l 1 H/π for the three values of l 1 a considered previously, namely l 1 a = 2, 5, 10. As noted in Figures 3-5, D 1 /D 0 displays maxima at about l 1 H/π = 0.5 + n, separated by minima at about l 1 H/π = n, where n is an integer. It can be shown from (31) that, in hydrostatic conditions, all maxima and minima have magnitudes of l 1 /l 2 and l 2 /l 1 , respectively (cf. Leutbecher, 2001). However, non-hydrostatic effects decrease these maxima and make them become smaller as l 1 H/π increases, for reasons pointed out earlier, related to downstream wave propagation. It can also be shown that non-hydrostatic effects, when relatively weak (as in the case with l 1 a = 10), can make the drag in the first (and even in the second) maximum exceed the hydrostatic limit, as is observed in Figure 6(c). Maxima of D 1 /D 0 in Figure 6(c) reach a value of ≈ 5.7 for l 1 a = 10, of ≈ 4.5 for l 1 a = 5 and of ≈ 1.7 for l 1 a = 2. The drag maxima for l 1 a = 5 and especially l 1 a = 2 are sharply peaked.
Figure 6(d) shows the trapped lee wave drag as a function of l 1 H/π . As noted in Figures 3-5, the maxima of D 2 /D 0 are more localized than those of D 1 /D 0 and are displaced to higher values of l 1 H/π . The drag also increases suddenly to its maximum as l 1 H/π increases (as a new lee wave mode is established), then decreasing much more gradually. The first maximum of D 2 /D 0 reaches ≈ 1.2 for l 1 a = 10, ≈ 2.9 for l 1 a = 5 and ≈ 2.2 for l 1 a = 2. The increase of these maxima from l 1 a = 10 to l 1 a = 5 may be attributed to the spread in parameter space of the region where the trapped lee wave drag is non-zero (see comments on Figures 3-5). The decrease from l 1 a = 5 to l 1 a = 2 may be attributed to weakening of the trapped lee waves by downstream propagation, as explained earlier -an effect that also explains why the second trapped lee wave maximum is lower than the first, and the third is lower than the second. Physically, the decrease of trapped lee wave drag at high l 1 a probably occurs because the wavelength of the trapped lee waves becomes too short compared to the width of the mountain, so the contributions of the corresponding positive and negative pressure anomalies to the drag tend to cancel out.
The total drag D/D 0 is shown as a function of l 1 H/π in Figure 6(e). The behaviour of D/D 0 resembles that of D 1 /D 0 , especially for l 1 a = 10, but departs from it progressively as l 1 a decreases. Maxima of D/D 0 reach ≈ 5.7 for l 1 a = 10, ≈ 5.0 for l 1 a = 5 and ≈ 2.5 for l 1 a = 2. More importantly, gaps on the right slopes of the maxima of D 1 /D 0 are filled in D/D 0 , making the corresponding peaks become wider and smoother, and extend to higher values of l 1 H/π . The trapped lee wave drag and the propagating wave drag are thus confirmed to have a complementary role, owing to their imperfect overlap in parameter space.
A comparison of some of these results with numerical simulations will be carried out in section 3.4.

Discontinuity in U
A situation where the discontinuity of l is purely due to the variation of U is now considered. N is assumed to be continuous at z = H (i.e. N 1 = N 2 ), so that l 2 /l 1 = U 1 /U 2 . This situation serves to illustrate leading-order effects of the wind variation with height on the drag using the same theoretical framework.
Because the number of trapped lee wave modes does not depend on U 1 /U 2 , it is unnecessary to repeat a graph of this quantity here, and Figure 1 remains valid. Figure 7 presents the wavelengths of the first three trapped lee wave modes, λ j l 2 /(2π ) (j = 1, 2, 3), normalized in the same way as previously. The behaviour of these quantities bears some resemblance to that shown in Figure 2, however, one important difference is that, instead of decreasing steadily with l 1 H/π , the normalized wavelengths take an approximately constant value of ≈ 1 from the l 1 H/π where the lee wave mode is established up to a l 1 H/π higher by about 0.5. (For relatively low l 2 /l 1 , this corresponds to l 1 H/π ≈ 1, 2, 3, respectively, for the first, second and third trapped lee wave modes.) Because this effect is most pronounced at l 2 /l 1 = 0, the wavelength is no longer a uniformly increasing function of l 2 /l 1 .
In Figures 8 and 9, the two components of the drag and their sum are presented for l 1 a = 10 and l 1 a = 2, respectively (the case l 1 a = 5 is omitted for brevity). In Figures 8(a) and 9(a), the main differences from Figures 3(a) and 5(a) are that the propagating wave drag D 1 /D 0 has maxima at l 1 H/π = n (where n is an integer) instead of at l 1 H/π = 0.5 + n, and these maxima are lower, never exceeding ≈ 4.6 and ≈ 0.9, respectively. The location of these maxima is easily explained using (31), by assuming that U 1 /U 2 = l 2 /l 1 . The trapped lee wave drag, shown in Figures 8(b) and 9(b), on the other hand, has non-zero values in approximately the same regions as in Figures 3(b) and 4(b), but the maxima occur at values of l 1 H/π slightly higher than l 1 H/π = n, rather than higher than l 1 H/π = 0.5 + n. Additionally, the maxima of D 2 /D 0 are much larger, attaining a value ≈ 179 for l 2 /l 1 = 0 in Figure 8(b) and ≈ 9.8 in Figure 9(b). The total drag, displayed in Figures 8(c) and 9(c), is dominated by these maxima, reaching values of ≈ 180 and ≈ 9.8 respectively (which are of limited meteorological relevance), but at higher l 2 /l 1 is dominated by the propagating wave drag contribution (particularly in Figure 8(c)). As in Figure 5, Figure 9(c) has a considerably smoother appearance than D 1 /D 0 in Figure 9(a), due to filling of some gaps in the D 1 /D 0 field by the trapped lee wave drag. Figure 10 presents results for l 2 /l 1 = 0.2, as in Figure 6, but where the discontinuity in l is purely due to U, namely U 1 /U 2 = 0.2. Figure 10(a) shows the normalized wavelengths of trapped lee waves λ j l 2 /(2π ) as a function of l 1 H/π . The wavelengths are decreasing functions of l 1 H/π , with the same asymptotic values as in Figure 6(a), but show a plateau after l 1 H/π = 0.5 + n, as noted in Figure 7. As pointed out in Figure 8, the maxima of the propagating wave drag D 1 /D 0 , shown in Figure 10(b), are now centred at l 1 H/π = n, where n is an integer. The magnitude of these maxima is ≈ 3.9 for l 1 a = 10, ≈ 2.3 for l 1 a = 5 and ≈ 0.7 for l 1 a = 2. Curiously, the drag maxima decrease very little as l 1 H/π increases, which means that the effect of downstream propagation on the decrease of the energy density in the lower layer seems relatively modest here. However, it probably only makes sense to analyze this effect for the total drag. The trapped lee wave drag, displayed in Figure 10(c), has a different configuration from the corresponding quantity displayed in Figure 6(d). Instead of extending primarily to the right of the maxima, as in Figure 6(d), the non-zero drag regions can now extend to both sides, since the condition for the existence of trapped lee waves remains the same as in section 3.1. The maxima of D 2 /D 0 take a modest value of ≈ 0.6 for l 1 a = 10, ≈ 2.2 for l 1 a = 5 and ≈ 2.3 for l 1 a = 2. The total drag, displayed in Figure 10(d), attains maxima of ≈ 4.3 for l 1 a = 10, of ≈ 4.4 for l 1 a = 5 and of ≈ 2.9 for l 1 a = 2.
Differences between the results presented in this section and those presented in section 3.1 may be attributed primarily to the effect of the boundary conditions at z = H, (10) and (11), since these depend on the ratio U 1 /U 2 , unlike the Taylor-Goldstein equation (1), which only depends on the Scorer parameter and not on N and U individually. Clearly, the phase of the wave reflection at z = H is altered by the jump in U, so that locations of the drag maxima are shifted, in agreement with (31). Concerning the drag, a unifying feature between Figures 6 and 10 is that drag maxima occur in regions of parameter space where the trapped lee wave wavenumber varies rapidly with l 1 H/π . This rapid variation results directly from the resonance condition (21), but the fact it coincides with drag maxima is a consequence of D 2 being proportional to m 2 1 and n 2 (24), both of which become zero for the extreme values of λ j l 2 /(2π ). Hence the drag is maximized for intermediate values of this quantity.
From a theoretical point of view, it is instructive to consider an intermediate situation, where both N and U are discontinuous at z = H. This is done next.

Discontinuities both in N and U
When the discontinuity in l at z = H is due simultaneously to discontinuities in N and U, many possible combinations exist. Since the situations treated in the previous two sections led to a different location of the drag maxima, a transition between these two cases is aimed at. Equations (31) and (32) suggest that this transition occurs for U 1 /U 2 = (l 2 /l 1 ) 1/2 , where, as pointed out earlier, D 1 /D 0 = 1 in the hydrostatic approximation. However, non-hydrostatic effects slightly modify this situation. For that reason, U 1 /U 2 = (l 2 /l 1 ) 0.6 , implying N 2 /N 1 = (l 2 /l 1 ) 0.4 , is assumed here instead, which better focuses on the transition. (The value of 0.6 for the exponent was found by trial and error.) Concerning the number of existing trapped lee wave modes, Figure 1 again remains valid, because this quantity does not depend on U 1 /U 2 . Figure 11 presents the normalized wavelengths of the trapped lee wave modes, λ j l 2 /(2π ) (j = 1, 2, 3). The behaviour of these quantities is qualitatively similar and  Figure  7, extending for l 2 /l 1 = 0 mostly up to l 1 H/π ≈ 0.75, 1.75 and 2.75 for the first, second and third trapped lee wave modes, respectively.
In Figures 12 and 13, the drag behaviour is displayed for l 1 a = 10 and l 1 a = 2, respectively. Figure 12(a) shows the propagating wave drag. For l 2 /l 1 ≈ 1, D 1 /D 0 approaches 1 in Figure 12(a) as expected, but for low l 2 /l 1 , and unlike in Figures 3-5 and 8-9, there are not very well-defined drag maxima and minima, despite the large contrast in l between the two atmospheric layers. On the contrary, D 1 /D 0 tends to decrease towards zero as l 2 /l 1 → 0, and does not exceed ≈ 1.1. This behaviour, excepting the decrease of the drag to zero at l 2 /l 1 = 0, is consistent with (31), where D 1 /D 0 = 1 for U 1 /U 2 = (l 2 /l 1 ) 1/2 . Clearly, for fairly similar values of U 1 /U 2 , a relative suppression of the drag extrema is preserved, even in non-hydrostatic conditions. In Figure  13(a), D 1 /D 0 is always lower than 1 and has a maximum of ≈ 0.8, generally decreasing as l 2 /l 1 decreases, without any clearly defined peaks.
In Figures 12(b) and 13(b), where the trapped lee wave drag is shown for l 1 a = 10 and l 1 a = 2, respectively, it can be seen that the D 2 /D 0 pattern, on the other hand, resembles somewhat that shown in Figures 8(b) and 9(b), with absolute maxima of D 2 /D 0 (reaching ≈ 45.7 in Figure 12(b) and ≈ 7.6 in Figure 13(b)) centred at l 2 /l 1 = 0 and near l 1 H/π = n, where n is an integer. However, the non-zero drag zones extend to lower values of l 1 H/π than in Figures 8  and 9. More importantly, the structure of the drag maxima in (l 2 /l 1 , l 1 H/π ) parameter space is slanted, with the maxima as a function of l 1 H/π shifting to lower values of this parameter as l 2 /l 1 increases. Figures 12(c) and 13(c) present the total drag D/D 0 . The corresponding field has maxima of ≈ 45.7 in Figure 12(c) and of ≈ 7.6 in Figure 13(c), like D 2 /D 0 and, again displays a smoother structure resulting from the merging of complementary contributions given by D 1 /D 0 and D 2 /D 0 .
The main conclusion from Figures 12-13 is that the propagating wave drag and the trapped lee wave drag are affected differently by the profiles of N and U prescribed here. Specifically, no cancellation of the drag maxima like the one predicted by (31) for D 1 /D 0 is observed for D 2 /D 0 .
In Figure 14, results are presented as a function of l 1 H/π for l 2 /l 1 = 0.2, as in the two previous sections, but keeping the same proportions for the discontinuities of N and U as in Figures 12-13. Hence U 1 /U 2 = 0.2 0.6 ≈ 0.38 and N 2 /N 1 = 0.2 0.4 ≈ 0.53. Figure 14(a) confirms the comment made about Figure 11, that the decrease of the normalized wavelengths of the trapped lee waves λ j l 2 /(2π ) from their maximum of 1 to their asymptotic value of 0.2 occurs initially at a rate intermediate between those of Figures 6 and 10.
The propagating wave drag D 1 /D 0 , shown in Figure 14(b), has a very different behaviour from that shown in Figures 6(c) or 10(b), since the maxima are much lower and the minima higher, which is consistent with Figures 12-13. For l 1 a = 10, drag maxima (of ≈ 1.1) occur at l 1 H/π = n, where n is an integer, however with small secondary maxima at l 1 H/π = 0.5 + n. For l 1 a = 5, the magnitude of D 1 /D 0 generally decreases, but the maxima (of ≈ 0.7) at l 1 H/π = 0.5 + n and those at l 1 H/π = n have almost the same magnitude, although the former have sharper peaks. For l 1 a = 2, the magnitude of D 1 /D 0 decreases further, but the maxima (of ≈ 0.2) at l 1 H/π = 0.5 + n become dominant. The trapped lee wave drag D 2 /D 0 , presented in Figure 14(c), has maxima which decrease in magnitude as l 1 H/π increases and become larger as l 1 a decreases, taking values of up to ≈ 0.3 for l 1 a = 10, ≈ 1.2 for l 1 a = 5 and ≈ 1.4 for l 1 a = 2. Unlike those in Figures 6(d) and 10(c), however, the values of l 1 H/π at which these maxima are centred shift from ≈ 0.75 + n for l 1 a = 10 to slightly lower than l 1 H/π = n for l 1 a = 2. This behaviour shows that the location of the trapped lee wave drag maxima in this  transitional case is not only a function of l 2 /l 1 (as shown in Figures 12-13), but also of l 1 a.
The total drag D/D 0 , presented in Figure 14(d), shows non-zero values everywhere like D 1 /D 0 , and well-defined maxima like D 2 /D 0 , but with a smooth variation, without any sharp peaks. The maxima, which are centred between slightly above l 1 H/π ≈ 0.75 + n (for l 1 a = 10) and slightly below l 1 H/π = n (for l 1 a = 2), take values reaching ≈ 1.2 for l 1 a = 10, ≈ 1.8 for l 1 a = 5 and ≈ 1.6 for l 1 a = 2.
The behaviour of the drag in the three preceding sections can be better understood using the reflection coefficient R. Figure 15 shows R as a function of l 2 /l 1 and the wavenumber normalized by l 2 , k/l 2 , the two dimensionless parameters on which this quantity depends. Note that k/l 2 takes values between 0 and 1 for the vertically propagating waves. When U 1 /U 2 = 1 (Figure 15(a)) or U 1 /U 2 = l 2 /l 1 (Figure 15(b)) (the cases treated in sections 3.1 and 3.2, respectively), it can be seen that R varies from 0 to 1 as l 2 /l 1 varies from 1 to 0 in hydrostatic conditions (k/l 2 ≈ 0). The same qualitative tendency is also roughly followed for higher k/l 2 (i.e. non-hydrostatic conditions), but the behaviour of R becomes asymmetric, increasing with k/l 2 for U 1 /U 2 = 1 and decreasing for U 1 /U 2 = l 2 /l 1 . This may account for the relatively smaller maxima of D 1 /D 0 seen in section 3.2 for the second case, compared with those seen in section 3.1 for the first case. Figure 15(c) presents R for the case U 1 /U 2 = (l 2 /l 1 ) 1/2 . As noted previously, R = 0 in hydrostatic conditions (k/l 2 ≈ 0), but becomes R > 0 when k/l 2 > 0. In Figure 15(d), R is presented for U 1 /U 2 = (l 2 /l 1 ) 0.6 , as assumed in the present section. Although in Figure 15(d) R is non-zero for k/l 2 = 0, it attains smaller values than in Figure 15(c) for higher k/l 2 . This may account for the near-absence of D 1 /D 0 maxima in Figures 12, 13 and 14, since for U 1 /U 2 = (l 2 /l 1 ) 0.6 the reflection coefficient is presumably minimized for the range of k contributing to this part of the drag.

Comparison with numerical simulations
In order to test the preceding theoretical predictions from linear theory, some numerical simulations were carried out using the FLEX numerical model (Arga´ın et al., 2009) for the case treated in section 3.1, namely a discontinuity in l at z = H due only to the discontinuity of N. It would not be appropriate to address the cases of sections 3.2 and 3.3, where U is discontinuous, because, as mentioned before, this would produce Kelvin-Helmholtz instability at z = H, which is incompatible with steady flow.
FLEX is a 2D micro-to-mesoscale model using orthogonal curvilinear coordinates, employed previously by the authors to address various problems of resonant mountain wave flow (Teixeira et al., 2005(Teixeira et al., , 2008(Teixeira et al., , 2012. The numerical simulations performed here used a horizontal resolution of x = 200 m and a vertical resolution of z = 24 m. The wind velocity was assumed to be U = 10 m s −1 , and the static stability in the lower and upper layers was assumed to be N 1 = 0.02 s −1 and N 2 = 0.004 s −1 , respectively. This gives N 2 /N 1 = 0.2 and l 2 /l 1 = 0.2. The mountain was assumed to have a bell-shaped form, (27), with height h 0 = 10 m, which gives l 1 h 0 = 0.02, i.e. strongly linear flow. The mountain width was assumed to be a = 1000 m, 2500 m, or 5000 m. This gives, respectively, l 1 a = 2, l 1 a = 5, or l 1 a = 10, as considered in Figure 6. The number of grid points in the vertical was always N z = 916, which makes the domain extend up to z = 21984 m, or approximately 7 hydrostatic vertical wavelengths 2π/l 1 , while the number of grid points in the horizontal was N x = 250 for l 1 a = 2, N x = 625 for l 1 a = 5, and N x = 1250 for l 1 a = 10. This gives a horizontal domain extent of 50 km, 125 km, and 250 km, respectively, which always corresponds to 50a (10a upstream of the mountain and 40a downstream of it).
The time step used in the model integration was t = 0.5 s for l 1 a = 2, t = 1.0 s for l 1 a = 5, and t = 2.0 s for l 1 a = 10. The number of time steps considered was that necessary to make the drag stabilize to an approximately constant value. It corresponded to an integration time of ≈ 500a/U, or ≈ 14 h for l 1 a = 2, ≈ 35 h for l 1 a = 5, and ≈ 69 h for l 1 a = 10. A sponge layer was applied at the top of the domain over a depth of ≈ 2.5 hydrostatic wavelengths. The radiation boundary condition of Raymond and Kuo (1984) was imposed at the downstream boundary, but sponge layers, spanning 15 and 30 grid points, respectively, were also applied at the upstream and downstream boundaries. Figure 16 shows the normalized wavelengths of the trapped lee waves λ j l 2 /(2π ) from numerical simulations with l 1 a = 2 (symbols), and from linear theory (the lines reproduce those of Figure 6(b)). Numerical simulations with l 1 a = 5 and l 1 a = 10 (not shown) gave essentially similar results, which is consistent with the prediction from linear theory that λ j l 2 /(2π ) does not depend on l 1 a. The wavelengths were determined in the numerical simulations by taking the power spectrum of the vertical velocity field downstream of the mountain at a given height (z = 0.7H) and determining the maximum of this spectrum. It may be seen that there is excellent agreement between the numerical and the analytical results, with the trapped lee wave wavelength in the numerical simulations being predicted very accurately by linear theory. Normalized trapped lee wave wavelengths as a function of l 1 H/π for l 2 /l 1 = 0.2 and U 1 /U 2 = 1 from linear theory (lines) and from numerical simulations (symbols). The solid line and circles denote the first trapped lee wave mode, the dashed line and squares the second trapped lee wave mode, and the dotted line and triangles the third trapped lee wave mode. Numerical simulation results are for l 1 a = 2, although the displayed quantities should be independent of this parameter. Figure 17 shows the normalized total drag D/D 0 from the numerical simulations (symbols) and from linear theory (solid lines), for l 1 a = 10 (Figure 17(a)), for l 1 a = 5 (Figure 17(b)) and for l 1 a = 2 (Figure 17(c)). Also shown for reference are the propagating wave drag (dashed lines) and the trapped lee wave drag (dotted lines) from linear theory (although possible, it is not easy to separate these two components of the drag in the numerical results). It can be seen that the total drag given by linear theory in Figure 17(a, b) is in very good agreement with that provided by the numerical simulations. In Figure 17(a), the total wave drag is not too different from the propagating wave drag, since the trapped lee wave drag is relatively small. However, in Figure 17(b), where the magnitude of the two drag components is more comparable, it becomes clear that both their contributions are necessary to obtain good agreement with the numerical results, since none of them in isolation accurately predicts those results. This can also be noticed in Figure 17(c), where the trapped lee wave drag component has become dominant. Although the agreement is not as good as in Figure 17(a, b), clearly both drag components are necessary to correctly predict the drag given by the numerical simulations.
As a brief illustration of the flow structure associated with the drag behaviour described above, Figure 18 presents the normalized vertical velocity perturbation w/(Uh 0 /a) obtained from numerical simulations essentially similar to those used to produce the results of Figures 16 and 17. Figure 18(a) shows results for l 1 H/π = 0.5, a case just before the establishment of the first trapped lee wave mode, and so where linear theory predicts D 2 /D 1 = 0 and D/D 0 = 1.4. In Figure 18(b), a case with l 1 H/π = 0.7 is presented, where linear theory predicts D 1 /D 2 = 0.05 and D/D 0 = 1.6. In Figure 18(a), the typical structure of vertically propagating mountain waves can be seen, with the flow tilted in the upstream direction and no substantial wave decay in the upper layer. By contrast, in Figure 18(b) the structure typical of trapped lee waves is visible, with a shorter horizontal wavelength, no flow tilting in the vertical, and wave decay in the upper layer. In each case, the flow structure is clearly consistent with the drag partition associated with it.
These comparisons give some confidence in the theoretical results calculated earlier, suggesting that the cases addressed in sections 3.2 and 3.3 are also correct.

Concluding remarks
The surface drag force associated with trapped lee waves and upward propagating internal gravity waves generated in nonhydrostatic stratified flow over a 2D ridge was investigated using linear theory. Following the original approach of Scorer (1949), a two-layer atmosphere was considered, with a piecewise-constant Scorer parameter. This is the simplest possible model set-up that is able to produce trapped lee waves. However, unlike in Scorer (1949), the possibility of discontinuities in both the static stability and the wind speed was considered. The main aim of this study was to compare the magnitude of trapped lee wave drag to the magnitude of propagating wave drag, and the magnitude of those two drags to that of the drag for a single-layer hydrostatic atmosphere with a constant Scorer parameter. The very idealized conditions adopted here seem justified, since they preserve essential physical processes, while allowing the problem to be analytically tractable.
It was found that: • As initially suggested by the numerical study of Stiperski and Grubišić (2011), the trapped lee wave drag may have a magnitude comparable to that of the propagating wave drag, and both drags may be substantially larger than the reference one-layer drag. • Both drag components are maximized for large contrasts in Scorer parameter between the two layers (small l 2 /l 1 ) and values of l 1 H corresponding to resonant wave amplification in the lower layer. However, while the propagating wave drag is larger in approximately hydrostatic conditions (i.e. large l 1 a), trapped lee wave drag is larger when l 1 a is of O(1).
This happens because, for the pressure perturbations associated with trapped lee waves to be optimally efficient, they must have a scale that matches the mountain width. This aspect deserves further study. • When the jump in the Scorer parameter is due to N only (i.e. U 1 /U 2 = 1), the maxima either of the propagating drag or of the trapped lee wave drag occur near l 1 H/π = 0.5 + n, where n is an integer. When the jump is only due to the discontinuity in U (i.e. U 1 /U 2 = l 2 /l 1 ), the drag maxima are located instead near l 1 H/π = n. In intermediate situations (U 1 /U 2 ≈ (l 2 /l 1 ) 1/2 ), the propagating wave drag is not substantially amplified, since the reflection coefficient at z = H takes values near zero. However, the modulation and possibility of amplification of trapped lee wave drag persists, since the corresponding waves always have a reflection coefficient of 1 at z = H.
The linear results calculated for the case with U 1 /U 2 = 1 were compared with similar results from numerical simulations, showing very good agreement.
The rich behaviour revealed by these calculations shows that there are many situations where the trapped lee wave drag may exceed the propagating wave drag, and the hydrostatic drag for a single-layer atmosphere. Since the latter is the reference value generally used in drag parametrizations, the present results emphasize the need to account for trapped lee wave drag in such parametrizations. As trapped lee waves are essentially non-hydrostatic, and therefore forced by relatively narrow mountains, their correct representation will remain a relevant issue even as the horizontal resolution of large-scale meteorological models increases. Obviously, the mean atmospheric profiles considered here may be viewed as no more than crude approximations to real profiles with sudden vertical variations. A more systematic study of trapped lee wave drag would require considering much more complex vertical variations of N and U, either continuous or comprising more layers. Depending on their local rates of variation, continuous atmospheric profiles should modify or suppress partial wave reflections, leading to changes in the propagating wave drag. Since the trapping height for trapped lee waves then necessarily depends on the wavenumber, the resonant trapped lee wave drag amplification process is likely to become weaker. However, its basic mechanism should remain essentially the same as considered here, relying on constructive interference of upward and downward propagating waves.
It would also be useful to take into account flow over threedimensional mountains and frictional effects, the latter of which are actually crucial for representing dissipative processes mediating the reaction force exerted by the mountains on the atmosphere. To address those situations using linear theory, more sophisticated mathematical techniques would be necessary.
As mountains with realistic heights typically generate nonlinear trapped lee waves, and the amplitude of these waves is known to be substantially larger than predicted by linear theory (Vosper, 2004), the concepts outlined here gain still higher practical relevance. 0 φ 2 φ 1 π/2 Figure B1. Schematic diagram of the integration path for the integral in (B2) along the real axis. The imaginary positive semi-plane is assumed to exist above this line. Two singularities, corresponding to two trapped lee wave modes, are denoted by φ 1 and φ 2 . Their order is consistent with Figures 2, 7 and 11, and (B1), i.e. wave modes of higher order correspond to a lower wavenumber, and thus to a lower value of φ.
From (B2), the condition for the location of singularities is just where j is any positive integer number. At the discrete points determined by (B4), the tangent function in (B2) becomes infinite. The roots of (B4), which will be called φ j (j being the corresponding integer index), must be found numerically. This method of finding the singularities has the advantage over that of Scorer (1949) that the roots need to be sought only within the fixed interval [0, π/2]. The next step is to evaluate the integral in (B2) at these singularities. It may shown that, with the addition of friction to this problem, the singularities move to the negative imaginary semi-plane. This means that the integration path in the inviscid problem under consideration must be indented above the singularities ( Figure B1). If the integration is carried out in this way, then D 2 = 4π 2 ρ 0 U 6 1 U 6 2 l 2 1 − l 2 2 3/2 × j h(φ j ) 2 sin φ j cos 2 φ j U 4 1 sin 2 φ j +U 4 2 cos 2 φ j 5/2 1− ∂m 1 ∂φ (φ j )H , where ∂m 1 ∂φ (φ) = − n 2 (φ)U 2 1 U 2 2 U 4 1 sin 2 φ + U 4 2 cos 2 φ .