Orographic gravity‐wave drag over multiple bell‐shaped mountains

Mountains exert a drag on the atmospheric flow through breaking gravity waves, known as orographic gravity‐wave drag (OGWD). This article presents a theoretical analysis of the impacts and limitations of missing processes in current parameterization schemes for OGWD and examines the significance of various approximations. We derive OGWD analytically for multiple three‐dimensional bell‐shaped mountains, while previous studies have derived OGWD for a single three‐dimensional elliptic mountain. We take into account the geometric overlap between the mountains, vertical shear of the horizontal wind, nonhydrostatic atmosphere, and Coriolis forces. Our findings show that OGWD depends on the direction of the incoming wind and the degree of overlap between two mountains. OGWD increases with the Richardson and Rossby numbers; however, it decreases with the Froude number. Even when the mountain is isotropic, an associated transverse force is generated if the sheared flow passes over the mountain on the f$$ f $$ ‐plane. By calculating OGWD analytically from more complex terrain and airflow, this study can contribute to the development of an advanced OGWD parameterization in the numerical models used for climate and weather prediction.


INTRODUCTION
Mountains exert a force on the atmosphere, which leads to the generation of vertically propagating gravity waves.This force redistributes horizontal momentum vertically, transports or dissipates energy, and affects the angular momentum budget of the atmosphere, along with decelerating atmospheric circulation (Sandu et al., 2019).Collectively, the drag from mountain waves influences weather and climate (Miller et al., 1989;Palmer et al., 1986) by changing temperature, chemical species, and also clouds in the polar stratosphere (Okamoto et al., 2011;Solomon, 1999) through complex feedback processes.However, these intricate interactions and effects are challenging to resolve directly in many atmospheric models, due to the models' coarse horizontal resolution.This limitation necessitates the use of certain parameterizations to represent the impact of mountain drag effectively on a larger scale.An accurate parameterization of mountain drag is important to improve the performance of weather and climate models because mountain drag exists in a near-persistent way.
Mountain drag is represented in numerical models through several distinct processes: low-level flow blocking (Kim & Doyle, 2005;Lott & Miller, 1997), orographic gravity-wave drag (OGWD; Phillips, 1984;Kim & Arakawa, 1995;Choi & Hong, 2015), and turbulent orographic form drag (TOFD; Beljaars et al., 2004).Low-level flow blocking is the mountain drag generated when the wind blowing into the mountains splits and flows around the terrain in a vertically stratified atmosphere, whereas OGWD and TOFD are associated with the wind blowing over the terrain.OGWD and TOFD are generated via gravity waves and turbulence, respectively.These representations are highly idealized and act only as an approximation to the multiple interacting dynamic processes around orography.Notably, TOFD encompasses turbulence-induced drag and drag from small-scale orography (horizontal scales less than ∼5 km) that cannot generate propagating gravity waves.Consequently, existing models, assuming a sharp cut-off scale for nonpropagating waves, overlook features under 5 km.This simplification fails to consider the more complex interplay of nonhydrostatic scales, influenced by the large-scale flow.Our study can rectify these oversights, enhancing mountain drag representation in atmospheric models.
Many theoretical studies for OGWD have used a set of linearized governing equations for the steady, inviscid, and adiabatic atmosphere.The OGWD generated by a single three-dimensional mountain in the hydrostatic nonrotating atmosphere is proportional to the product of wind speed, horizontal length scale, and the square of the mountain's height (e.g., Phillips, 1984;Smith, 1980).In a constant shear flow over a single three-dimensional mountain, Teixeira et al. (2004) and Teixeira and Miranda (2006) demonstrated that OGWD decreases as the wind shear increases.Several studies have derived asymptotic expressions for OGWD, including the nonhydrostatic, rotating effect over two-dimensional bell-shaped mountains (Grisogono et al., 1993;Smith, 1979;Teixeira et al., 2008), as well as over a single three-dimensional bell-shaped mountain (Miranda & James, 1992;Xu et al., 2021).For a flow with constant wind and static stability, nonhydrostatic effects tend to reduce OGWD (Gill, 1982).Hydrostatic numerical simulations demonstrate that OGWD decreases with increasing rotating effects, for low mountains (Ólafsson & Bougeault, 1997).Additionally, several review articles have documented the physical mechanisms (e.g., Teixeira, 2014) and parameterizations (e.g., Alexander et al., 2010;Kim et al., 2003) of OGWD, as well as the associated mountain meteorology (e.g., Smith, 2019).
Nevertheless, previous studies have primarily been restricted to two dimensions or assumed a single three-dimensional orography.In contrast, using multiple bell-shaped orographies to represent the given topography is more realistic than employing a single elliptic orography.To our knowledge, no previous study has derived analytical solutions for the OGWD over multiple three-dimensional mountains, even in the simplest hydrostatic nonrotating conditions.Perhaps this is the reason why most state-of-the-art parameterizations for OGWD either assume a single elliptic orography derived from subgrid heights of the mountains (Lott & Miller, 1997;Phillips, 1984;Scinocca & McFarlane, 2000) or estimate the anisotropic parameters in a bulk form (Kim & Arakawa, 1995;Kim & Doyle, 2005;Xie et al., 2020).Moreover, the Coriolis effect on OGWD in the nonhydrostatic atmosphere in three dimensions or in an atmosphere with vertical wind shear has not been explored.
In our study, we derived analytical solutions for OGWD driven by randomly distributed multiple three-dimensional bell-shaped mountains using a set of linearized governing equations and explored the impacts of various effects that have not been examined previously by studies.In Section 2, starting from the equations of motion for the inviscid atmosphere, we formulate the wave drag over any three-dimensional orography.Then, we discuss how different effects, such as the overlapping effect, wind-shear effect, Coriolis effect, and nonhydrostatic effect, influence the drag force in Section 3. Finally, Section 4 contains the main conclusions of this study.

