# Mountain Waves Produced by a Stratified Shear Flow with a Boundary Layer. Part III: Trapped Lee Waves and Horizontal Momentum Transport

Clément Soufflet, François Lott, Bruno Deremble

► To cite this version:

Clément Soufflet, François Lott, Bruno Deremble. Mountain Waves Produced by a Stratified Shear Flow with a Boundary Layer. Part III: Trapped Lee Waves and Horizontal Momentum Transport. Journal of the Atmospheric Sciences, 2022, 79 (6), pp.1601-1614. 10.1175/JAS-D-21-0263.1. insu-03726891

HAL Id: insu-03726891

https://insu.hal.science/insu-03726891v1

Submitted on 21 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Mountain Waves Produced by a Stratified Shear Flow with a Boundary Layer. Part III: Trapped Lee Waves and Horizontal Momentum Transport

CLÉMENT SOUFFLET,ª FRANÇOIS LOTT,ª AND BRUNO DEREMBLEª

ª Laboratoire de Météorologie Dynamique, PSL Research Institute, Ecole Normale Supérieure, Paris, France

(Manuscript received 5 October 2021, in final form 11 February 2022)

**ABSTRACT:** The boundary layer theory for nonhydrostatic mountain waves presented in Part II is extended to include upward-propagating gravity waves and trapped lee waves. To do so, the background wind with constant shear used in Part II is smoothly curved and becomes constant above a “boundary layer” height *d*, which is much larger than the inner layer scale *δ*. As in Part II, the pressure drag stays well predicted by a gravity wave drag when the surface Richardson number *J* > 1 and by a form drag due to nonseparated sheltering when *J* < 1. As in Part II also, the sign of the Reynolds stress is predominantly positive in the near-neutral case (*J* < 1) and negative in the stable case (*J* > 1) but situations characterized by positive and negative Reynolds stress now combine when *J* ~ 1. In the latter case, and even when dissipation produces positive stress in the lower part of the inner layer, a property we associated with nonseparated sheltering in Part II, negative stresses are quite systematically found aloft. These negative stresses are due to upward-propagating waves and trapped lee waves, the first being associated with negative vertical flux of pseudomomentum aloft the inner layer, the second to negative horizontal flux of pseudomomentum downstream the obstacle. These results suggest that the significance of mountain waves for the large-scale flow is more substantial than expected and when compared to the form drag due to nonseparated sheltering.

**KEYWORDS:** Atmosphere; Downslope winds; Gravity waves; Kelvin-Helmholtz instabilities; Mountain waves; Orographic effects; Topographic effects; Boundary layer; Idealized models

## 1. Introduction

Low-level orographic drag which results from the interaction between mountain waves and the atmospheric boundary layer has a significant impact on the general circulation of the atmosphere (Pithan et al. 2016; Elvidge et al. 2019). However, this interaction is still not well understood and, in this regard, the example of climate models is instructive (see Lott et al. 2020, 2021, Parts I and II henceforth). In them the impact of mountains on (i) the boundary layer and (ii) the mountain gravity waves dynamics is actually handled by two distinct parameterizations: one for neutral flows (or small mountains), and one for stably stratified flow (or big mountains) (Beljaars et al. 2004; Lott and Miller 1997; Parts I and II). In this three-part study, we are trying to unify the theory of flow-topography interaction in the different regimes in a simple case where the slopes are small and eddy diffusivity represented by constant viscosity coefficients. With this framework, we can do a thorough investigation of the interaction and of the transition from neutral to stratified conditions and we view this approach as mandatory before considering larger slopes and more sophisticated turbulent closure.

In Parts I and II, we focused on the case where the background wind vanishes at the surface, and where the background wind shear $u_0$, and stratification $N^2$ are constant. In

this context, dissipation controls the dynamics over an inner layer which thickness is about 5 times the “inner” layer scale

$$ \delta = \left( \frac{\nu L}{u_0} \right)^{1/3}, \qquad (1) $$

with *L* the characteristic length of the obstacle, and *ν* the constant viscosity coefficient.

In Part I, we analyzed the wave-boundary layer interaction in the hydrostatic case and showed that for small mountains the wave stress is extracted from the inner layer instead of the ground surface as in the inviscid case: the large-scale flow is accelerated near the surface within the inner layer to balance the gravity wave drag. We also showed that the surface pressure drag and the Reynolds stress amplitude are well predicted using linear inviscid gravity wave theory as long as we take for the incident wind its value around the inner layer scale.

In Part II we examined the nonhydrostatic case, and more precisely, we studied the transition from stratified conditions to neutral conditions (small Richardson number). In the neutral case, we found that surface drag is well predicted by a form drag due to nonseparated sheltering. Henceforth, we refer to this situation as the “form drag regime” (see Parts I and II). In this case, the Reynolds stress profile is also maximum near the top of the inner layer indicating that the mean flow is decelerated in the lower part of the inner layer and accelerated in the upper part. For more stable flows (larger Richardson number), we recover the results from Part I for which internal waves control the dynamics: the surface pressure drag is well predicted by a wave drag, and the Reynolds stress accelerates the large-scale flow at the bottom of the

<footnote>ª Denotes content that is immediately available upon publication as open access.</footnote>

<footnote>Corresponding author: François Lott, flott@lmd.ens.fr</footnote>
inner layer. A major difference between the hydrostatic case and the nonhydrostatic case though is that all the upward gravity waves are reflected back toward the surface in the nonhydrostatic case whereas gravity waves travel only upward in the hydrostatic case. So in the hydrostatic case, gravity waves decelerate the far field flow (generally referred to as the gravity wave drag regime), whereas in the nonhydrostatic case, gravity waves decelerate the upper part of the inner layer.

We showed in Part II that the transition from form drag to wave drag regimes occurs for values of the Richardson number $J \sim 1$. Indeed, the turning height of the dominant wave-number $\bar{k} = 1/L$ is around $\sqrt{J}L$, such that wave dynamics can develop aloft and over the mountain when $J$ is large, whereas it is somehow inhibited when $J$ is small. Hence, we observed that when $J \approx 1$, the waves are reflected at altitudes about the length of the hill, and so, they are still close to the mountain when they return to the surface. It results destructive and constructive interactions between the wave induced pressure fields and the orography which produces low and high drag states, respectively.

However, a limitation of Parts I and II is that we excluded trapped lee waves from our analysis. Indeed, trapped lee waves cannot develop in constant shear flow, in part because pure trapped modes are related to neutral modes of Kelvin-Helmholtz (KH) instability (Lott 2016, hereafter L16; Soufflet et al. 2019), and so to emerge, such modes require that the Richardson number $J$ varies in the vertical according to the Miles-Howard theorem (Miles 1961; Howard 1961). Trapped lee waves are important because they can transport momentum in the horizontal direction only (Bretherton 1969), and this horizontal transport can be as significant as the one due to upward-propagating mountain waves (Teixeira et al. 2013). To reconcile such an horizontal transport of momentum and the noninteraction Eliassen-Palm theorem (Eliassen and Palm 1961), one simply has to translate this momentum transfer into pseudomomentum fluxes (Lott 1998; Georgelin and Lott 2001; see also Broad 2002; Héreil and Stein 1999).

Again, and now for trapped lee waves, the example of the coarse resolution models (e.g., resolutions ranging from around 50 to 200 km; Sandu et al. 2019) is instructive. Tsiringakis et al. (2017) showed that trapped lee waves can impact synoptic systems and near-surface meteorology as much as the blocked-flow drag (Lott and Miller 1997) or the turbulent orographic form drag (Beljaars et al. 2004). From observations, Steeneveld et al. (2009) also estimated that trapped waves induced substantial drag in the boundary layer.

The purpose of the present paper (Part III) is to study the impact of trapped lee waves when they coexist with upward-propagating waves. To permit these two type of waves we introduce a curvature in the background wind. Because boundary layer winds are generally small near the surface and present significant curvature near the top of the boundary layer, we will use this curvature to define a boundary layer height $d$ (which should not be confused with $\delta$, the inner layer scale over which waves are affected by dissipation). We will only consider boundary layers that are thicker than the inner layer ($d > \delta$). In this configuration, we will analyze how the