FORMULATION
We extend the work of Teixeira et al. (2004) to derive the mountain drag in an inviscid atmosphere, where only gravity waves contribute to the drag and no turbulence exists.This section presents the general formulation, and then the following section examines how the wind-shear effect, nonhydrostatic effect, and Coriolis effect influence the wave drag over multiple bell-shaped mountains.
The linearized equations of motion with the inviscid, adiabatic, steady condition and Boussinesq approximation take the form Here, the superscript hats denote the Fourier transform in the horizontal wavenumber space, (k, l) =∶ k is the horizontal wavenumber vector, the primes denote differentiation with respect to height, z, and i is the imaginary unit.(U, V) =∶ U is the mean incoming wind velocity (which is assumed to depend only on height z), (u, v, w) is the velocity perturbation associated with the waves, p is the pressure perturbation,  0 is a reference constant density, b = g∕ 0 is the buoyancy (where g is the acceleration of gravity,  is a potential temperature perturbation, and  0 is a reference constant potential temperature), and N 2 = (g∕ 0 )dΘ∕dz is the static stability (where Θ is the mean potential temperature, assumed to depend only on z).In Equations ( 1a), (1b), the Coriolis parameter is expressed as f , which is considered to be a constant.These five equations can be combined in order to yield something akin to the Taylor-Goldstein equation for the Fourier transform of the vertical velocity perturbation, ŵ: where The derivation of Equations ( 2)-( 3) is given in Appendix B. Note that U ′′ = V ′′ = 0 is required when the rotation effect of Earth (herein rotating effect) is considered, in order not to contradict the steady-state flow condition (Inverarity & Shutts, 2000;Teixeira, 2014).Further, the wavelike solution satisfying the radiation boundary condition, that is, ŵ → 0 as z → +∞, is where m(k, l, z) is the vertical wavenumber.Also, there is a surface boundary condition where the subscript 0 attached to the velocity (or its derivatives) denotes values of the velocity (or its derivatives) taken at z = 0. Therefore, Equations ( 2)-( 5) form a complete boundary-value problem of a second-order ordinary differential equation.If both N and (U, V) = (U 0 , V 0 ) are constant, m does not vary with height.However, if a vertical shear of the horizontal wind (herein wind shear) exists, the Wentzel-Kramers-Brillouin (WKB) approximation is required to solve the equation (see Bender et al., 1999;Teixeira et al., 2004).The WKB approximation is valid when the solutions vary slowly in the vertical.Thus, the vertical wavenumber m(z) is assumed to slowly vary in the vertical, for example,  ∶= d(1∕m)∕dz ≪ 1.If taking the zeroth-order of the vertical wavenumber, m 0 , we get where F = (U 0 k + V 0 l)∕N and R = (U 0 k + V 0 l)∕f are the local Froude number and the local Rossby number, respectively.Therefore, the magnitude of Equation ( 6) should be small for the validity of the WKB approximation.Subsequently, the admissible domain of the wavenumber vector should be mean wind speed and a is a horizontal length scale.Also, the model is observed to be formally valid for high Richardson number, low Froude number, and high Rossby number.
In general, the linearized surface pressure drag force exerted on the mountain by the atmosphere is given by where ∇ H is the horizontal gradient and the asterisk denotes a complex conjugate.Evidently, the key is to determine the pressure perturbation from the equations of motion when the model orography is given and evaluate the double integrals.The associated pressure perturbation p at the surface can be calculated by combining the equations of motion (Equations 1a,1b): If the solution of Equations ( 2)-( 5) is inserted into Equation (8), we get the pressure perturbation for each order in a small parameter , where the expressions are given in Appendix C. Subsequently, by substituting this pressure with Equation (7), we get F I G U R E 1 Schematic diagrams for three bell-shaped mountains (Equation 10). (a) The topography of three-dimensional bell-shaped mountains.The mountains are labeled 1, 2, and 3. (b) Density plot of the given topography.The black circles are the contour, and a n , where n = 1, 2, 3, is the radius of the contour with height 2 −3∕2 times the height of each individual peak of the mountain.(c) The cross-section between the centers of mountains 1 and 3. d 1,3 is the distance between the peaks of the mountains divided by a 1 + a 3 (Equation 12).
where each term on the right-hand side corresponds to the zeroth, first, and second-order term, respectively.The drag force is a function only of the mean horizontal wind (and its vertical gradient) at the surface and the geometry of the given topography.If wind shear does not exist, the zeroth-order term alone is sufficient to describe the total drag.The first-order term is nonzero if and only if the flow is rotating and a wind shear exists.When f ≠ 0, the imaginary part of the first-order term of the Fourier transform of the pressure perturbation at the surface does not vanish.This can generate drag, which is represented as the first-order term of the drag in Equation (9).To the authors' knowledge, the first-order term has never been expressed before.It plays an important role in the presence of the wind shear on the f -plane.The second-order term is simplified for the hydrostatic and nonrotating case.If either a nonhydrostatic or rotating effect is included, the expression is too complicated to integrate.We calculated the drag up to first order when a nonhydrostatic or rotating effect is considered.Otherwise, we use a second-order WKB approximation.This is the reason we place the round bracket in Equation ( 9) and integrate over R 2 .
In this work, multiple bell-shaped mountains are considered: for example, where # is the number of bell-shaped mountains.Here a n is the width from the center of each mountain, (x n , y n ); therefore, the height is 2 −3∕2 times the height of the top of the mountain, h n .Schematic diagrams of the mountains for # = 3 are shown in Figure 1.The result should have terms from each peak and from the interaction of every pair of bell-shaped mountains, because the drag force is quadratic with ĥ, that is, i≠ Re{ ĥi ĥ *  }, where ĥn denotes each bell-shaped mountain.The following section presents the drag formulae for just two bell-shaped mountains.This discussion can be extrapolated to any number of bell-shaped mountains.The formulae for an arbitrary number of bell-shaped mountains are detailed in Appendix A.

Hydrostatic, nonrotating flow
Generally, the pressure drag increases as the correlation between p and ∇ H h increases (Equation 7). Figure 2 illustrates the surface pressure fields caused by one hill and two hills, when hydrostatic, nonrotating flow is assumed.It shows horizontal cross-sections of the normalized pressure perturbation near the mountains when U(z) = (U 0 + z, 0).The value of  is irrelevant, because the pressure perturbation is normalized, but here the sign is  > 0. If a meridional wind blows, the pressure perturbation rotates (Teixeira et al., 2004).The contribution of the pressure field from each order is shown separately to highlight its different roles.This is done by taking the inverse Fourier transform of the present model (Equation 8) analytically, where the expressions for the surface pressure perturbation are given in Appendix C.
Here, we summarize each order and the total pressure perturbation for a single identical bell-shaped mountain; see Teixeira et al. (2004) for more details.The F I G U R E 2 Normalized pressure perturbation at the surface for a hydrostatic nonrotating flow as given by U = U 0 + z, V ≡ 0 (thick solid circle).Terrain elevation equal to h max ∕2.(a,b) Zeroth-order pressure; (c,d) first-order pressure; (e,f) second-order pressure; (g,h) total pressure when  > 0 and Ri = 1.The pressure perturbations are normalized by  0 NU 0 h, , and  0 NU 0 h for zeroth-, first-, and second orders and total pressure, respectively.asymmetry of the pressure perturbation field induces or reduces drag, because the mountain is isotropic.Figure 2a shows the zeroth-order pressure, which is antisymmetric.The upwind maximum and downwind minimum characteristics cause a positive drag force on the mountain (Smith, 1980).Only zeroth-order pressure appears if wind shear does not exist ( = 0); however, higher-order pressure terms appear and complicate the total drag if wind shear is present ( ≠ 0).The first-order pressure perturbation (Figure 2c) is symmetric for Boussinesq flow (Tang et al., 2007); hence, producing the drag force is irrelevant.The second-order pressure perturbation (Figure 2e) is antisymmetric, similar to the zeroth-order; however, it has the opposite sign, which contributes to decreasing drag.Figure 2g shows the sum of these three pressure contributions where  > 0 (positive shear) and Ri = 1.Teixeira et al. (2004) verified that the linear approach with the WKB approximation gives good agreement for Ri = (1) compared with the nonlinear numerical model simulation (Grubiši ć & Smolarkiewicz, 1997).
If more hills are added, the pressure field becomes more complicated.Since the pressure perturbation is linear to the orography, it can be presumed that the surface pressure field is caused by one hill sitting on top of the other hill.The four right panels in Figure 2 display each order and the total pressure perturbation for double identical bell-shaped mountains.As the pressure perturbations change intricately and the orography is no longer isotropic, the drag cannot be reduced to merely twice that of a single mountain; instead, we anticipate it will exhibit distinct characteristics.To substantiate these expectations, we will proceed by deriving the solution for OGWD directly.Notably, the traditional OGWD scheme, as described by Lott and Miller (1997), has assumed that the drag from several ellipses in a line is additive.From our work, however, one can see that this assumption is not correct.

Flow without wind shear
The closed-form formula of the OGWD can be obtained by integrating Equation ( 9) with the model orography from Equation (10).For hydrostatic, nonrotating flow without wind shear over two bell-shaped mountains, for example, F = 0, R −1 = 0, f = 0, and only a zeroth-order term survives in Equation ( 9) and # = 2 in Equation ( 10); the drag is then given by where On the right-hand side, the drag caused by each peak and the drag caused by the overlap (or interaction) of two bell-shaped mountains are shown.The first and second term, D peak , has been demonstrated in many studies over the past 50 years.The third term, D ovlp , is one of the main results in this work.D ovlp is a function of d, , and , as well as the geometry of orography and incoming mean wind speed.The existence of D ovlp indicates that, when multiple mountains are superposed, the total force is not simply the sum of the forces exerted by each peak.Additionally, D ovlp plays a significant role in generating a transverse force.By applying the change of coordinate in Equation ( 11), we can directly obtain a closed-form expression for the drag D and transverse force T on two hills: ) .
(13) Here, D > 0 indicates that the air experiences a force in the direction opposite to the incoming uniform flow.T > 0 represents a force acting in the right perpendicular direction to the incoming flow.As shown in Equation ( 13), a single-axis symmetrical mountain does not produce a transverse force.A transverse force only exists if d ≠ 0 and  −  ≠ n∕2, where n is an integer.The term "drag" has a dual meaning.First, it refers to the drag vector, denoted by D. This vector is represented by (D x , D y ) in the x-y coordinate system, or by (D, T) in the transformed coordinate system.Second, it specifically denotes a force that acts only in the direction opposite to the wind, represented by D in the transformed coordinate system.
Figure 3 shows the normalized drag and transverse force of OGWD for hydrostatic, nonrotating flow over two identical bell-shaped mountains.The normalization is performed using the magnitude of D peak .Consequently, the factor of 2 in Figure 3a represents the sum of the drag caused by each peak, and the deviation from 2 is due to the effect of D ovlp .Figure 3b   accurate.The numerical simulation from Vosper (1996) shows that the OGWD exerted on two three-dimensional mountains by a unidirectional wind is qualitatively similar to the two-dimensional result by Grisogono et al. (1993).Our results when  −  = 0 agree with the results of Vosper (1996), and we investigate further the effect of wind blowing from any direction on drag.
To summarize Figure 3, both forces are periodic and the corresponding period is both .Moreover, drag and transverse force have a reflection symmetry for every  −  = n∕2 and  −  = ∕4 + n∕2, respectively.Drag force is an even function, while transverse force is an odd function.D and T have extreme values at  −  = n∕2 and  −  = ∕4 + n∕2, respectively.For instance, when two mountains are aligned parallel with the wind ( −  = 0), the total drag is less than two distant mountains, for example, 71.This is because the air passes over the mountain and through the next mountain before returning to the original altitude, implying a decrease in the height gradient between the mountains.However, when two mountains are aligned perpendicular with the wind ( −  = ∕2), the total drag is greater than two distant mountains, for example, D > 2||D peak ||.This is due to the extra cross-sectional area created by the overlap of the two mountains, which blocks the air.Hence, it gives rise to an increase in the pressure gradient force.Note that D is always positive, which means that mountains always slow the speed of the wind and never accelerate it.In contrast, T can be any sign, depending on the direction in which the mountains are aligned.The wind can be bent to any side with respect to the incoming direction.
In summary, mountain anisotropy with various wind directions is important to improve the accuracy of OGWD, as many works agree with this rationale (e.g., Phillips, 1984;Smith & Kruse, 2018;Xie et al., 2020).Equation ( 13) and Figure 3 have the same physical features as eq.6.4 and fig.8 of Phillips (1984), which considered the flow over a mountain with an elliptical horizontal cross-section.This fact implies that one can create anisotropy by overlapping several isotropic mountains.When complex terrain is constructed by overlapping those basis mountains, the drag due to the superposed mountains can be obtained by adding all D peak and D ovlp .