boundary layer depth influences the transition between the form drag regime and the wave drag regime, and also how the presence of the boundary layer impacts the Reynolds stress vertical profiles. We will also point out the role of trapped lee waves in this transition, and quantify their contribution to the wave drag.

The framework of this paper is close to the one used in L16. However, it is important to underline two major differences. First, in L16, the dynamics is inviscid and does not take into account the viscous dissipation in the boundary layer. Second, the influence of the boundary layer height ($d$) and the stability of the flow ($J$) will be here investigated independently which was not the case in L16 where the static stability was kept constant and the Richardson number was changed by varying the value of $d$.

The remainder of this paper is organized as follows. In section 1, we adapt the theoretical model from Part II to include an incident wind profile with a variable shear. In section 2, we study the impact of variable shear on the wave field and drags induced by the mountain. In section 3, we explain the onset of lee waves in the model. In section 4, we explain how lee waves contribute to the interaction between the mountain and the large-scale flow. Last, we analyze pseudomomentum budgets in section 5.

## 2. Theoretical framework

### a. Background profiles

The setup used here is close to the one used in Parts I and II, so we only recall here the salient features and emphasize the differences. For instance, the background wind and density profiles are now given by

$$u_0(z) = u_{0,z} \tanh(z/d), \quad \rho_0(z) = \rho_r + \rho_{0,z}, \qquad (2)$$

where the surface wind shear $u_{0,z}$ and stratification $\rho_{0,z}$ are both constant. $u_0(z)$ is displayed in Fig. 1a) where the values $d = 1$ km and $u_{0,z} = 10^{-3}$ s$^{-1}$ correspond to those chosen in the dimensional simulations (see section 1c). We choose this particular profile to represent the mean wind in the boundary layer because it is solution of the viscous equations near the surface but becomes constant above $d$, allowing a fraction of the mountain waves to propagate upward without being reflected. In the remainder of this analysis, we will refer to $d$ as the boundary layer height. As we shall see, this wind profile supports the existence of pure trapped lee waves, at least when $J < 0.25$ and in the inviscid limit, a dynamic that was completely absent in Parts I and II. Topography is still represented by a 2D Gaussian ridge of characteristic length $L$

$$h(x) = He^{-x^2/(2L^2)}, \qquad (3)$$

and is shown in Fig. 1b for the characteristic height $H = 150$ m and length $L = 1$ km used in the dimensional simulations. To interpret our results and clarify the differences in dynamics we also show in Fig. 1b a bulk representation of the turning layers calculated as the altitudes where
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic presenting (a) the flow profiles and (b) mountain parameters used in the dimensional simulations done with the MITgcm and when $d = 1$ km. (b),(c) Also shown is a bulk representation of the turning layers (TL), which we qualitatively defined as the layers in which the dominant harmonics forced by the mountain in Eq. (3) can potentially encounter a turning height and be reflected toward the ground [see Eq. (4)].</figcaption></figure>

$$ \frac{0.04}{L^2} < -g \underbrace{\frac{\rho_{0z}}{\rho_r u_0^2}}_{S_c(z)} - \underbrace{\frac{u_{0z}}{u_0}}_{S_c(z)} < \frac{5}{L^2}. \quad (4) $$

In (4), $S_c(z)$ is the Scorer (1949) parameter [see also Eq. (12)], and the bounds $0.04/L^2$ and $5/L^2$ have been chosen because in the spectral band $0.2/L < k < \sqrt{5}/L$ the mountain slope $[(kHL/\sqrt{2\pi})e^{-k^2L^2/2}]$ is above 20% of its maximum value at $k = 1/L$. In other words, we can expect the wave forcing to be substantial in this band. As we have seen in Part II, if the situation is near neutral, as is the case in Fig. 1b for $J = 0.1$, the reflection region is near the surface, and if the flow is stratified, as is the case when $J = 4$ in Fig. 1b, then the turning layer is well above the surface and the waves return to the surface far downstream. In the presence of nonconstant wind shear (Fig. 1c) we can use the same definition, except that at the altitude where the wind becomes almost constant, typically above $z = 2d$, very few additional turning levels occur above, so we can take $z = 2d$ as an upper bound of the turning layer. Hence, the situations for small $J$ (Fig. 1c) is comparable to the constant wind case, with near-neutral dynamics near the surface, the main difference being that the curvature of the wind also favor the emergence of trapped waves. For large $J$ the turning layer is narrow: the dynamics is stratified near the surface and many modes can propagate vertically in the far field without being reflected downward.

## b. Linear model

We now recall the main equations of the linear model. As in Part II, we use the nondimensional scalings

$$ \begin{aligned} (x, z) &= L(\bar{x}, \bar{z}) ; (u', w') = u_{0z}L(\bar{u}, \bar{w}) ; p' = \rho_r u_{0z}^2 L^2 \bar{p} ; \\ b' &= g \frac{\rho'}{\rho_r} = u_{0z}^2 L \bar{b}, \end{aligned} \quad (5) $$

where the prime denotes eddy flow with respect to the background profile, and the overbars are used for nondimensional variables, $x$ and $z$ are the horizontal and vertical dimensions and $u, w, p, \rho$, and $b$ are the horizontal and vertical velocities, the pressure, the density, and the buoyancy, respectively. With this scaling, the 2D Boussinesq linear equations, under the Prandtl approximation, are

$$ \bar{u}_0 \partial_{\bar{x}} \bar{u} + \bar{u}_{0z} \bar{w} = - \partial_{\bar{x}} \bar{p} + \bar{v} \partial_{\bar{z}} \bar{z}^2 \bar{u}, \quad (6a) $$

$$ \bar{u}_0 \partial_{\bar{x}} \bar{w} = - \partial_{\bar{z}} \bar{p} + \bar{b} + \bar{v} \partial_{\bar{z}} \bar{z}^2 \bar{w}, \quad (6b) $$

$$ \bar{u}_0 \partial_{\bar{x}} \bar{b} + J \bar{w} = P^{-1} \bar{v} \partial_{\bar{z}} \bar{z}^2 \bar{b}, \quad (6c) $$

$$ \partial_{\bar{x}} \bar{u} + \partial_{\bar{z}} \bar{w} = 0, \quad (6d) $$

in which

$$ \bar{u}_0(\bar{z}) = D \tan h(\bar{z}/D). \quad (7) $$

In that context no slip boundary conditions are

$$ \begin{aligned} \bar{h}(\bar{x}) + \bar{u}(\bar{x}, \bar{h}) &= 0, \quad \bar{w}(\bar{x}, \bar{h}) = 0, \quad \text{and} \quad J\bar{h}(\bar{x}) + \bar{b}(\bar{x}, \bar{h}) = 0 \quad \text{at} \\ \bar{h} &= Se^{-\bar{x}^2/2}. \end{aligned} \quad (8) $$

In Eqs. (6)–(8),

$$ J = -\frac{g\rho_{0z}}{\rho_r u_{0z}^2}, \quad P = \frac{\nu}{\kappa}, \quad S = \frac{H}{L}, \quad D = \frac{d}{L}, \quad \text{and} \quad \bar{v} = \frac{\nu}{u_{0z}L^2} \quad (9) $$

are a Richardson number, a Prandtl number, a slope parameter, a nondimensional boundary layer depth, and an inverse Reynolds number, respectively. With this new background flow profile the action budget is of form
$$ \frac{\partial}{\partial \bar{x}} \left[ \bar{u}_0 \underbrace{\left( \frac{\bar{\zeta} \bar{b}}{J} - \frac{\bar{u}_{0\bar{z}\bar{z}} \bar{b}^2}{2J^2} \right)}_{A} + \frac{\bar{b}^2}{2J} + \frac{\bar{u}^2 - \bar{w}^2}{2} \right] + \frac{\partial}{\partial \bar{z}} \left( \frac{\bar{u} \bar{w}}{F^{\bar{z}}} \right) = \frac{\bar{v}}{J} \underbrace{\left[ \bar{b} \partial_{\bar{z}}^2 \bar{\zeta} + P^{-1} \partial_{\bar{z}}^2 \bar{b} \left( \bar{\zeta} - \bar{b} \frac{\bar{u}_{0\bar{z}\bar{z}}}{J} \right) \right]}_{Q}, \quad (10) $$

with $\bar{\zeta} = \partial_{\bar{z}}\bar{u} - \partial_{\bar{x}}\bar{w}$ the vorticity, A the pseudomomentum, $F^x$ and $F^z$ the horizontal and vertical fluxes of pseudomomen- tum, and Q the pseudomomentum production/destruction by dissipative processes.

As in Part I [Eqs. (10) and (11)], we search inflow solutions that are linear, and express them in Fourier space in the hori- zontal direction. For instance, Eq. (6a) here transforms into

$$ \bar{i} \bar{k} \bar{u}_0 \bar{u} + \bar{u}_{0\bar{z}} \bar{w} = -i \bar{k} \bar{p} + \bar{v} \partial_{\bar{z}}^2 \bar{u}, \quad (11) $$

where the boldface notation is used for variables in the Fou- rier space.

For high Reynolds number $\bar{v} \ll 1$, the dynamics is inviscid at leading order, each harmonics $\bar{w}$ satisfies a Taylor-Gold- stein equation of the form,

$$ \bar{w}_{\bar{z}\bar{z}} + \left[ \bar{S}_c(\bar{z}) - \bar{k}^2 \right] \bar{w} = 0, \quad \text{where} \quad \bar{S}_c(\bar{z}) = \frac{J}{\bar{u}_0^2} + \frac{2}{D^2} \left( 1 - \frac{\bar{u}_0^2}{D^2} \right), \quad (12) $$

is the nondimensional expression of the Scorer parameter in (4): $(\bar{S}_c = L^2 S_c)$. We find the solution of Eq. (12) using appro- priate change of variables (see appendix and Lott et al. 1992) and we get

$$ \bar{w}_I = 2^{-m} r^{1/4+i\mu/2} (1-r)^{-m/2} W_{2(1)} \underset{\bar{z}\to\infty}{\approx} e^{-m\bar{z}/D}, \quad (13) $$

where $r = \tanh^2(\bar{z}/D)$. In (13),

$$ \mu = \sqrt{\left|J - \frac{1}{4}\right|} \quad \text{and} \quad m = \sqrt{J - D^2 k^2}, \quad (14) $$

where $m$ is the vertical wavenumber. Note that $\mu$ and $m$ are changed in $i\mu$ and/or $-im$, when $J < 1/4$ and/or $k^2 D^2 - J < 0$, respectively. Note also that the hydrostatic approximation is simply derived by omitting the horizontal wavenumber $\bar{k}$ in Eqs. (12) and (14).

In Eq. (13), $W_{2(1)}$ can be expressed in terms of hypergeo- metric functions, and the solutions for $k < 0$ are constructed by using the complex conjugate of the solutions obtained with $k > 0$. Near the surface the inviscid solution has an asymptotic behavior of the form,

$$ \bar{w}_I(\bar{k}, \bar{z}) \underset{\bar{z} \to 0}{\approx} \bar{w}_M(\bar{k}, \bar{z}) = \bar{a}_1(\bar{k}) \bar{z}^{1/(2-i\mu)} + \bar{a}_2(\bar{k}) \bar{z}^{1/(2+i\mu)}, \quad (15) $$

where $\bar{w}_M$ is a matching function and $\bar{a}_1(k)$ and $\bar{a}_2(k)$ are coef- ficients given in the appendix (they are independent of $\bar{k}$ in the hydrostatic approximation).

Because the background wind profile near the surface is close to the one used in Part II, we derive the viscous solution in the boundary layer in a similar way as in Part II: we define a nondimensional inner layer depth

$$ \bar{\delta} = \left( \frac{\bar{v}}{\bar{k}} \right)^{1/3}, \quad (16) $$

which represents the scale over which waves are affected by dissipation. In this region, a viscous solution $\bar{w}_V$ is derived numerically that satisfy the lower boundary condition Eq. (8) and that matches $\bar{w}_M$ when $\bar{z}/\bar{\delta} \to \infty$:

$$ \bar{w}_v [(\bar{k}, \bar{z}) / \bar{\delta}] \underset{\bar{z} / \bar{\delta} \to \infty}{\approx} f_{12}(\bar{k}) \bar{w}_M [(\bar{k}, \bar{z}) / \bar{\delta}]. \quad (17) $$

In (17) $f_{12}(\bar{k})$ are proportionality coefficients imposed by the lower boundary condition and that control the disturbance amplitude in the outer region. From these three solutions $(\bar{w}_I, \bar{w}_M, \bar{w}_V)$ we construct a uniform approximation for $\bar{w}$,

$$ \bar{w}(\bar{k}, \bar{z}) = f_{12}(\bar{k}) [\bar{w}_I(\bar{k}, \bar{z}) - \bar{w}_M(\bar{k}, \bar{z})(\bar{k}, \bar{z})] + \bar{w}_V(\bar{k}, \bar{z}), \quad (18) $$

with similar expressions for the horizontal wind, buoyancy, and pressure.

### c. Nonlinear model

As in Parts I and II, we will compare the theoretical model against nonlinear simulations using the MITgcm model (Mar- shall et al. 1997). The configuration of the model is almost the same as in Parts I and II, the values for the surface wind shear ($u_{0\bar{z}} = 1 \times 10^{-3}$ s$^{-1}$), boundary layer depth ($d = 1$ km), moun- tain length ($L = 1$ km), and height ($H = 150$ m) are those pre- sented in Fig. 1, which yield $S = 0.15$ (see Fig. 1). We impose a constant vertical temperature gradient, the temperature is then related to the density via a linear equation of state and we adjust the vertical stratification $N^2$ to match the nondimen- sional values of $J$: from $N^2 = 10^{-8}$ s$^{-2}$ ($J = 0.01$) to $N^2 = 1.6 \times 10^{-5}$ s$^{-2}$ ($J = 16$). The horizontal size of the domain is extended to 100 km and the total height of the domain is set to 50 km. This is a bigger domain than in Parts I and II to allow lee waves to propagate downstream and avoid numeri- cal instability. We initialize the model with the background flow and run it forward in time for 24 h (until we reach a steady state) with a time step of 0.2 s. We use a sponge layer active above 15 km and at the lateral boundaries to relax the dynamic variables to the prescribed upstream profiles [Eq. (2)]. We use a stretched grid to have maximum resolu- tion near the topography. The finest grid point has dimension
<figure><img src="image_3.png" /><figcaption>FIG. 2. Vertical velocity field for all simulations, $S = 0.15$. Each row corresponds to a value of $J$. (left three columns) Different values of $D$, and (right) the hydrostatic case with $D = 1$. In all panels, the contour interval CI = 0.004 and the negative values are dashed.</figcaption></figure>

of 205 m (horizontal) and 11 m (vertical) near the topography whereas the resolution coarsens to 715 m (horizontal) and 830 m (vertical) at the edge of the domain.

### 3. Upper-level and trapped waves

We plot in Fig. 2 the vertical velocity field for different values of the boundary layer depth $D$ and of the surface Richardson number $J$. In each simulation $S = 0.15$, $Pr = 2$, and the height of the inner layer for the dominant harmonic $\bar{k} = 1$ is $\bar{\delta} = 0.1$. We also plot the hydrostatic results for $D = 1$ to emphasize the significance of the reflected waves and of the trapped lee waves.

The top four panels in Fig. 2 show the vertical velocities in the stratified case ($J = 4$). We choose to present the $J = 4$ case first because it corresponds to the first figures shown in Parts I and II, e.g., the hydrostatic case with constant shear in Part I (Fig. 1) and the nonhydrostatic case with constant shear in Part II (Fig. 1). For the smallest value $D = 0.5$, Fig. 2a shows a train of upward-propagating waves with a small downstream signal at low level. At upper levels the wave field extends downstream in comparison with the hydrostatic case (Fig. 2d) indicating that nonhydrostatic effects make a difference at high altitudes. To understand why the low-level signal is small in this case, we recall that the square of the vertical wavenumber is given by [see Eq. (14)]

$$ \bar{m}^2 = J - D^2 \bar{k}^2. \qquad (19) $$

So the only modes that encounter a turning altitude are those for which $\bar{k} > \sqrt{J/D}$. In Fig. 2a as the Gaussian mountain forces harmonics with dominant wavenumber around $\bar{k} \approx 1$, one sees that the evanescent modes are for $\bar{k} > \sqrt{J/D} = 4$, most harmonics are free to propagate in the far fields. As $D$ increases nevertheless the number of reflected waves increases (Figs. 2b,c) and the wave signal near the surface becomes more pronounced downstream. An interesting aspect is that when these waves return to the surface on the lee side, their phase lines tilt significantly in the direction of the shear. This is consistent with the fact that for large $J$, the mountain waves are absorbed at the surface in the stable cases (Lott 2007): the signal is dominated by downward-propagating waves being absorbed.

The second and third rows in Fig. 2 are for the two values of the surface Richardson number that characterize the transition between the stratified and neutral case in Part II (i.e., when $D = \infty$). At $J = 1.7$ we found in Part II that there is a resonant interaction between reflected waves and the surface that yields a very strong wave signal aloft and immediately downstream, whereas at $J = 0.7$, we found that the interaction is destructive and the disturbance field is evanescent. The fact that some gravity waves can now propagate upward to $z = \infty$ when $D$ is finite profoundly changes the response. The cases with $J = 1.7$ (second row in Fig. 2) differ little from the cases with $J = 4$ (first row) except that the overall direction of
<figure><img src="image_4.png" /><figcaption>FIG. 3. Horizontal profiles of nondimensional vertical $\bar{w}$ velocity at $\bar{z} = D$ for $S = 0.15$. Each panel corresponds to a value of $D$ and each line to a value of $J$; the profiles from the MITgcm are represented in gray. The curves are shifted vertically for clarity, and the amplitude scale for all curves are given by the arrow.</figcaption></figure>

propagation is more horizontal, consistent with the fact that far aloft, more modes are impacted by nonhydrostatic effects. When $J = 0.7$ (third row in Fig. 2), we still visualize a system of gravity waves, which was not the case in Part II (Fig. 3c). Most gravity waves are propagating up when $D = 0.5$ (Fig. 2i) but there is now a system of downstream and horizontally propagating waves near the surface. For these waves, the phase line are more vertical, which indicate that fewer waves are absorbed at the surface compared to the previous cases. The signature is very much like that of a trapped lee wave. When $D$ increases (Figs. 2j,k) these near-surface waves become more and more prominent, which is again consistent with the fact that less modes can propagate far aloft according to (19). Interestingly, when $D$ increases, the horizontal wavelength near the surface increases as well. Finally, for $J = 0.1$ (Figs. 2m-p), there are few upward waves: the near-surface signal dominates but remains overall small.

## 4. Lee waves

As shown in Fig. 2, a significant difference between Parts I and II and this study is the presence of trapped lee waves for small values of the Richardson number $J$. In this section we analyze the impact of the boundary layer height $D$ and stability $J$ on the onset of these trapped lee waves and compare the results with the nonlinear model.

We plot in Fig. 3 the horizontal profiles of vertical velocity at $\bar{z} = D$ for $S = 0.15$ in the theory (black) and in the MITgcm (gray). Each panel corresponds to a different value of $D$, and $J$ is decreasing from top to bottom. In Fig. 3, we see that weakly stratified flows ($J < 1$) favor the onset of trapped lee waves regardless of the value of $D$, due to weaker wave absorption at the ground. This result with the linear model (also well corroborated by nonlinear simulation) extends the quasi-inviscid theoretical framework of L16, who showed that near-surface critical-level absorption is an active dissipation

process. Hence, in the present study, the same mechanism is still at play even when viscous dissipation acts in the inner layer near the ground. Note also that, due to dissipation, the downstream extent of lee waves is reduced compared to the quasi-inviscid results in L16 (even for $J < 0.25$).

We also observe that the trapped lee wave signal is small when $J$ is small (for instance, when $J = 0.01$ and $J = 0.1$). This is actually in agreement with L16 who showed that trapped lee waves are also near-neutral modes of KH instability. Hence for the vertical profile of horizontal wind given in (7), these modes satisfy the dispersion relation:

$$ \bar{k}_T^2 = \frac{1 - \sqrt{1 - 4J}}{2D^2} \quad \text{when } J < 1/4. \qquad (20) $$

It follows that for near-neutral flow ($J \ll 1$), and for $D \ge 1$, the trapped lee waves have predominant wavelength $\bar{k}_T \ll 1$. Such wavelength corresponds to quite long disturbances which are not efficiently excited in our case. To illustrate this point, we can consider that for mountain waves the vertical velocity forcing scales as $\bar{k}_T \bar{s} h (\bar{k}_T) e^{-\bar{k}_T / 2}$ in Fourier space (the Fourier transform of the slope). If we take the case $J = 0.2$ and $D = 1$, then the resonant mode prediction (20) gives $\bar{k}_T \approx 0.5$ and the vertical velocity forcing of this wavenumber scales as $\bar{k}_T S e^{-\bar{k}_T / 2} \approx 0.5S$, whereas for $J = 0.1$, $\bar{k}_T \approx 0.1$, the forcing $\bar{k}_T S e^{-\bar{k}_T / 2} \approx 0.1S$, it is about 5 times smaller then at $J = 0.2$.

To support this interpretation, we plot in Fig. 4 the lee waves wavelength for different value of $J$ and $D$ as calculated with the dispersion relation (20), the theoretical model and the MITgcm. One sees a good agreement between the different values (cf. each line style of the same color). We also see that the increase in boundary layer height systematically increases the lee waves wavelength whereas the increase in stability tends to reduce it, which is consistent with (20). The theoretical model (solid line) slightly overestimates the wavelength for small values of $D$.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Lee waves wavelength function of stability calculated from Eq. (20) (dashed), linear model (solid), and MITgcm (dash-dotted). Each grayscale line stands for a different value of D.</figcaption></figure>

This difference might be because as the value of D decreases (and approaches $5\delta$), the boundary and inner layers overlap, questioning the validity of our asymptotic matchings. This also explains why we limit our study to $D \ge 0.5$. The above results corroborate the observational study of Ralph et al. (1997), where the increase of boundary layer height during daytime induces an increase of lee waves wavelength.

If we now return to Fig. 3, another interesting point is that low-level oscillations can be found when D and J are large (see, for instance, Fig. 3c for $J = 6$ and 9). This is because when J and D are large some modes with $\bar{k} \approx 1$ can be reflected back to the surface [see Eq. (19)]. However, since this reflection occurs at high altitude (D is large), they return to the surface further downstream (we already noticed that in Figs. 2b,c). In this case the lee waves signal near the surface results from waves reflected downward in the lee side (referred as trapped waves or reflected waves in the remainder of this analysis) and do not correspond to trapped lee waves in the sense that they are not related to free modes of oscillation that exist in the inviscid case.

## 5. Pressure drag and Reynolds stress

To evaluate the effects of the wave field on the mean flow, we plot in Figs. 5a and 5b the surface pressure drag Dr along with the minimum and maximum of the mountain wave stress $\overline{F}^z$:

$$Dr = - \int_{-\infty}^{+\infty} \bar{p}(\bar{x}, \bar{h}) \frac{\partial \bar{h}}{\partial \bar{x}} d\bar{x}, \quad \bar{u} \bar{w} = \overline{F}^z = \int_{-\infty}^{+\infty} \bar{u} \bar{w} d\bar{x}. \quad (21)$$

These diagnostics are scaled using the drag predictor derived in Part II,

$$D_{rP} = \text{Max}(1, \sqrt{J})\delta(1)S^2/2. \quad (22)$$

We recall that the idea behind this formulation is to scale mountain drag as a form drag in weakly stratified cases ($J < 1$) where enhanced diffusion reduces pressure when the air parcels pass

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Surface pressure drag and (b) minimum and maximum of Reynolds stress, for different boundary layer depths $D = 0.5, 1, 2$, and 4, and for $S = 0.15$. All values are normalized by $D_{rP}$ [Eq. (22)].</figcaption></figure>

over the obstacle, and as a wave drag, when the flow is more stratified ($J > 1$), due to vertical propagation of gravity waves.

In Fig. 5a, we see that the drag predictor gives a rather good estimate of the surface pressure drag in a large range of flow stability J and boundary layer depth D. The best performance of the predictor is for $D = 1$ (black dashed line). For smaller value (for instance, $D = 0.5$) the form drag predictor overestimates the drag when $J < 0.1$. This is consistent with the fact that for small D, only few harmonics are confined near the surface. So for small D, these “long” harmonics contribute less to the near-surface dynamics responsible of the form drag than for larger D. When $D > 1$, we recover the behavior found in Part II where the transition zone around $J = 1$ presents strong variations in pressure drag. For instance, for $D = 4$ in Fig. 5a we recover the behavior found in Fig. 2 of Part II ($D = \infty$), with a pronounced low drag amplitude near $J = 0.7$ and a large drag amplitude near $J = 1.7$.

Interestingly, the transition from neutral to stratified flow when D is large occurs more smoothly when $D \approx 1$ (less amplitude between the lowest and highest value of the drag during the transition). To understand this behavior, we recall again that the variations in drag around $J = 1$ and in the constant shear case are caused by the fact that (i) all the reflected waves return to the surface, (ii) all harmonics encountering turning levels affect the surface pressure near the mountain downstream. For large J the waves are also reflected but the turning levels are sufficiently high that the wave “packet” return to the surface too far downstream to impact the drag (Fig. 1b). In the variable shear case, a good fraction of the harmonics excited by the mountain can propagate vertically
<figure><img src="image_7.png" /><figcaption>FIG. 6. (a),(d) Streamfunction, (b),(e) total wind vector and contours of $\bar{u}\bar{w}$ (negative values are dashed), and (c),(f) profiles of horizontally averaged Reynolds stress $\bar{\bar{u}}\bar{w}$, $S = 0.175$, $D = 4$; (a)–(c) $J = 9$, (d)–(f) $J = 0.1$.</figcaption></figure>

without being reflected. For instance, the dominant wavenumber $\bar{k} = 1$ only encounter turning levels when $J/D^2 < 1$. It means that when $D = 1$ the harmonics with $\bar{k} < 1$ no longer reaches turning level and that the turning levels of the harmonics with $\bar{k} > 1$ are located at higher altitude than in the constant shear case (see Fig. 1c). For these two reasons, the interaction between the reflected waves and the mountain are less significant compare to the case $D = \infty$ in Part II.

Figure 5b shows the minimum and maximum of horizontally averaged Reynolds stress $\bar{\bar{u}}\bar{w}$ (normalized by the predictor). These normalized extrema indicate how the wave field interacts with the mean flow. When the vertical profile of Reynolds stress presents a minimum at a given height, the mean flow is accelerated below that height, and decelerated above (the so-called gravity wave drag) and this situation corresponds to a wave drag regime. On the contrary, when the vertical profile of Reynolds stress is maximum at a given height, the mean flow is decelerated below that height, this situation corresponds to a form drag regime. Before discussing these regimes in detail, it is worth recalling that these changes in sign of the Reynolds stress have a profound dynamical origin. To illustrate it qualitatively, we show in Fig. 6 two cases with $D = 4$ and $S = 0.175$ (strong slope). In the first case, the flow is strongly stratified ($J = 9$) and is characterized by upstream blocking and downslope winds (Figs. 6a,b, respectively). In the downslope wind region where $\bar{w} < 0$ the disturbance in horizontal wind $\bar{u} > 0$, yielding the correlation $\bar{\bar{u}}\bar{w} < 0$ (see contours in Fig. 6b). When we average $\bar{\bar{u}}$ horizontally, we get a negative Reynolds stress (Fig. 6c). In the second case shown in Fig. 6, the flow is near neutral ($J = 0.1$) the dynamics is characterized by upslope winds upstream

and nonseparated sheltering downstream as illustrated by the streamfunction and the wind fields in Figs. 6d and 6e. In the sheltered zone the horizontal wind is weaker, so the disturbance wind is mostly $\bar{u} < 0$. Because this zone is located on the lee of the mountain, we also have a negative vertical velocity $\bar{w} < 0$, so the product $F^z = \bar{\bar{u}}\bar{w} > 0$ in a large sector behind the hill top, as shown in Fig. 6b. Averaged horizontally, this yields a positive Reynolds stress (see Fig. 6c).

If we now return to the extrema in $F^z$ in Fig. 5b, one sees that positive and negative extrema can occur simultaneously in the near-neutral cases ($J < 1$), at least when $D \le 1$. This strongly contrasts with what we found in Part II (or for $D = 4$ here) where form drag and wave drag do not occur simultaneously. For values of $D < 4$, one sees that form drag and wave drag are no longer exclusive of each other, clearly here the presence of trapped waves and the fact that more waves can propagate aloft when $D$ is small extent the domain over which the gravity waves dynamics contribute to the interaction between the orography and the large-scale flow. In Fig. 5b we also see that positive extrema only occur for $J < 0.4$ when $D = 0.5$, which means that in the presence of a thinner boundary layer the transition from neutral to stratified flow occurs for smaller values of the surface Richardson number $J$.

To assess the validity of this result, we now compare the linear model with the fully nonlinear model (MITgcm). For conciseness, we summarize this comparison in Fig. 7 using again the diagnostics of the extrema of the Reynolds stress (Fig. 7a), and also the index constructed in Part II [see Eq. (28) there]: we recall that this index discriminates between the regime of downslope sheltering versus the regime of upstream blocking. Since the results for $D = 4$
<figure><img src="image_8.png" /><figcaption>FIG. 7. Diagnostics from the MITgcm runs for S = 0.15 and for different values of the boundary layer depths D. In all panels, the corresponding results from the theory are shown in thin solid lines. (a) Normalized extrema in momentum flux. (b) Downslope sheltering vs upstream blocking index defined as the ratio between the max downslope wind amplitude and the max upslope wind amplitude [Eq. (28) from Part II]:</figcaption></figure>

$$ \underbrace{\text{Max}}_{\bar{z}<2\bar{h}/3, 0<\bar{x}<2} \sqrt{(\bar{z}+\bar{u})^2+\bar{w}^2} / \underbrace{\text{Max}}_{\bar{z}<2\bar{h}/3, -2<\bar{x}<0} \sqrt{(\bar{z}+\bar{u})^2+\bar{w}^2}. $$

correspond to D = ∞ in Part II, we only present the aforementioned diagnostics for D = 0.5, 1, and 2. For all these indicators, we see that the nonlinear model is in good agreement with the linear theory. Last we also observe that the sheltering versus blocking index does not seem to depend on the value of D.

## 6. Pseudomomentum budget

We have shown that in the presence of an inner layer and a boundary layer, form drag and wave drag coexist. This coexistence directly impacts the structure of the vertical profile of the Reynolds stress. We have also seen that the presence of a finite boundary layer depth enriches the inviscid dynamics, with trapped waves developing downstream the topography. We now provide more insight on the way these waves redistribute momentum not only in the vertical but also in the horizontal direction. To visualize this redistribution of momentum, we plot in Fig. 8 contours of the vertical action flux component [F^x defined in Eq. (10)] along with the total action flux vector for different value of the Richardson number J and boundary layer depth D.

In all panels in Fig. 8 one sees near the ground a downward flux on the upstream side of the ridge (dashed lines) and an upward flux on the downstream side (solid lines). This dipole

structure in the lower part of the inner layer is characteristic of the dynamics at work in our three-part paper and that we could refer to as linear dissipative, or weakly nonlinear dissipative. The key point is that when the mountain is in the inner layer, waves pseudomomentum is extracted from the inner layer rather from the surface as it occurs in the inviscid case. Near the top of the inner layer (i.e., around $\bar{z} = 5\delta$) and above, one sees in Fig. 8a that for small D and J = 1.7 the pseudomomentum flux vector points downward, as in the hydrostatic case in Part I, such that trapped waves (present, for instance, in Fig. 2e and to less extent in Fig. 3a) contribute little to the action flux. For larger D in Fig. 8b the reflected waves downstream produce an upward pseudomomentum flux, also slightly oriented upwind on the lee side, as if trapped waves where transferring momentum laterally rather than vertically. This larger contribution of trapped waves to the pseudomomentum budget is consistent with the fact that for J = 1.7 and D = 2 in Fig. 2g the low-level wave signal is quite substantial. For smaller J (Figs. 8c,d), trapped lee waves seem to contribute further in the far field, at least when J = 0.7, consistent with the fact that for small J, the mountain waves are less absorbed at the surface.

To provide a more quantitative estimate of the lateral fluxes due to the reflected and/or the trapped lee waves we next evaluate pseudomomentum fluxes through horizontal and vertical boundaries that encapsulate well the entire ridge. More specifically, we calculate the pseudomomentum fluxes outgoing from the top hat defined by the three segments,

$$ (-\bar{X}, 0) \times (-\bar{X}, \bar{Z}), (-\bar{X}, \bar{Z}) \times (+\bar{X}, \bar{Z}), (+\bar{X}, \bar{Z}) \times (+\bar{X}, 0), \quad (23) $$

where the capitals letter are used to distinguish the bound of the integration domains from the coordinates (see an illustration of this box in Fig. 8a for $\bar{X} = 3$, and $\bar{Z} = 0.3$). We always take $\bar{Z} > \bar{h}$ and $\bar{X} > 3$, the latter condition guarantees that $\bar{h}(\pm\bar{X}) \approx 0$. The integral of the pseudomomentum fluxes across the boundaries writes

$$ P^x(\bar{X}, \bar{Z}) = \int_0^{\bar{Z}} F^x(\bar{X}, \bar{z}) d\bar{z} \text{ and } P^z(\bar{X}, \bar{Z}) = \int_{-\bar{X}}^{\bar{X}} F^z(\bar{x}, \bar{Z}) d\bar{x}, \quad (24) $$

$$ P^{\text{out}}(\bar{X}, \bar{Z}) = P^x(\bar{X}, \bar{Z}) - P^x(-\bar{X}, \bar{Z}) + P^z(\bar{X}, \bar{Z}). \quad (25) $$

As the in-going flux is always small, we will only discuss the fluxes along the upper and downstream sides of the box. The solid lines in Fig. 9 are the vertical profiles of the total outgoing momentum fluxes, $P^x(\bar{X}, \bar{Z}) + P^z(\bar{X}, \bar{Z})$ for three different downstream locations: one near the mountain $\bar{X} = 3$, one further downstream $\bar{X} = 5$, and one very far downstream $\bar{X} = 25$. We selected the first two positions to illustrate the large erosion of the emitted pseudomomentum fluxes ($P^{\text{out}}$) that occur just downstream the hill (i.e., between $\bar{X} = 3$ and $\bar{X} = 5$). And we selected $\bar{X} = 25$ to measure the total erosion occurring in the boundary layer (for $\bar{X} = 25$ we found that the lateral
<figure><img src="image_9.png" /><figcaption>FIG. 8. Contours of vertical action component ($F^z$; negative values are dashed), along with total action vector for $S = 0.15$. For illustration in (a) are the limits of a characteristic box used to calculate the emitted PM fluxes.</figcaption></figure>

pseudomomentum flux is almost always null, see the thick gray dotted lines in Fig. 9). The first thing to notice is that for all values of $D$ and $J$, the total flux of pseudomomentum $P^{\text{out}}(\bar{X}, \bar{Z}) = \text{const.}$ on the vertical when $\bar{Z} > 5\bar{\delta}$, i.e., when the upper bound $\bar{Z}$ is in the inviscid region. This can be viewed as a generalization of the Eliassen-Palm theorem in the presence of trapped lee waves. Also, and when $D$ is small ($D = 0.5$ in Figs. 9a,c), almost all the pseudomomentum flux is transmitted vertically through the boundary layer: the lateral fluxes of pseudomomentum are always small, at least below $\bar{Z} = D$ (see the dotted lines). Above the boundary layer, the lateral propagation of the gravity waves in the inviscid region produces substantial horizontal fluxes when the downstream distance is not too large ($\bar{X} = 3$ and $\bar{X} = 5$, thick and thin dotted black lines, respectively; see also Fig. 2). When $D = 2$ in Fig. 9b one sees that the total pseudomomentum fluxes diminishes in amplitude when $\bar{X}$ increases and in the inviscid zone $\bar{Z} > 5\bar{\delta} = 0.5$. This diminution is due to the fact that for large values of $D$, there are more reflected waves returning into the inner layer than when $D$ is smaller. Moreover, these reflected waves are associated with positive vertical fluxes of pseudomomentum $F^z$. Therefore when the horizontal extension of the upper bound of the box increases, these reflected waves cancel the negative contribution of the upward waves in the integral flux $P^z$. This mechanism combines with a substantial contribution of the trapped lee waves propagating horizontally and at lower level when $J = 0.7$ in Fig. 9d. In this case one sees that the amplitude of the vertical flux first increase between $5\bar{\delta} < \bar{Z} < D$ (above the inner layer but inside the boundary layer) when $\bar{X}$ increases, consistent with the fact that the reflected waves are less absorbed when $J$ decreases.

As seen in Parts I and II, and repeated here, it is quite difficult to pin the location of extraction of pseudomomentum. It is not entirely extracted from the surface as in the inviscid case (Durran 1995; Lott 1998), and it is not extracted from the viscous fluid in the inner layer as in the case with $S \ll \bar{\delta}$ (Part I). Because of these difficulties, we instead propose to diagnose the largest amount of pseudomomentum that is produced by the interaction between the mountain and the inner layer. We call it the emitted pseudomomentum, and evaluate it as the total pseudomomentum flux going out of the top hat defined by (23) with lateral boundary near the downstream foot of the hill $\bar{X} = 3$ and upper boundary at the altitude $Z_T$ that minimizes the outgoing flux:

$$P^{\text{Emit}} = P^x(3, \bar{Z}_T) + P^z(3, \bar{Z}_T) = \min_{S<\bar{Z}<\infty} [P^x(3, \bar{Z}) + P^z(3, \bar{Z})]. \qquad (26)$$

A typical box through which we measured the emitted PM flux is shown in Fig. 8a for illustration. To measure how much of this emitted flux goes in gravity wave drag far aloft and to estimate how much stays at low level, we compare $P^{\text{Emit}}$ to the far field Reynolds stress, and to the minimum in Reynolds stress,

$$\overline{F^z}(\bar{Z} = \infty) \text{ and } \min_{S<\bar{Z}<\infty} \overline{F^z}, \qquad (27)$$

respectively. As we defined the Reynolds stress using the vertical component of action $F^z$ [see Eq. (10)], and because its minimum is always found well inside the inner layer, we can conservatively consider the difference between the emitted
<figure><img src="image_10.png" /><figcaption>FIG. 9. Vertical profiles of $P^x$ (dashed), $P^z$ (dotted), and their sum $P^{\text{out}} = P^z + P^x$ (solid). Shades indicate different horizontal box sizes: $\bar{X} = 3$ (thick black), $\bar{X} = 5$ (thin black), and $\bar{X} = 25$ (thick gray).</figcaption></figure>

pseudomomentum flux and the minimum in Reynolds stress as the part due to horizontal propagation yielding to wave large-scale flow interactions occurring inside the inner layer and downstream exclusively. Figure 10 shows that when $D$ increases the emitted pseudomomentum fluxes and the minimum in Reynolds stress are quite different. More specifically, we find that for narrow boundary layers ($D = 0.5$ in Fig. 10a) lateral fluxes are small. We note also that the minimum in Reynolds stress is larger than the far field Reynolds stress when $J > 1$, the upward-propagating waves are dissipated in the upper part of the inner layer (see also Part I). This effect occurs for all boundary layer depth $D$ but is rather weak for small $J$ (the far field and minimum Reynolds stress almost coincide for $J < 1$). As $D$ increases up to $D = 1$ and for moderate stability $0.1 < J < 1$, the presence of lee waves induces a lateral flux: the emitted flux $P^{\text{Emit}}$ substantially exceeding the Reynolds stresses in magnitude (see Fig. 10b). This lateral contribution decreases as $J$ increases, the waves being absorbed at the surface. For larger $D$ in Figs. 10c and 10d, the contributions of the lateral fluxes become more pronounced, and for large $J$ a good part of the lateral fluxes are due to the fraction of the emitted waves, which are reflected downward (the surface absorption being large).

## 7. Conclusions

In this paper we have analyzed how a background wind curvature, which mimics a boundary layer of depth $D$, modulates

<figure><img src="image_11.png" /><figcaption>FIG. 10. Emitted pseudomomentum, minimum, and emitted value of the Reynolds stress. In all simulations $S = 0.15$.</figcaption></figure>

the impact of small-scale mountains on the large-scale flow while staying in the weakly nonlinear and dissipative regime used in Parts I and II. We found that trapped lee waves develop much more than in the constant shear case, they resemble to Kelvin–Helmholtz modes of instability, at least when the surface Richardson number $J < 0.25$. This corroborates the results in L16 and Soufflet et al. (2019) but using another boundary layer parameterization and another fully nonlinear model [the MITgcm here versus WRF in Soufflet et al. (2019)]. We also found that for large $J$ and $D$, low-level waves are related to modes that have been reflected at turning levels and that return to the surface downstream where they are absorbed. For small $J$, the trapped lee waves may not be efficiently excited, simply because the corresponding modes of KH instability have small horizontal wavenumber compared to the characteristic scale of the mountain (a condition that writes $\bar{k} \ll 1$ in dimensionless form).

As in the constant wind shear case we recover the transition from the form drag regime to the wave drag regime when the flow stability near the surface increases. The wave drag regime is associated with downslope winds and upstream blocking and is characterized by a negative Reynolds stress which mostly radiates in the far field (see Figs. 6a–c). The form drag regime is associated with upslope winds and downstream sheltering and is always associated with positive
Reynolds stress, confined to the inner layer (see Figs. 6d–f). One key result of this part is that there exists a transition zone for which these two regimes coexist. For this intermediate situations, the Reynolds stress is positive in the lower part of the inner layer and negative in the upper part and aloft: the interaction between the boundary layer flow and the mountain produces deceleration near the surface, acceleration in the middle of the inner layer, and deceleration (gravity wave drag) near the top of the inner layer and above. As a direct consequence, we can measure the relative importance of the form drag regime and of the wave drag regime by comparing the minima and maxima of the Reynolds stress.

The nature of this transition is controlled by the number of reflected waves that return to the surface and by the absorption properties of the surface. In this paper, we controlled the reflected waves with the nondimensional boundary layer depth D, while we controlled the absorption with the surface Richardson number J. When D is small, most harmonics are free to propagate in the far field, and upward-propagating gravity waves control the dynamics for values of J > 0.1. When D increases the background wind curvature starts supporting horizontally propagating trapped lee waves when J ≈ 1. For larger values of J, these trapped lee waves do not develop well (the ground absorption is too large), but there can be vertically propagating waves returning from the far field to the surface where they are absorbed. We showed that, when they exist, trapped lee waves and reflected waves can produce significant lateral fluxes of momentum downstream the mountain. Pseudomomentum budget near the topography indicates that lateral and vertical momentum flux are on the same order of magnitude for intermediate values of J. These downstream fluxes remain substantial up to 5 times the mountain width, the associated lee wave drag being applied in the inner layer.

As said in the introduction, our formalism is still too simple for direct use to change subgrid-scale orography parameterizations in large-scale models. Nevertheless, some directions and extensions to the 3D case are detailed in the conclusions section in Part II. For instance, the fact that we have “form” and wave drag predictors [Eq. (22)] and some indications of where the drags should be deposited could be helpful for low hills. Here, we also have learned that with trapped waves the drag predictors are still accurate, and that the trapped waves interact with the large-scale flow in the inner layer and not below the turning layer. In the future, we also plan to combine in GCMs the so-called turbulent orographic form drag scheme (Beljaars et al. 2004) and the subgrid-scale orography schemes (Lott and Miller 1997) by calculating a separation scale between them. For this, and based on our results, we could measure explicitly the scale L at which the dynamics at the top of the inner layer passes from neutral to stratified, e.g., when

$$ \frac{u_0[\delta(L)]}{N[\delta(L)]L} = 1. \quad (28) $$

This equation is an extrapolation of our results since it translates into $\bar{\delta} = J$, the value $J = 1$ being a good measure

of the neutral to stratified transition in our case. In such an equation, we should also adapt the evaluation of the depth $\delta$. For instance, we could choose for it the altitude at which disturbance advection equilibrates disturbance dissipation (which is the definition of the inner layer depth in the viscous case). As an example and for completeness, we repeat Eq. (30) in conclusions of Part II, where we considered a turbulent closure based on mixing length theory. If $\Lambda(z)$ is the mixing length profiles of the undisturbed flow, $\delta$ could be defined as

$$ \frac{u_0(\delta)}{L} = 2 \frac{\Lambda(\delta)^2}{\delta^2} \left\| \frac{du_0}{dz}(\delta) \right\|. \quad (29) $$

*Acknowledgments.* This work was supported by the VESRI Schmidt Future project “DataWave,” and by the Laboratoire de Recherche Conventionné “Yves Rocard,” a collaborative unit between CEA and Ecole Normale Supérieure. The MITgcm used in the present study can be downloaded at mitgcm.org. The theoretical model is written in Fortran, and is available on F. Lott web page: www.lmd.jussieu.fr/flott/homepage.html.

# APPENDIX

## Calculation of the Outer Solution

If we change variables and take $r = \tanh^2(\bar{z}/D)$, Eq. (12) transforms into

$$ \frac{d^2 \bar{w}}{dr^2} + \left( \frac{1}{2r} - \frac{1}{1-r} \right) \frac{d \bar{w}}{dr} + \left[ \frac{J}{4r^2(1-r)^2} + \frac{1}{2r(1-r)} - \frac{D^2 \bar{k}^2}{4r(1-r)^2} \right] \bar{w} = 0. \quad (A1) $$

This equation has three regular singular points in $r = 0, 1, \infty$, when $k^2 D^2 - J > 0$ and $J > 1/4$ their exponent pairs are

$$ \begin{aligned} r = 0: \quad & \alpha_1 = \frac{1}{4} + i \frac{\mu}{2}, & \alpha_2 = \frac{1}{4} - i \frac{\mu}{2}; \\ r = 1: \quad & \gamma_1 = -\frac{m}{2}, & \gamma_2 = +\frac{m}{2}; \\ r = \infty: \quad & \beta_1 = 1, & \beta_2 = -\frac{1}{2}. \end{aligned} \quad (A2) $$

In (A2),

$$ \mu = \sqrt{|J - \frac{1}{4}|}, \quad \text{and} \quad m = \sqrt{|J - D^2 k^2|}, \quad (A3) $$

they are changed in $i\mu$ and/or $-im$, when $J < 1/4$ and/or $k^2 D^2 - J < 0$, respectively. Introducing the change of variable,

$$ \bar{w} = r^{\alpha_1} (1 - r)^{\gamma_1} W. \quad (A4) $$

Equation (A1) transforms into the hypergeometric equation, and the inviscid solution
$$ \bar{w}_I = 2^{-m} r^{1/4+i\mu/2} (1-r)^{-m/2} W_{2(1)} \underset{\bar{z} \to +\infty}{\approx} e^{-m\bar{z}/D}, \quad (A5) $$

behaves like a pure exponential function in the far field, for instance, like a unit amplitude upward-propagating gravity wave when $D^2\bar{k}^2 < J$. In (A5), the solution

$$ W_{2(1)} = (1-r)^m F(c-b, c-a; c-a-b+1; 1-r) \quad (A6) $$

is expressed using the hypergeometric function $F$, and the coefficients

$$ \begin{aligned} a &= \alpha_1 + \beta_1 + \gamma_1 = \frac{5}{4} + i\frac{\mu}{2} - \frac{m}{2}, \\ b &= \alpha_1 + \beta_2 + \gamma_1 = -\frac{1}{4} + i\frac{\mu}{2} - \frac{m}{2}, \\ c &= 1 + \alpha_1 - \alpha_2 = 1 + i\mu. \end{aligned} \quad (A7) $$

To evaluate $\hat{w}_c$ near the surface, the transformation (15.3.6) in Abramowitz and Stegun (1964) is used to express (A5) in terms of the solutions (15.5.3) and (15.5.4) in AS, e.g.,

$$ W_{1(0)} = F(a, b; c; r) \text{ and } W_{2(0)} = r^{1-c} F(a-c+1, b-c+1; 2-c; c; r): \quad (A8) $$

$$ W_{1(0)} = A_1 W_{1(1)} + A_3 W_{2(1)}, \quad W_{2(0)} = A_2 W_{1(1)} + A_4 W_{2(1)}, \quad (A9) $$

where

$$ \begin{aligned} A_1 &= \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}, & A_3 &= \frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}, \\ A_2 &= \frac{\Gamma(2-c)\Gamma(c-a-b)}{\Gamma(1-a)\Gamma(1-b)}, & A_4 &= \frac{\Gamma(2-c)\Gamma(a+b-c)}{\Gamma(a-c+1)\Gamma(b-c+1)}. \end{aligned} $$