Flow with wind shear
When wind shear is introduced into the incoming mean flow profile, the magnitude of the pressure perturbation around the mountain is reduced compared with the case with no wind shear (Figure 2e).Additionally, the effects of wind shear lead to a more horizontal dispersion of the waves.Therefore, the drag is expected to decrease slightly, and this effect must be included in the formula.The formula can be obtained by calculating Equation ( 9) up to the second order and the analytical result for the general flow over randomly arranged multiple mountains is given in Equation (A1).The wind curvature can lead to either an increase or a decrease in drag.As suggested by Equation (A1), the drag increases when the sign of U 0 U ′′ 0 , U 0 V ′′ 0 , V 0 U ′′ 0 , or V 0 V ′′ 0 becomes negative, due to an increase in the magnitude of the pressure perturbation around the mountain (Equation C8).In this section, we discuss the case of double bell-shaped mountains along the x-axis and let the flow have a linear variation, for example, U ′′ 0 = V ′′ 0 = 0. Subsequently, the zonal drag by peak D peak and interaction D ovlp reduces from Equation (A1) to respectively.The Richardson number for each x and y direction is defined as Ri and the inverse of the conventional Richardson number Ri is the sum of the inverse of Ri x and Ri y .The functions g 1 and g 3 are defined in Equation ( A2f)-(A2h).Equation ( 14) was already derived by Teixeira et al. (2004), and the wind-profile shear and curvature were found to cause notable corrections (Xu et al., 2020).The newly introduced Equation ( 15) refers to the additional drag by a linearly sheared flow when the mountains are superposed.Note that the sign of the wind shear does not matter for drag, because the inverse Richardson number is proportional to the square of the wind shear.
Figure 4 shows how D peak and D ovlp change with the inverse Richardson number for the x and y directions.D peak is normalized to the zonal drag generated by a constant flow over a single peak; therefore, the value at , where a moderate overlap is chosen (d = 1); hence, the value at 18.Note that the normalizing factor is not the interaction term induced by a constant flow shown in Equation ( 13).The intention is to emphasize that the magnitude of D peak is larger than that of D ovlp when d = 1, and we also wanted to specify that D ovlp is less than zero when d = 1.Evidently, the magnitude of the normalized D peak and D ovlp gets smaller as the inverse Richardson numbers increase.This implies that the normalized total drag decreases as wind shear increases, since the magnitude of D peak is larger than that of D ovlp , and D ovlp becomes ineffective if the shear is strong.Additionally, the zonal drag D peak is more sensitive to Ri −1 x than Ri −1 y , and D ovlp hardly changes as Ri −1 y varies.The change of D peak with wind shear is highly dependent on the incoming mean wind direction (V 0 ∕U 0 = tan ), while the variation of D ovlp with V 0 ∕U 0 is less significant.This is attributed to the fact that the magnitude of D ovlp is one order smaller than that of D peak .Drag changes depending on the wind direction in the presence of shear, because the pressure perturbation shown in Figure 2 rotates, and the degree of rotation varies across different orders.

Nonhydrostatic, nonrotating flow
For our model orography (Equation 10), we define the horizontal Froude number, Fr n =  0 ∕[Na n ] and Fr i =  0 ∕[N(a i + a  )] for each mountain and their interactions, respectively, where 1 ≤ n, i,  ≤ #.They are used to estimate nonhydrostatic effects.For example, if F ≠ 0, the calculation of Equation ( 9) requires much effort and creates additional functions of Fr n and Fr i that are multiplied by the hydrostatic, nonrotating drag.Note that, while the horizontal Froude number relates to whether a gravity wave is generated, the vertical Froude number,  0 ∕[Nh n ], relates to whether the air stream splits or goes over the mountain, which is not addressed herein.Additionally, the linear theory of mountain drag is known to be invalid for low vertical Froude number (Teixeira, 2014); thus, nonlinear drag becomes important.Some characteristics of nonlinear flow include wave breaking and flow blocking (Garner, 1995).
When the flow over a single bell-shaped mountain is nonhydrostatic and nonrotating (R −1 = 0), the asymptotic drag D peak can be obtained by approximating 9).Then, the drag is given by where P is a bunch of special functions represented in Equation (A7).The function P(Fr) is smaller than unity and takes the value 1 in the hydrostatic limit (Fr → 0).P(Fr) can be seen as the nonrotating drag over a single bell-shaped mountain normalized by its hydrostatic and nonrotating value, D 0 ∶=  4  0 N 0 ah 2 , that is, P(Fr) = ||D peak ||∕D 0 .Note that F 2 < 1 does not necessarily imply Fr < 1, allowing for the possibility of Fr > 1.This is because Figure 5a plots P(Fr) (pink solid curve) along with the exact drag (black solid curve) obtained from numerical integration (Equation 9) and the corresponding relative error (blue solid curve).The relative error converges to zero as Fr n goes to zero, increases as Fr n goes up, and has the value 0.15 at Fr n = 5.Our result is compared with another asymptotic expression (pink dashed curve) from the recent previous result, eq.18 of Xu et al. (2021).They presented the OGWD excited in a nonhydrostatic flow approximated by elementary functions, which is computationally efficient.The relative error (blue dashed curve) from their result also converges to zero at a low Fr n limit; and has the larger error when Fr n ≳ 1.Moreover, as seen in Figure 5a, the normalized drag decreases drastically when Fr n ≳ 1, which implies that the drag in a nonhydrostatic atmosphere should be much smaller than hydrostatic drag.This is reasonable, because a high Froude number implies that the time an air parcel takes to pass over the terrain (= a∕ ) is less than it takes for a vertical oscillation due to the buoyancy force (= 1∕N) (Lin, 2007).Hence, the air parcel is less affected by orography, or mountain waves are barely generated.Another reason that nonhydrostatic drag is weaker than hydrostatic drag is because waves disperse more horizontally in nonhydrostatic flow, causing the energy to spread more horizontally rather than vertically.When two bell-shaped mountains overlap, the additional term D ovlp appears as where Kc (x) and Ks (x) are given in Equation (A5).Equation ( 17) is given by the semi-analytical expression, and the coordinate system with the direction of the inflow and its perpendicular direction is chosen.Figure 5b,c displays how the normalized drag and transverse force according to the D ovlp change in Fr i and  − , respectively.A moderate overlap is chosen (d = 1), and they are divided by 8 0 N 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 to check only the properties of the result of the integral in Equation ( 17), not the wind speed or the geometry of the mountain.The pattern with the exact drag (Equation 9) is identical, and the corresponding absolute error is one or two orders less.D ovlp and T ovlp are periodic with period  as in the hydrostatic, nonrotating case (Figure 3).However, the reflection symmetry for T ovlp is broken if Fr i > 0. The  −  at which T ovlp has extreme values recedes from ∕4 + n∕2 as Fr i increases.If Fr i ≳ 1, D ovlp becomes less dependent on incoming wind direction, and the magnitude of normalized D ovlp becomes small.This is because mountain waves are hardly generated when the horizontal length scale is small enough.Considering normalized D peak is also small when Fr n ≫ 1 (Figure 5a), we calculated the ratio of the magnitude of D ovlp and T ovlp to the magnitude of D peak , and found that the ratio approaches 2 in drag and 0 in transverse at high Fr n limit.Therefore, a highly nonhydrostatic atmosphere rarely bends the mean wind despite the mountain anisotropy.

Hydrostatic, rotating flow
We define the Rossby number, Ro n =  0 ∕[fa n ] and Ro i =  0 ∕[f (a i + a  )], for each mountain and their interactions, respectively.They are used to estimate rotating effects.For example, if R −1 ≠ 0, the integration in Equation ( 9) creates an additional function of Ro that is multiplied by the hydrostatic, nonrotating drag.When the flow over a single mountain is hydrostatic (F = 0) and rotating, the drag is given by where 9); hence, Equation ( 18) is an exact drag representation of hydrostatic drag as a function of Rossby number.Also, R −2 < 1 does not necessarily imply Ro −1 < 1, allowing for the possibility of Ro −1 > 1.The first term (zeroth-order term) of Equation ( 18) is equivalent to Equation (A15) of Miranda and James (1992).However, if wind shear is present, the second term (first-order term) appears to be given by a semi-analytical expression.The terms multiplied by U 0 and U ⊥ 0 refer to drag and transverse force, respectively.The first-order term shows a production of transverse force if (U 0 , V 0 ) and (U ′ 0 , V ′ 0 ) are not perpendicular.This is remarkable because the transverse force can be created without the mountain anisotropy when the Coriolis force and the wind shear are present.
Figure 6a shows the normalized drag D peak for hydrostatic flow over a single mountain as a function of Ro −1 n , where the normalizing factor is D 0 =  4  0 N 0 ah 2 .The drag decreases drastically when Ro −1 n ≳ 1, which is equivalent to horizontal scales of the mountains a n exceeding the Rossby radius of deformation  ∕f .This implies that the drag in a rotating atmosphere should be much smaller than the nonrotating drag.Hence, the Coriolis force suppresses mountain wave drag, or the inertia-gravity wave barely produces a drag force.Note that this result is specific to conditions with low vertical Froude numbers.In contrast, at high vertical Froude numbers, the drag is observed to increase with rotation, attributable to nonlinear regimes associated with flow blocking (Ólafsson & Bougeault, 1997).
When perpendicular wind shear is observed, the normalized drag moves away from the zeroth-order drag (Figure 6a).The normalized first-order term is proportional to Ri −1∕2 ; therefore, the stronger the shear, the more effective first-order drag is.The changes in drag depend on  −  ′ .The mountain produces more drag if  <  −  ′ < 2 and less drag if 0 <  −  ′ < .All the results mentioned have suggested that an isotropic mountain does not create a transverse force when only one of wind shear, nonhydrostatic, or rotating effects exists.However, in the presence of combined Coriolis and wind-shear effects, the transverse force T peak is shown to act on the incoming mean flow (Figure 6b).The order of T peak ∕D 0 is (10 −2 ).The transverse force is positive if ∕2 <  −  ′ < 3∕2, and negative if −∕2 <  −  ′ < ∕2.
By examining the surface pressure field around a mountain in the presence of the Coriolis force, it becomes evident why zeroth-order drag is smaller than in a nonrotating situation and why first-order drag should occur in the presence of wind shear.Figure 7 shows the normalized pressure perturbation by the Coriolis force and the vertical wind shear of meridional wind, where the normalizing factor is the same as in Figure 2. The zeroth-order pressure perturbation is not affected by the wind shear; however, it changes on the f -plane (Figure 7a,b).Compared with Figure 2a,b, the pressure perturbation is generally reduced when the Coriolis effect is present.The upwind maximum and downwind minimum near the mountain are also noticed, like that shown in Figure 2a,b.However, the sign is changed far from the mountain, unlike Figure 2a,b.Consequently, the horizontal distribution of the pressure perturbation makes the drag smaller.The first-order pressure does not produce drag force without a Coriolis effect (Figure 2c,d).Only when the Coriolis effect acts on the perpendicularly sheared flow is the first-order pressure perturbation antisymmetric (Figure 7c,d); therefore, first-order drag appears.At large a, the additional normalized first-order pressure perturbation (Figure 7c,d) is much smaller than the reference normalized first-order pressure perturbation (Figure 2c,d).
When two bell-shaped mountains overlap, the additional term D ovlp appears as where the first term corresponds to the zeroth-order term and the second term corresponds to the first-order term.Each order term is given by the semi-analytical expression, and the coordinate system with the direction of the inflow and its perpendicular direction is chosen.As can be seen, the first order depends on  −  ′ as well as  − .5b,c, a moderate overlap is chosen (d = 1), and they are divided by 8 0 N 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 .D ovlp and T ovlp are periodic with period  as in the hydrostatic, nonrotating case (Figure 3).However, the reflection symmetry of T ovlp is broken if Ro −1 i > 0. The  −  at which T ovlp has extreme values recedes from ∕4 + n∕2 as Ro −1 i increases.The sign still changes frequently depending on  −  when the rotating effect is strong, which is in contrast to the sign changing less depending on  −  when the nonhydrostatic effect is strong (Figure 5b,c).The magnitude of normalized D ovlp becomes small if Ro −1 i ≳ 1.This is because the Coriolis force suppresses the wave drag when the horizontal length scale is large enough.Considering that normalized D peak is also small when Ro −1 n ≫ 1 (Figure 6a), we calculated the ratio of the magnitude of D ovlp and T ovlp to the magnitude of the zeroth order of D peak , and found that D peak converges much faster than D ovlp or T ovlp in the high Ro −1 n limit.Figure 8b,d displays the normalized drag and transverse force up to the first order; hence, it depicts the wind-shear effects.The degree of overlap and the normalizing factor are the same as above, and  −  ′ = ∕4 is set.Evidently, the vertical wind shear changes the value of drag and transverse force slightly.The reflection symmetry of D ovlp is broken if Ro −1 i > 0 and  ′ ≠ ∕2 + n, and rotational symmetry of T ovlp is broken if Ro −1 i > 0 and  ′ ≠ n.