This yields

$$ \begin{aligned} \bar{w}_I &= r^{\alpha_1} (1-r)^{\gamma_1} (b_1 W_{2(0)} + b_2 W_{1(0)}), \quad \text{where} \\ b_j &= (-1)^{j-1} \frac{2^{-m} A_j}{A_1 A_4 - A_2 A_3} \quad \text{for } j=1,2. \end{aligned} \quad (A10) $$

When approaching the surface, this inviscid solutions behaves as the matching function,

$$ \bar{w}_I(\bar{k}, \bar{z}) \underset{\bar{z} \to 0}{\approx} \bar{w}_M(\bar{k}, \bar{z}) = \bar{a}_1(\bar{k})\bar{z}^{1/2-i\mu} + \bar{a}_2(\bar{k})\bar{z}^{1/2+i\mu}, \quad (A11) $$

providing that

$$ \bar{a}_1 = D^{-1/2+i\mu} b_1 \quad \text{and} \quad \bar{a}_2 = D^{-1/2-i\mu} b_2. \quad (A12) $$

REFERENCES

Abramowitz, M. J., and I. A. Stegun, 1964: Handbook of Mathematical Functions, Dover, 1045 pp.

Beljaars, A. C. M., A. R. Brown, and N. Wood, 2004: A new parametrization of turbulent orographic form drag. Quart. J. Roy. Meteor. Soc., 130, 1327–1347, https://doi.org/10.1256/qj.03.73.

Breherton, F. P., 1969: Momentum transport by gravity waves. Quart. J. Roy. Meteor. Soc., 95, 213–243, https://doi.org/10.1002/qj.49709540402.

Broad, A., 2002: Momentum flux due to trapped lee waves forced by mountains. Quart. J. Roy. Meteor. Soc., 128, 2167–2173, https://doi.org/10.1256/003590002320603593.

Durran, D. R., 1995: Pseudomomentum diagnostics for two-dimensional stratified compressible flow. J. Atmos. Sci., 52, 3997–4009, https://doi.org/10.1175/1520-0469(1995)052<3997:PDFTDS>2.0.CO;2.

Eliassen, A., and E. Palm, 1961: On the transfer of energy in stationary mountain waves. Geofys. Publ., 22, 1–23.

Elvidge, A. D., and Coauthors, 2019: Uncertainty in the representation of orography in weather and climate models and implications for parameterized drag. J. Adv. Model. Earth Syst., 11, 2567–2585, https://doi.org/10.1029/2019MS001661.

Georgelin, M., and F. Lott, 2001: On the transfer of momentum by trapped lee waves: Case of the IOP3 of PYREX. J. Atmos. Sci., 58, 3563–3580, https://doi.org/10.1175/1520-0469(2001)058<3563:OTTOMB>2.0.CO;2.