Nonhydrostatic, rotating flow
In this subsection, the drag is calculated when both nonhydrostatic and rotating effects are considered.The drag for the general condition is obtained by approximating 9).The asymptotic expression of nonhydrostatic, rotating drag up to zeroth order is given by On the right-hand side, the first and second terms correspond to the drag by each peak and the interaction of two bell-shaped mountains, respectively, and are given by analytical and semi-analytical expressions, respectively.
Here, the bunch of special functions is named Kw, as defined in Equation (A4).The function Kw may be viewed as an inclusive expression of linearized drag force by all mesoscale waves (gravity waves and inertia-gravity waves) generated by a single bell-shaped mountain.In order to understand our results and estimate errors for peak terms, we have reproduced figs 1 and 2 of Teixeira et al. (2008) using our results, as below; this differs only in the dimension of the orography and thus the expression.Figure 9 plots the drag normalized by its hydrostatic and nonrotating value D 0 as a function of hill width, a, and shows our results for the geophysical range of variation in the values of N, f , and  .The nonrotating drag (pink curves) corresponds to Equation ( 16).The hydrostatic drag (green curves) corresponds to the first term of Equation ( 18).To apply the hydrostatic, nonrotating approximation to OGWD, it should be that the horizontal scale of the terrain a is not too small (i.e., nonhydrostatic effect), while the horizontal scale of the terrain a is not too large (i.e., rotating effect).The nonhydrostatic, rotating drag is shown by blue dotted curves (the exact drag; Equation ( 9)) and black curves (the approximate drag; Equation ( 20)).The blue and black curves pass both the green and pink lines; hence, the nonhydrostatic, F I G U R E 9 Drag normalized by its hydrostatic and nonrotating value as a function of the hill width, a. Upper curves: Solid line: asymptotic drag for nonhydrostatic and rotating flow (Equation 20).Blue: exact drag (Equation 9); green: hydrostatic drag (Equation 18); pink: nonrotating drag (Equation 16).The upper solid curve is cut off due to a machine underflow.
rotating drag can be seen as a smooth transition from nonhydrostatic, nonrotating drag to hydrostatic, rotating drag.The result emphasizes the need to consider both nonhydrostatic and rotation effects.The range of a that can be applied for D 0 is just about 2-30 km (upper curve), and even completely vanishes (lower curve) depending on  , N, and f .This fact explains the rationale for why general circulation models (GCMs) treat OGWD on the horizontal length scales of at least a few km.For example, the ECMWF model does not treat horizontal scales less than 5 km in the model because these are approximated to be nonhydrostatic (ECMWF, 2021).However, our work suggests that the ECMWF's OGWD scheme, by not accounting for the rotating effect, inevitably overestimates the drag force on mountains larger than ∼30 km in horizontal scale, particularly in low-resolution simulations.The approximate drag slightly overestimates the exact drag, where the error is largest in the lower curves (N∕f = 10).Hence, the error analysis is done for a wide range of Fr −1 n and Ro −1 n .Figure 10 compares the variation of the exact drag (Equation 9; Figure 10a) for flow over a single bell-shaped mountain with the approximate drag (Equation 20; Figure 10b), which is virtually indistinguishable.The dashed lines represent the corresponding values of log 10 N∕f .The ratio of N and f can be seen as N∕f = Ro n ∕Fr n = Bu, where Bu is the (horizontal) Burger number.Drag exists if and only if Bu > 1 (Equation 6); buoyancy plays more of a role than rotation in generating wave drag.The absolute and relative errors are displayed in Figure 10c,d.If 1 < Bu ≲ 10 2 , the drag tends to be overestimated.The maximum absolute error, slightly larger than 0.050, occurs when the inverse Rossby number is slightly below 1 and the inverse Froude number is slightly above 1.However, this is not a vital issue, as the relevant drag values are of order unity.Additionally, the relative error is lower than 5% in most regions of the parameter space where the drag takes appreciable values (Bu ≳ 10 2 ), which is a good precision.
The drag from the overlap of multiple mountains is more sensitive to  −  (not shown).Figure 11 shows the same error analysis as Figure 10; however, it is different for the overlapping term when nonhydrostatic and rotating effects are included, where the overlapping term is normalized by 8 0 N 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 to compare the effect of approximation in the integrand in Equation ( 20).The exact drag (Figure 11a) and the approximate drag (Figure 11b) are virtually indistinguishable in most regions.The absolute and relative errors are displayed in Figure 11c,d.A maximum of the absolute error slightly larger than 0.020 occurs for log 10 Ro −1 i a little below 0.5 or near log 10 Fr −1 i = 0.5.This error is of minor consequence, because the relevant drag value is much smaller than the peak term.In addition, the relative error is lower than 5% in most regions of the parameter space where the drag takes appreciable values, which is a good precision.The relative error diverges near log 10 Fr −1 i = 0.5 or the high Ro −1 i region because the exact drag vanishes while the absolute error is nonzero; nevertheless, it is ∼ 0.02.This also has a minor consequence, because, near log 10 Fr −1 i = 0.5, the relative error becomes lower than 5% if considering drag by both peak and overlap, and for the high Ro −1 i region the drag is small enough to be neglected.