Hérel, P., and J. Stein, 1999: Momentum budgets over idealized orography with a non-hydrostatic anelastic model. I: Two-dimensional flows. Quart. J. Roy. Meteor. Soc., 125, 2019–2051, https://doi.org/10.1002/qj.49712555806.

Howard, L. N., 1961: Note on a paper of John W. Miles. J. Fluid Mech., 10, 509–512, https://doi.org/10.1017/S0022112061000317.

Lott, F., 1998: Linear mountain drag and averaged pseudo-momentum flux profiles in the presence of trapped lee waves. Tellus, 50A, 12–25, https://doi.org/10.3402/tellusa.v50i1.14509.

—, 2007: The reflection of a stationary gravity wave by a viscous boundary layer. J. Atmos. Sci., 64, 3363–3371, https://doi.org/10.1175/JAS4020.1.

—, 2016: A new theory for downslope windstorms and trapped lee waves. J. Atmos. Sci., 73, 3585–3597, https://doi.org/10.1175/JAS-D-15-0342.1.

—, and M. Miller, 1997: A new subgrid scale orographic drag parameterization: Its formulation and testing. Quart. J. Roy. Meteor. Soc., 123, 101–127, https://doi.org/10.1002/qj.49712353704.

—, H. Kelder, and H. Teitelbaum, 1992: A transition from Kelvin-Helmholtz instabilities to propagating wave instabilities. Phys. Fluids, 4A, 1990–1997, https://doi.org/10.1063/1.858368.

—, B. Deremble, and C. Soufflet, 2020: Mountain waves produced by a stratified boundary layer flow. Part I: Hydrostatic case. J. Atmos. Sci., 77, 1683–1697, https://doi.org/10.1175/JAS-D-19-0257.1.

—, —, and —, 2021: Mountain waves produced by a stratified shear flow with a boundary layer. Part II: Form drag, wave drag, and transition from downstream sheltering to upstream blocking. J. Atmos. Sci., 78, 1–35, https://doi.org/10.1175/JAS-D-20-0144.1.

Marshall, J., A. Adcroft, C. Hill, L. Perelman, and C. Heisey, 1997: A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers. J. Geophys. Res., 102, 5753–5766, https://doi.org/10.1029/96JC02775.

Miles, J. W., 1961: On the stability of heterogeneous shear flow. J. Fluid Mech., 10, 496–508, https://doi.org/10.1017/S0022112061000305.

Pithan, F., T. G. Shepherd, G. Zappa, and I. Sandu, 2016: Climate model biases in jet streams, blocking and storm tracks
resulting from missing orographic drag. *Geophys. Res. Lett.* **43**, 7231–7240, https://doi.org/10.1002/2016GL069551.

Ralph, F. M., P. J. Neiman, T. L. Keller, D. Levinson, and L. Fedor, 1997: Observations, simulations, and analysis of nonstationary trapped lee waves. *J. Atmos. Sci.*, **54**, 1308–1333, https://doi.org/10.1175/1520-0469(1997)054<1308:OSAAON>2.0:CO2.