SUMMARY AND CONCLUSIONS
Mountain drag is determined by the model orography and the associated pressure perturbation.This study derived an analytical expression for the orographic gravity-wave drag exerted on randomly distributed multiple three-dimensional bell-shaped mountains in the Boussinesq and linearized framework.The pressure perturbation is obtained using a wave solution to equations governing an inviscid atmosphere with the WKB approximation.We have described the impacts of various effects, such as anisotropic effects, wind-shear effects, Coriolis effects, and nonhydrostatic effects.The total force differs from the sum of the individual forces by each mountain when multiple mountains exist.Specifically, this discrepancy is due to the drag resulting from the overlap of every pair of mountains, D ovlp (Equations 11,A1).D ovlp is influenced by the degree of overlap between every two mountains and the incoming wind direction (Figure 3).D ovlp becomes negligible when the overlap is minimal.In three dimensions, a transverse force redistributes the momentum between the x-and y-directions (Equation 13, Figure 3) if the direction of the incoming air and the alignment of the mountains are neither parallel nor perpendicular.
Hydrostatic, nonrotating drag, in the absence of wind shear, is directly proportional to the horizontal scale of the terrain, denoted as a.The effects of wind shear lead to mountain-wave dispersion and a reduction in the pressure perturbation, which tends to weaken the OGWD (Figures 2e-h, 4).This weakening is controlled by both the Richardson numbers and the incoming wind direction (Equations 14,15,A1).When accounting for nonhydrostatic or rotating effects, the drag decreases drastically with an increase in the Froude number or the inverse Rossby number (Figures 5,6a, 8a,c).Moreover, on the f -plane, additional drag and transverse forces can arise due to wind shear even without mountain anisotropy (Figure 6).This results from the antisymmetric term in the first-order WKB approximation of the pressure perturbation (Figure 7), which occurs if and only if perpendicular wind shear is present and f ≠ 0. Subsequently, we derived the drag for a nonhydrostatic, rotating flow and compared the derived approximate expressions (Equations 16-20) with the exact drag (Figures 5a,10,11).We found that the derived and exact drags are very similar, specifically for large Burger number.The derived expressions indicate that the range of a is limited under the hydrostatic, nonrotating approximation.This limitation underpins why GCMs, such as the ECMWF model, address OGWD only at horizontal scales larger than a few km, typically excluding scales under 5 km as nonhydrostatic (ECMWF, 2021).Our research suggests that the ECMWF's OGWD scheme potentially overestimates mountain drag for scales over 30 km due to the unaccounted-for rotating effect, especially in low-resolution simulations.Further improvement and exploration are necessary.First, we could not obtain an analytical solution to OGWD for multiple three-dimensional elliptic mountains because the integral involving the product of Bessel functions and nonradially symmetric functions was too challenging to handle.Second, our study is based on linear theory, which performs effectively under conditions with a low vertical Froude number; thus, nonlinear mountain drag, flow blocking, and wave breaking over multiple orography require further investigation (e.g., Garner, 1995;Lilly & Klemp, 1979).Third, the mountain drag associated with critical levels, trapped lee waves, and rotors caused by multiple bell-shaped mountains emerges as an important topic for future research (e.g., Teixeira et al., 2013).Fourth, the phenomenon of mountain waves returning to the surface in a nonhydrostatic and rotating atmosphere with wind shear also requires exploration.Finally, this study may offer some insights into the linearized theory of orographic gravity-wave-planetary boundary layer interactions and orographic precipitation over multiple bell-shaped mountains, which needs to be verified.
Current OGWD schemes in general circulation models or numerical weather prediction models do not consider nonhydrostatic or rotating effects in full but truncate them with a sharp length-scale cut-off.This work proposes a more advanced orographic drag scheme that accounts for nonhydrostatic, rotating flow.It allows for a more dynamic procedure for truncating admissible length scales, since the scales being removed depend on the flow itself.It also incorporates the effects of mountain overlap and vertical shear of the horizontal wind.The next step is to apply the new scheme to a general circulation model and examine the impacts on the global simulation.Moreover, to develop a new OGWD scheme based on this study, one should either decompose the high-resolution topography into several bell-shaped mountains or adjust our expressions in terms of the standard deviation of subgrid orography terrain height.

APPENDIX A. WAVE DRAG OVER AN ARBI-TRARY NUMBER OF BELL-SHAPED MOUN-TAINS
The wave drag over multiple bell-shaped mountains is given.If # = 2 and there is no wind shear, the results in the main text are derived.For hydrostatic nonrotating flow, the analytical drag up to second order is given by The notations are the following: ) where  ′ = U ′ k + V ′ l.Then, we differentiate Equation (B4) with respect to z: where  ′′ = U ′′ k + V ′′ l.Finally, we substitute Equation (B5) into Equation (B1) to eliminate p.Then, the final equation, an expression in terms of ŵ, is identical to Equations ( 2)-(3).