Sandu, I., and Coauthors, 2019: Impacts of orography on large-scale atmospheric circulation. *npj Climate Atmos. Sci.*, **2**, 10, https://doi.org/10.1038/s41612-019-0065-9.

Scorer, R. S., 1949: Theory of waves in the lee of mountains. *Quart. J. Roy. Meteor. Soc.*, **75**, 41–56, https://doi.org/10.1002/qj.49707532308.

Soufflet, C., F. Lott, and F. Damiens, 2019: Trapped mountain waves with a critical level just below the surface. *Quart. J. Roy. Meteor. Soc.*, **145**, 1503–1514, https://doi.org/10.1002/qj.3507.

Steeneveld, G.-J., C. J. Nappo, and A. A. Holtslag, 2009: Estimation of orographically induced wave drag in the stable boundary layer during the CASES-99 experimental campaign. *Acta Geophys.*, **57**, 857–881, https://doi.org/10.2478/s11600-009-0028-3.

Teixeira, M. A. C., J. L. Argain, and P. M. A. Miranda, 2013: Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere. *Quart. J. Roy. Meteor. Soc.*, **139**, 964–981, https://doi.org/10.1002/qj.2008.

Tsiringakis, A., G.-J. Steeneveld, and A. Holtslag, 2017: Small-scale orographic gravity wave drag in stable boundary layers and its impact on synoptic systems and near-surface meteorology. *Quart. J. Roy. Meteor. Soc.*, **143**, 1504–1516, https://doi.org/10.1002/qj.3021.