APPENDIX C. ANALYTICAL EXPRESSION FOR PRESSURE PERTURBATION
As shown in Equation ( 8), the Fourier transform of the pressure perturbation, p, is a function of the Fourier transform of the vertical velocity perturbation, ŵ, and its derivative.ŵ can be obtained by solving Equations ( 2)-( 5).The WKB solution to Equations ( 2)-( 5) in a scaled vertical coordinate Z = z is In this equation, the vertical wavenumber of the gravity waves m(Z) has been expanded as a power series of .The subscripts in this series denote the order of each respective term.Considering only up to second-order term, the expressions for m 0 , m 1 , and m 2 can be determined when Equation (C1) is introduced into Equation (2) within a scaled coordinate: ) Definition for  1 ,  2 ,  3 , and  4 are provided in Equation (3).Then, by substituting the solution (Equations C1, C2) into Equation 8, p is expressed as a power series of , with the coefficients corresponding to each order given as follows: Here k, F = ∕N, and R = ∕f .Note that U ′′ = V ′′ = 0 is required, that is,  ′′ =  ′ = 0 in Equation (C5) when the rotating effect is considered.One can verify that Equations (C3)-(C5) simplify to eqs.36-38 in Teixeira et al. (2004) in the case of nonrotating flow, namely when f = R −1 = 0 in Equations ( C3)-(C5).Additionally, it is important to note that in a nonrotating flow p0 and p2 are purely imaginary, while p1 is purely real.
Assuming hydrostatic, nonrotating flow, Equations (C3)-(C5) reduce to p0 (k, l) = i 0 N ĥ  ||k|| , (C6) Therefore, when the flow and the mountain function are specified, the surface pressure perturbation in position space can be calculated.Analytical expressions can be derived regardless of the incoming wind's direction, shear, and curvature.For simplicity, however, we present only the zonal wind case here, that is, U(z) = (U 0 + z, 0).For a single bell-shaped mountain, defined as h(x, y) = h 1 ∕(1 + x 2 + y 2 ) 3∕2 , the surface pressure perturbation is obtained as follows by applying the inverse Fourier transform to Equations (C6)-(C8): p 0,#=1 (x, y) = − 0 NU 0 h 1 x (1 + x 2 + y 2 ) 3∕2 , (C9) In Equations ( C10)-(C11), B1 , B2 , B3 , and B4 are defined as The pressure field resulting from multiple bell-shaped mountains (Equation 10) is simply the sum of the pressure fields resulting from each individual bell-shaped mountain.The pressure field generated by a single bell-shaped mountain is given by Equations (C9)-(C11).However, adjustments need to be made to the position of each mountain's center and its horizontal width, which are straightforward modifications.Equations (C9)-(C12) are plotted in Figure 2.
illustrates the relative size of the transverse force compared with the magnitude of D peak .Both forces vary with d.D and T reach their extreme values at d = 4∕3 ≈ 1.33 and d = √6∕3 ≈ 0.82, respectively.If the overlap is too small (e.g., d ≥ 3), D ovlp is small enough to be neglected.If the two mountains overlap completely, this can be seen as a single mountain with double height.This happens when d = 0, D ovlp = 2D peak : then the total drag becomes 4D peak , equivalently, D∕||D peak || = 4 and T = 0.If a moderate overlap occurs, the role of D ovlp becomes important to make the total mountain drag

F
I G U R E 3 Normalized (a) drag and (b) transverse force due to two identical bell-shaped mountains.Both forces are normalized by D 0 =  4  0 N 0 ah 2 , the magnitude of the drag force over a single mountain.

F
Normalized zonal drag as a function of the inverse Richardson number of x and y directions for a hydrostatic, nonrotating flow as given by U = U 0 + z, V = V 0 + z.Double mountains are aligned on the x-axis, and d = 1 for simplicity.The upper panels are the normalized drag due to a single peak, and the lower panels are the normalized drag due to the overlap of double hills.The normalizing factors are D 0 =  4  0 NU 0 ah 2 and 4 0 NU 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 for the upper and lower panels, respectively.(a,d) V 0 ∕U 0 = −1; (b,e) V 0 = 0; (c,f) V 0 ∕U 0 = 1.

F
I G U R E 5 (a) Normalized drag over a single bell-shaped mountain for nonrotating flow as a function of Fr, where the normalizing factor is D 0 =  4  0 N 0 ah 2 .(b) Normalized drag and (c) transverse force due to D ovlp .d = 1, where the normalizing factor is 8

F
I G U R E 6 Normalized (a) drag D peak and (b) transverse force T peak for hydrostatic flow as a function of Ro −1 , where the normalizing factor is D 0 =  4  0 N 0 ah 2 .(a) Solid curve: drag when no wind shear exists.Green region: drag in the presence of wind shear, where Ri = 1.The upper green curve corresponds to  −  ′ = −∕2, and the lower green curve corresponds to  −  ′ = ∕2.(b) Solid curve (T peak = 0): transverse force when no wind shear exists.Green region: transverse force in the presence of wind shear, where Ri = 1.The upper green curve corresponds to  −  ′ = , and the lower green curve corresponds to  −  ′ = 0.

F
Normalized pressure perturbation at the surface for sheared flow on an f -plane as given by U = U 0 + z, V = z when Ro = 1 (thick solid circle).Terrain elevation equal to h max ∕2; (a,b) zeroth-order pressure perturbation; (c,d) difference between the first-order pressure when a n = 10 4 and that from Figure 2c,d.

Figure
Figure 8a,c shows how the normalized drag and transverse force according to D ovlp change in Ro −1 i and  − , respectively.As in Figure5b,c, a moderate overlap is chosen (d = 1), and they are divided by 8 0 N 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 .D ovlp and T ovlp are periodic with period  as in the hydrostatic, nonrotating case (Figure3).However, the reflection symmetry of T ovlp is broken if Ro −1 i > 0. The  −  at which T ovlp has extreme values recedes from ∕4 + n∕2 as Ro −1 i increases.The sign still changes frequently depending on  −  when the rotating effect is strong, which is in contrast to the sign changing less depending on  −  when the nonhydrostatic effect is strong (Figure5b,c).The magnitude of normalized D ovlp becomes small if Ro −1 i ≳ 1.This is because the Coriolis force suppresses the wave drag when the horizontal length scale is large enough.Considering that normalized D peak is also small when Ro −1 n ≫ 1 (Figure6a), we calculated the ratio of the magnitude of D ovlp and T ovlp to the magnitude of the zeroth order of D peak , and found that D peak converges much faster than D ovlp or T ovlp in the high Ro −1 n limit.Figure8b,d displays the normalized drag and transverse force up to the first order; hence, it depicts the wind-shear effects.The degree of overlap and the

F
I G U R E 8 Normalized drag and transverse force due to D ovlp .d = 1, where the normalizing factor is 8 0 N 0 h 1 h 2 a 2 1 a 2 2 ∕(a 1 + a 2 ) 3 .(a) Normalized drag force up to zeroth order; (b) normalized drag force up to first order ( −  ′ = ∕4); (c) normalized transverse force up to zeroth order; (d) normalized transverse force up to first order ( −  ′ = ∕4).

F
I G U R E 10 Normalized drag as a function of Fr −1 and Ro −1 for a bell-shaped mountain and constant wind.Solid lines with color: labeled contours; dashed lines values of log 10 N∕f .(a) Exact drag; (b) approximate drag; (c) absolute error; (d) relative error.

F
I G U R E 11 Same as Figure 10 but for the overlap of bell-shaped mountains.(a) Exact drag; (b) approximate drag; (c) absolute error; (d) relative error.