An Empirical Model for Non-Linear Pressure Drag Across Non-Hydrostatic Flow Regimes with Trapped Lee Waves

1. Introduction

The hydrostatic approximation has been adopted in many classic studies because mesoscale non-hydrostatic effects are often modest, and this assumption greatly simplifies the analysis. Consequently, most modern orographic drag parameterization schemes were developed for hydrostatic flows, with closed-form expressions for total drag derived in several foundational studies [1,2,3,4].

Since many flows over mountains are nonlinear, studying the impact of nonlinearity on gravity-wave drag is essential. In the hydrostatic regime, this impact is well-documented: Lilly and Klemp [5] showed how finite amplitude and asymmetry amplify drag, while Smith [6] and Durran [7] clarified the limits of linear theory and scaling laws. Furthermore, Ólafsson and Bougeault [8] and Miranda and James [9] examined drag enhancement and the transition to wave-breaking for non-dimensional mountain heights $N{h}_0/U\sim O\left(1\right)$, where drag substantially surpasses linear values.

On the other hand, as wind speed increases or static stability decreases, flow becomes increasingly non-hydrostatic. Unlike the hydrostatic regime, where all forced waves propagate vertically, non-hydrostatic scenarios render high-wavenumber components evanescent, halting their upward momentum transport while longer waves continue to contribute to drag. If an evanescent layer overlies a propagating one, wave trapping and reflection occur, leading to resonance and amplified drag—phenomena also possible via partial reflection under hydrostatic conditions [10,11]. Notably, non-hydrostatic drag can match or exceed its hydrostatic counterpart [12,13].

In contrast to hydrostatic flows, nonlinear drag in non-hydrostatic flows has received considerably less attention. Among the existing literature, Peltier and Clark [14] documented amplification regimes and transitions to wave breaking, while Lott [15] analysed how nonlinear intensification of TLW modes and critical-level interactions modify momentum flux. Vosper [16] demonstrated that linear theory underestimates wave amplitudes in boundary-layer inversions for short horizontal wavelengths; large amplitudes can trigger flow separation and rotors, leading to significant regime shifts in drag. Doyle and Coauthors [17], using T-REX campaign simulations, reported how wave intensification, TLW, and rotor patterns modulate pressure drag. Moreover, Teixeira et al. [12] used a two-layer atmosphere with piecewise-constant parameters following [18] to show that TLW drag can equal or exceed that of vertically propagating waves, with significant implications for drag parameterisations [19]. Similarly, Teixeira et al. [13] found comparable contributions from TLW and propagating waves for lee waves trapped at an inversion capping a neutral layer. Despite these foundational studies, a comprehensive mapping of drag as a function of nonlinearity in non-hydrostatic flows featuring TLWs is still lacking. Recently, Argain [20] verified that nonlinearity is comparably important for the drag produced by both upward-propagating waves and TLWs. Therefore, further studies are essential to better understand this drag behaviour and to establish a robust basis for more accurate parameterisation schemes in weather forecasting models.

For the reasons explained above, the primary motivation of this work is to develop baseline, easily implementable empirical model that can serve as a base tool for building viable alternatives to parameterize total pressure drag under nonlinear conditions. Based on the foundational two-layer configuration of Scorer [18], the proposed model is evaluated using diverse TLW cases from the literature, including aircraft observations over the Blue Ridge Mountains [21], satellite imagery of Macquarie Island [22], and wind profiler data from Hong Kong [23]. Additionally, reference profiles derived from high-resolution simulations and reanalyses [24] are used to represent upper-level evanescence regimes. The model is shown to agree well with numerical simulations for moderate to strong non-hydrostatic and nonlinear regimes featuring TLW generation.

"The article is organized as follows: Section 2 and Section 3 describe the linear and numerical models, respectively. Section 4 introduces the empirical model, its theoretical basis, and its practical implementation procedure. Section 5 provides the model validation and discussion, and finally, Section 6 offers concluding remarks."

2. Linear Model

We consider a linear model for a steady, inviscid, non-rotating, and stratified flow over a 2D mountain ridge of relatively small amplitude, oriented perpendicularly to the incident flow. The spatial scale of the flow is large enough that viscous effects are negligible (i.e., a free-atmosphere flow rather than a boundary-layer flow), yet sufficiently small that the Earth’s rotation can be safely ignored (implying a large Rossby number, $Ro\gg 1$).

By linearising the equations of motion under the Boussinesq approximation about a reference mean state and expressing the dependent variables as Fourier integrals in x, the Fourier transform of the vertical-velocity perturbation, $w_k$, is shown to satisfy the Taylor–Goldstein equation:

$${\displaystyle \frac{{\partial }^2{w}_k}{\partial z^2}}+\left(l^2-k^2\right)w_k=0,\mathrm{with}l={\left({\displaystyle \frac{N^2}{U^2}}-{\displaystyle \frac{1}{U}}{\displaystyle \frac{d^{2}U}{d{z}^2}}\right)}^{1/2},$$

(1)

where l is the Scorer parameter of the atmosphere. Here, k is the horizontal wavenumber in the x direction, $N^2\left(z\right)>0$ is the static stability of the reference state, and $U\left(z\right)$ is the incoming wind speed perpendicular to the ridge. Equation (1) describes how the amplitude of the vertical velocity perturbation for each wavenumber varies with altitude, depending on the vertical structure of the atmosphere through $l\left(z\right)$. The existence of TLWs fundamentally requires that the squared Scorer parameter, $l^2\left(z\right)$, decreases with height sufficiently for high-wavenumber components to become evanescent aloft while remaining vertically propagating below. Specifically, there must exist wavenumbers k satisfying $l_2<k<l_1$. This condition is a strict prerequisite for the formation of trapped modes within the lower tropospheric waveguide.

In this model, the atmosphere is divided into vertical layers, each with a constant Scorer parameter. Within each layer, the vertical velocity solution takes a simple mathematical form: either a combination of oscillatory functions (if the flow is propagating) or exponential functions (if it is evanescent), depending on the relationship between the horizontal wavenumber and the local Scorer parameter.

Boundary conditions are imposed to determine the arbitrary constants in each layer. At the surface ($z=0$), $w_k$ is determined by the orography ($h_k$), as the incident flow is vertically displaced by the underlying topography. At the interfaces between layers, continuity of vertical velocity and pressure (or its vertical derivative) is required to ensure smooth transitions. Finally, a radiation condition is applied at the top of the uppermost layer to ensure that waves propagate upward or decay without artificial reflection. Once the system is solved for each wavenumber, the complex amplitude $w_k(k,z)$ is determined for any altitude.

From the linearised horizontal momentum equation, the pressure perturbation $p_k$ can be expressed in terms of the vertical velocity and its vertical derivative:

$$p_k={\displaystyle \frac{{\rho }_0}{ik}}\left(U\left(z\right){\displaystyle \frac{\partial w_k}{\partial z}}+w_k{\displaystyle \frac{\partial U}{\partial z}}\right).$$

(2)

Evaluating this expression at $z=0$ yields the amplitude of the surface pressure perturbation for each wavenumber.

The spatial distribution of the surface pressure, $p(x,z=0)$, is reconstructed via an inverse Fourier transform. Finally, the pressure drag force across the ridge, per unit length in the span-wise direction, is defined as:

$$D={\int }_{-\infty }^{+\infty }p(x,z=0){\displaystyle \frac{\partial h}{\partial x}}dx.$$

(3)

The model described above (hereafter the Drag Linear Model, DLM) is primarily used to compute the resonance modes, characterised by their wavenumber $k_{tlw}$ (or wavelength ${\lambda }_{tlw}=2\pi /k_{tlw}$), and the total linear pressure drag, $D_L$. Additionally, the linear theoretical models of Teixeira et al. [12] and Teixeira et al. [13] are employed. Developed under the same assumptions, these models are applied to idealised atmospheric configurations and allow D to be partitioned into TLW drag ($D_{tlw}$) and vertically propagating wave drag ($D_{pw}$): $D=D_{pw} + D_{tlw}$.

Finally, the orography is assumed to be a symmetric, bell-shaped ridge:

$$h\left(x\right)={\displaystyle \frac{h_0}{1+{(x/a)}^2}}\Rightarrow h_k\left(k\right)={\displaystyle \frac{h_{0}a}{2}}{e}^{-a\left|k\right|},$$

(4)

where $h_0$ is the maximum height and a is the half-width.

3. Numerical Model

This study employs the FLEX model, a 2D micro-to-mesoscale numerical model in orthogonal curvilinear coordinates, previously used to investigate resonant mountain waves [11,12,13,25,26].

The vertical spacing is $\Delta z=25$ m. The horizontal resolution $\Delta x$ is the minimum of ${\lambda }_{\mathrm{tlw},min}/40$ and $a/5$, ensuring at least 40 points per minimum TLW wavelength and 5 points per mountain half-width. The computational domain extends vertically to a height of six hydrostatic vertical wavelengths (${\lambda }_{z0}=2\pi /l_1$). Horizontally, it spans $7a$ upstream and $15{\lambda }_{tlw,max}$ downstream. Here, ${\lambda }_{tlw,min}$ and ${\lambda }_{tlw,max}$ denote the shortest and longest wavelengths, respectively, for configurations where multiple TLW modes coexist; in single-mode cases, both default to the unique wavelength, ${\lambda }_{tlw}$. The integration time step is set to $\Delta t=1$ s.

A $2.5{\lambda }_{z0}$ top sponge layer is applied, while a radiation condition and lateral sponges are used at the boundaries. FLEX is calibrated using $D_L$ values from the DLM, as the total drag D is sensitive to boundary absorbing sponges. Accurate $D_L$ calculation is vital, as it is used to normalize D in the empirical model.

4. Empirical Model

Figure 1 illustrates the various atmospheric configurations considered for the development of the empirical model. Figure 1a depicts the two-layer Scorer atmosphere which, as previously noted, forms the basis of this empirical model (hereafter DEM, for Drag Empirical Model). In this baseline configuration, the DLM and FLEX models employ a uniform background wind of $U=10$ ms⁻¹ (so the curvature term vanishes). Static stability is prescribed as $N_1=0.02$ s⁻¹ in the lower layer and $N_2=0.004$ s⁻¹ in the upper layer, resulting in $N_2/N_1=0.2$ and $l_2/l_1=0.2$.

For the two-layer Scorer profile (Figure 1a), the dimensionless mountain height and half-width are determined as ${\widehat{h}}_0=l_1{h}_0$ and $\widehat{a}=l_{1}a$. For more complex profiles, such as those illustrated in Figure 1b–e, ${\widehat{h}}_0$ and $\widehat{a}$ are defined using the maximum value of $l\left(z\right)$, closest to the surface, that produces TLW. For instance, in the case of the profile in Figure 1c, ${\widehat{h}}_0=l_H{h}_0$ and $\widehat{a}=l_{H}a$, while for Figure 1d, ${\widehat{h}}_0=l_{H_1}{h}_0$ and $\widehat{a}=l_{H_1}a$.

In all cases, a bell-shaped mountain (4) is used to control a and $h_0$ variations. Nonlinearity is modulated via $h_0$, from the linear regime (${\widehat{h}}_{0,L}=0.02$) up to ${\widehat{h}}_0\approx 1$ (e.g., $h_0$ ranging from ∼10 m to ∼500 m, for typical $l_1$ values. Larger ${\widehat{h}}_0$ values are omitted because, as shown by Argain [20], TLW effects on drag variability—the dominant driver in this flow regime—become secondary for ${\widehat{h}}_0>1$. Consistent with that study, no wave breaking was observed in our simulations. Indeed, when atmospheric conditions strongly favor wave trapping in the lower layers, vertical energy propagation is severely restricted, effectively suppressing classic mid-tropospheric wave breaking even at higher non-linearities (e.g., Doyle and Reynolds [27]). Therefore, the associated drag enhancement typical of breaking waves is absent.

4.1. DEM for a Two-Layer Scorer Atmosphere

Figure 2 shows the normalised linear drag $D_L/D_{L,0}$ and its components ($D_{L,tlw}$ and $D_{L,pw}$) as a function of $\widehat{H}\in [0.5,3.0]$, calculated following Teixeira et al. [12] for a Scorer atmosphere (Figure 1a). Results are presented for moderate ($\widehat{a}=5,a=2500$ m; Figure 2a) and strong ($\widehat{a}=2,a=1000$ m; Figure 2b) non-hydrostatic effects, both in the linear regime (${\widehat{h}}_{0,L}=0.02$). Here, $D_{L,0}=0.25\pi {\rho }_0{l}_1{\left(h_{0,L}U\right)}^2$ is the linear hydrostatic reference.

Notably, in non-hydrostatic regimes (Figure 2b), $D_{L,tlw}$ can equal or exceed $D_{L,pw}$, representing a substantial fraction of the total drag. In Figure 2a (moderately non-hydrostatic), $D_{L,tlw}$ remains below $D_{L,pw}$ but still accounts for $\approx 50\%$ of the total drag. This contrasts with classical hydrostatic cases where the TLW contribution is negligible ([12] e.g., Figure 17a). Generally, flows where $D_{L,tlw}\sim D_L$ exhibit moderate to strong non-hydrostaticity ($\widehat{a}\approx 1$–5), which is the focus of this study. A key insight for the DEM, derived from Figure 2a,b, is the similarity in the functional form of drag across different trapped modes.

This insight inspired the formulation of a generic drag function, $D\left(\widehat{H}\right)$, calibrated for the first mode ($n=1$) within $\widehat{H}\in [0.5,1.5]$, which is then scaled to reconstruct the behavior of higher-order modes ($n>1$).

Figure 3a illustrates $D/D_L\left({\widehat{h}}_0\right)$ profiles for $\widehat{H}=0.54$, $0.7$, and $0.9$ (within the $n=1$ range, $\widehat{H}\in [0.5,1.5]$). Figure 3b presents analogous profiles for different modes: $\widehat{H}=0.8$ ($n=1$), $\widehat{H}=1.8$ ($n=2$), and $\widehat{H}=2.8$ ($n=3$). Furthermore, for comparison, both plots include a linear extrapolation of the hydrostatic single-layer total drag, $D_L=0.25\pi {\rho }_0{l}_1{\left(h_{0}U\right)}^2$, normalized by the linear value $D_{L,0}$.

Consistent with Argain [20], the $n=1$ profiles (Figure 3a) exhibit a quadratic behaviour ($D/D_L\propto {\widehat{h}}_0^2$), suggesting that the linear scaling is preserved and can be fitted by an expression of the form $\Phi \left(\tilde{H}\right){\tilde{h}}_0^2+1$. In Figure 3b, the $n=1$ and $n=2$ profiles coincide up to ${\widehat{h}}_0\approx 0.5$ before $n=2$ diverges toward higher drag. A similar trend occurs for $n=3$, with divergence starting earlier at ${\widehat{h}}_0\approx 0.3$. These results indicate that both the quadratic exponent and $\Phi$ must be adjusted (increased) as the mode number n increases. In this formulation, this is achieved by multiplying $\Phi \left(\tilde{H}\right)$ by a function of the form ${\left[{\tilde{h}}_0^{\beta }\Psi \left(\tilde{H}\right)\right]}^{n-1}$.

Based on a detailed analysis of all numerical profiles, the following expression is proposed to estimate D for the two-layer Scorer model (using the $l\left(z\right)$ profile shown in Figure 1a), for arbitrary values of n and $\widehat{H}$:xxx

$${\displaystyle \frac{D}{D_L}}\left({\tilde{h}}_0\right)=\Phi \left(\tilde{H}\right){\tilde{h}}_0^2{\left[{\tilde{h}}_0^{1/5}\Psi \left(\tilde{H}\right)\right]}^{n-1}+1,$$

(5)

where ${\tilde{h}}_0={\widehat{h}}_0/{\widehat{h}}_{0,1}-1$ and $\tilde{H}=\widehat{H}/0.5-1$. In this formulation, ${\widehat{h}}_{0,1}$ represents the minimum dimensionless mountain height considered in the numerical simulations. (It should be noted that the fractional exponent 1/5 is empirically derived to optimally capture the drag amplification). Furthermore, $\widehat{H}$ is adjusted such that $\widehat{H}=\widehat{H}+0.5$ if $\widehat{H}<0.5$, and $\widehat{H}=\widehat{H}-\mathrm{nint}\left(\widehat{H}\right)+1$, where $\mathrm{nint}$ denotes the nearest integer function. The parameters $D_L$ and n are determined using the DLM framework.

The functions $\Phi \left(\tilde{H}\right)$ and $\Psi \left(\tilde{H}\right)$ are defined by the following expressions:

$$\Phi \left(\tilde{H}\right)=a_{1}exp[+a_2{\tilde{H}}^{a_3}]+a_4{\tilde{H}}^{a_5}+a_6\mathrm{for}\widehat{H}\le 1.032$$

(6)

$$\Phi \left(\tilde{H}\right)=b_{1}exp[-b_2{\tilde{H}}^{b_3}]+b_4{\tilde{H}}^{b_5}+b_6\mathrm{for}\widehat{H}>1.032$$

(7)

$$\Psi \left(\tilde{H}\right)=c_1+{\displaystyle \frac{c_2-c_1}{1+exp\left[c_3(\tilde{H}-c_4)\right]}}\mathrm{for}\widehat{H}\le 0.825$$

(8)

$$\Psi \left(\tilde{H}\right)=d_2{\tilde{H}}^2+d_1\tilde{H}+d_0\mathrm{for}\widehat{H}>0.825$$

(9)

To determine the empirical coefficients for the functions $\Phi \left(\tilde{H}\right)$ and $\Psi \left(\tilde{H}\right)$, optimisation methods based on the least-squares criterion were applied. Specifically, the polynomial components were fitted using standard linear regression (via MATLAB’s polyfit routine). For the non-linear components, the fit was achieved by minimising a non-linear sum-of-squared-errors cost function. It is important to emphasise that the mathematical form chosen for these non-linear components—a logistic (sigmoid) function—was deliberately selected to guarantee a stable asymptotic behaviour. This formulation ensures that the drag corrections remain bounded, preventing numerical divergences or unrealistic drag amplification if the Drag Empirical Model (DEM) is extrapolated to flow configurations slightly outside the strictly tested parameter domain.

The sets of constants $a_i,b_i,c_i,$ and $d_i$, derived from fitting, are as follows:

a)

$a_1=-0.0007272$, $a_2=7.0221$, $a_3=1.0645$, $a_4=4.61$, $a_5=2.3724$, and $a_6=0.56966$;

b)

$b_1=2.8863$, $b_2=0.031476$, $b_3=19$, $b_4=0.056042$, $b_5=-1.1833$, and $b_6=2.1822$;

c)

$c_1=0.642$, $c_2=1.0026$, $c_3=20.7235$, and $c_4=0.29845$;

d)

$d_2=0.76594$, $d_1=-2.31083$ and $d_0=1.748$.

These functions are specifically designed to capture the regime transitions and modal dependence inherent in moderate to strong non-hydrostatic flows. By distinguishing between different ranges of the dimensionless duct depth $\widehat{H}$, the model accounts for the shifting dominance of trapped lee wave modes. Furthermore, the high precision of the reported constants is essential for the numerical reproducibility of the DEM, ensuring consistent drag estimations across the diverse atmospheric configurations explored in this study.

Modes are limited to $n=3$, as higher-order configurations are rare in nature and increase model complexity without proportional practical gains. Furthermore, numerical simulations show that achieving a quasi-stationary state for D becomes increasingly difficult as n increases, requiring greater computational effort to ensure convergence.

4.2. DEM for an Arbitrary Atmosphere

A key feature of the DEM is the normalization by $D_L$, which ensures $D/D_L=1$ in the linear limit and standardises results across distinct $l\left(z\right)$ structures. As demonstrated below, simulations using various idealized configurations (Figure 1a–d) with identical $\widehat{H}$ show that $D/D_L\left({\widehat{h}}_0\right)$ profiles remain remarkably similar, often nearly coincident.

Figure 4a contrasts the $D/D_L\left({\widehat{h}}_0\right)$ profiles for $\widehat{H}=0.6$ (no symbols) and $\widehat{H}=1.07$ (symbols) to illustrate how the geometry of the $l\left(z\right)$ profile influences drag. For $\widehat{H}=0.6$, the dotted line represents an extreme case where waves are trapped by a thin inversion capping a neutral layer ($l_1=0$ m⁻¹; Figure 1b), a flow regime investigated by Teixeira et al. [13]. This configuration features an inversion at $H=500$ m ($l_H=1.25\times {10}^{-3}$ m⁻¹, $\Delta z_1=\Delta z_2=25$ m) with $l_2=1\times {10}^{-3}$ m⁻¹ aloft.

Notably, the non-normalized Scorer (solid) and thin-inversion (dotted) profiles diverge significantly despite sharing the same $\widehat{H}$.

For $\widehat{H}=1.07$, the dotted line with symbols corresponds to an $l\left(z\right)$ profile of the type shown in Figure 1c with the following parameters: $H=1000$ m, $l_H=2.8\times {10}^{-3}$ m⁻¹, $l_1=1.4\times {10}^{-3}$ m⁻¹, $l_2=2.8\times {10}^{-4}$ m⁻¹, $\Delta z_1=200$ m, and $\Delta z_2=400$ m. The solid line with symbols represents the baseline Scorer profile. A comparison of these two $D\left({\widehat{h}}_0\right)$ profiles reveals that, consistent with the previous case, they differ considerably.

Using the same symbology, Figure 4b shows the same profiles as Figure 4a, but normalized by $D_L$. The high degree of coincidence observed for both $\widehat{H}=0.6$ and $1.07$ confirms that $D_L$-normalization enables the two-layer Scorer model to serve as a universal baseline for the DEM across distinct $l\left(z\right)$ structures with identical $\widehat{H}$.

Additional tests (not shown here) with $l\left(z\right)$ profiles featuring multiple $l_H$ maxima, such as that in Figure 1d, confirm that $D_L$-normalization is equally effective for these configurations. In such cases, the total surface drag $D\left({\widehat{h}}_0\right)$ is effectively the sum of the contributions from each maximum. Thus, following Eq. (5), the total drag is:

$${\displaystyle \frac{D}{D_L}}\left({\tilde{h}}_0\right)={\displaystyle \frac{1}{D_L}}{\displaystyle \sum _{i=1}^m}{D}_i({\tilde{H}}_i,{\tilde{h}}_0,n_i),$$

(10)

where m is the number of contributing maxima, while $n_i$ and $D_i$ denote the number of modes and the drag contribution associated with each ${\tilde{H}}_i$, respectively. It should be noted that using the DLM is crucial for determining, through exclusion, which maxima effectively generate TLW.

Equation (10) defines the final formulation of the Drag Empirical Model (DEM), which is, in principle, applicable to any $l\left(z\right)$ profile.

4.3. Final Remarks on the DEM Formulation

Although the DEM intentionally smooths the abrupt drag transitions associated with mode shifts, this simplification is robustly justified. Firstly, the regions of abrupt transition are parametrically narrow. For instance, in the non-hydrostatic linear reference regime (Figure 3b), the transition from $n=1$ to $n=2$ occurs over a small interval of $\Delta \widehat{H}\approx 0.2$ (between $\widehat{H}\approx 1.4$ and $1.6$), and a similar width is observed for the $n=2$ to $n=3$ shift (between $\widehat{H}\approx 2.4$ and $2.6$). Consequently, the statistical error introduced by this smoothing over a broad climatological distribution of atmospheric states is minimal. Secondly, and more importantly from a physical standpoint, the sharp resonant discontinuities predicted by idealized piecewise-constant theory are rarely realized in the actual atmosphere. Natural variability in atmospheric profiles, transient flow effects, and three-dimensional dispersion inherently "smear" these resonant peaks. Therefore, the continuous functions of the DEM arguably provide a more realistic representation of bulk atmospheric drag than a strictly discontinuous formulation.

Building upon this rationale of physical smoothing, a specific geometrical adjustment is introduced for regimes with an extremely shallow lower atmospheric layer, mapping $\widehat{H}$ to $\widehat{H}+0.5$ when $\widehat{H}<0.5$. This adjustment serves as a physical and numerical regularisation anchored in the behaviour of the linear reference regime. As established for the fundamental mode ($n=1$), the non-linear drag largely inherits the linear scaling ($D/D_L\propto {\widehat{h}}_0^2$); hence, the DEM anchoring for $n=1$ is guided by the reference behaviour shown in Figure 2a,b. For shallow lower layers ($\widehat{H}<0.5$), the flow exhibits one of two distinct states. In extremely shallow cases (e.g., $\widehat{H}\lesssim 0.4$), the atmosphere lacks the physical waveguide thickness required to sustain efficient wave trapping, rendering TLW effects negligible. Conversely, as $\widehat{H}$ approaches $0.5$, the flow undergoes a sharp transition from $n=0$ to $n=1$ over a narrow interval ($\Delta \widehat{H}\approx 0.2$, between $\approx 0.4$ and $0.55$). As previously discussed, the DEM intentionally smooths these abrupt mode-shift transitions. Therefore, the $+0.5$ shift effectively maps this entire shallow domain—encompassing both the near-zero drag region and the narrow resonant onset—directly into the fundamental pre-resonant interval $\widehat{H}\in [0.5,1.0]$. This transformation ensures the model scales predictably within a well-calibrated zone, remaining safely clear of narrow, abrupt mode-switching transitions that could induce spurious numerical resonances. Consequently, the continuous functions of the DEM successfully yield stable, bounded, and physically realistic drag predictions that accurately reflect the background flow as $\widehat{H}\to 0$.

4.4. Practical Implementation Procedure of the DEM

To facilitate the practical application of the proposed model, the following clear and concise workflow is recommended for its implementation:

1.

DLM Profile Representation: Use the DLM to represent the vertical profile of the Scorer parameter, $l^2\left(z\right)$, through a discretization into constant-value layers. This step establishes the theoretical basis for identifying wave modes.

2.

Filtering and Mode Selection: Based on the identified maxima of $l^2\left(z\right)$ in the profile, apply an exclusion criterion to determine which layers effectively support trapped lee waves. Identify only the maxima that generate genuine modes, thereby defining the number of interfaces, their locations (${\widehat{H}}_i$), and the respective horizontal wavenumbers ($k_n$) and eigenvalues (${\lambda }_n$).

3.

Numerical Model Calibration and Validation: Before final DEM application, use the numerical model (FLEX) to compute the total drag (D) under reference linear conditions. We suggest that this calibration be performed using the linear case ($D_L$), as this is typically the only robust reference value available. However, if estimated drag values for non-linear situations are available, derived from local observational data or site-specific campaigns, the calibration process will be significantly more robust and should be prioritized. During this step, adjust the sponge layer characteristics (including the vertical and lateral damping coefficients and the respective layer thicknesses) and the domain dimensions (including grid resolution, horizontal and vertical domain extent) to ensure the computed drag is physically consistent and free from spurious reflections.

4.

Empirical Model (DEM) Application: With the validated modes and the determined drag D, apply the final DEM formulation. This model allows for a robust estimation of non-linear pressure drag, integrating the contribution of the identified modes while avoiding the need for computationally expensive numerical integrations.

5. Model Validation and Discussion

To establish a real-world benchmark for validation, atmospheric profiles derived from observational TLW field campaigns are utilised to compare the $D/D_L\left({\widehat{h}}_0\right)$ results computed by the FLEX model against those predicted by the DEM. To facilitate ${\widehat{h}}_0$ variation, a bell-shaped mountain (4) replaces actual orography. Profiles and DLM parameters (Figure 5) generally exhibit moderate non-hydrostatic effects, except for the strongly non-hydrostatic case in Volkert and collaborators [24]. All data were digitized from the original publications.

Sachsperger et al. [28] analyzed lee waves trapped at a boundary-layer inversion. Applying the DLM to the $l\left(z\right)$ profile reveals a single $l_H$ maximum and one mode ($m=1,n=1$) with: ${\lambda }_{tlw}=4974$ m, $l_H=4.87\times {10}^{-3}$ m⁻¹, $H=1153$ m, $\widehat{H}=1.79$, $a=1000$ m, and $\widehat{a}=4.9$ (moderate non-hydrostaticity).

Chan [23] studied stationary lee waves over Hong Kong ($\lambda \approx 2.8$ km). While the author suggested a single mode, the DLM identifies two TLW modes ($m=1,n=2$) from a single $l_H$ maximum. Parameters include: ${\lambda }_{tlw1}=2740$ m, ${\lambda }_{tlw2}=942$ m, $l_H=8.9\times {10}^{-3}$ m⁻¹, $H=750$ m, $\widehat{H}=2.12$, $a=450$ m, and $\widehat{a}=4$ (moderate non-hydrostaticity).

Smith [21] reported aircraft observations of lee waves over the Blue Ridge Mountains, where vertical trapping was driven by an $l\left(z\right)$ profile decaying aloft. Similar to Sachsperger et al. [28], this case features $m=1,n=1$. Parameters: ${\lambda }_{tlw}=7450$ m, $l_H=3.9\times {10}^{-3}$ m⁻¹, $H=2520$ m, $\widehat{H}=3.13$, $a=1000$ m, $\widehat{a}=3.9$ (moderate non-hydrostaticity).

Volkert and collaborators [24] presented a phase diagram distinguishing between trapped lee waves (TLWs) and vertically propagating mountain waves based on the $l\left(z\right)$ profile. Unlike the direct observational data from other authors discussed herein, the profiles in Volkert and collaborators [24] represent reference cases derived from high-resolution numerical simulations (such as MESO-NH) and meteorological reanalyses. For the specific reference case selected from Volkert and collaborators [24], similarly to the case of Chan [23], the DLM reveals two TLW modes associated with a single $l_H$ maximum ($m=1$, $n=2$). The corresponding flow and terrain parameters are $l_H=7\times {10}^{-4}$,m⁻¹, ${\lambda }_{\mathrm{tlw},1}=14615$,m, ${\lambda }_{\mathrm{tlw},2}=8820$,m, $H=7609$,m, $\widehat{H}=1.69$, $a=2000$,m, and $\widehat{a}=1.4$ (strong non-hydrostaticity).

Finally, Mitchell et al. [22] presented satellite observations of TLW over Macquarie Island. We utilize profiles C and E; both contain local $l\left(z\right)$ maxima that do not produce waves ($7.7\times {10}^{-4}$ m⁻¹ at 4600 m for C, and $1.1\times {10}^{-3}$ m⁻¹ at 4000 m for E). Profile E reveals one mode ($m=1,n=1$) with: $l_H=1.52\times {10}^{-3}$ m⁻¹, ${\lambda }_{tlw}=8650$ m, $H=1095$ m, $\widehat{H}=0.53$, $a=1000$ m, $\widehat{a}=1.5$ (strong non-hydrostatic effects). Profile C exhibits two maxima and two modes ($m=2,n_1=1,n_2=1$) with: $l_{H_1}=1.47\times {10}^{-3}$ m⁻¹, ${\lambda }_{tlw,1}=8920$ m, $H_1=1536$ m, ${\widehat{H}}_1=0.72$, $l_{H_2}=7.93\times {10}^{-4}$ m⁻¹, ${\lambda }_{tlw,2}=13250$ m, $H_2=7104$ m, ${\widehat{H}}_2=1.79$, $a=2000$ m, $\widehat{a}=2.94$ (moderate to strong non-hydrostaticity).

Figure 6 and Figure 7 compare $D/D_L\left({\widehat{h}}_0\right)$ results from FLEX and DEM for the observational cases. Except for [28] and [22] C, both models show high consistency. The [22] C case is the most complex ($m=2,n_1=n_2=1$), with discrepancies appearing only for ${\widehat{h}}_0>0.7$ while maintaining robust overall agreement. For [28], numerical scatter complicates validation, though the DEM follows the general trend of the simulations. Overall, these results confirm the reliability of the DEM across diverse atmospheric conditions, ranging from moderate to severe non-hydrostatic regimes.

Limitations and Applicability of DEM

The DEM assumes a 2D, inviscid, non-rotating flow, where the boundary layer (BL) is either negligible or the waves of interest reside above the friction-dominated region. While the real atmosphere features a BL up to 1 km, an inviscid free-slip lower boundary is physically consistent for scenarios where its influence on wave generation is secondary.

Following the analysis of Argain [20], this inviscid configuration precludes wave breaking and classical rotors, which are intrinsically linked to friction and no-slip conditions. In their absence, the system responds via internal mechanisms: duct resonance detuning, harmonic generation, and leakage into propagating waves. These processes typically lead to drag saturation or redistribution rather than sharp monotonic increases.

However, since trapped lee waves (TLWs) and rotors are often interconnected, neglecting friction removes a key momentum transfer pathway. Consequently, by suppressing friction, the DEM may either overestimate wave drag in rotor-prone regimes, or conversely, underestimate total drag in cases where friction-driven flow separation creates additional form drag (e.g., wake effects). Despite this, the model remains a robust tool for a wide range of realistic atmospheric scenarios. The relevance of such internal-wave dynamics in the free atmosphere, far from surface frictional effects, is exemplified by recent observations and simulations of wave trapping occurring even at high altitudes, such as near the tropopause (e.g., [29]).

Note that the accuracy of the DEM is contingent upon the correct identification of the dominant Scorer parameter maxima; thus, it is not recommended for flows where the Scorer parameter vertical structure is poorly defined or highly transient.

6. Conclusions

This study introduces and validates the Drag Empirical Model (DEM), a parametric framework for estimating total pressure drag, as a function of nonlinearity ($D\left({\widehat{h}}_0\right)$), in moderate to strong non-hydrostatic flows, characterized by trapped lee waves (TLW). By leveraging the two-layer Scorer atmosphere as a physical baseline, the DEM provides a computationally inexpensive alternative to high-resolution numerical simulations, offering a physically grounded approach to drag estimation in realistic atmospheric conditions.

The primary contributions of this work are threefold:

• Normalization Strategy: Normalizing total drag by its linear counterpart ($D_L$, obtained via Drag Linear Model) provides a robust unifying framework. Since $D_L$ is efficiently calculated with negligible computational cost, this approach effectively collapses the complex dependence of drag on nonlinearity ($D\left({\widehat{h}}_0\right)$) into a set of non-linear correction functions suitable for real-time applications. These functions are defined by the mode number n, the duct depth H, and the Scorer parameter at the interface ($z=H$).
• Modal Dependence: Non-linear drag enhancement in the presence of TLW is critically dependent on the mode number n. While the fundamental mode ($n=1$) largely maintains the quadratic mountain-height scaling characteristic of linear theory, higher-order modes exhibit marked divergence as nonlinearity increases. The DEM successfully accounts for this behaviour through its mode-dependent formulation.
• Validation and Robustness: The model demonstrates strong agreement with FLEX numerical simulations across a diverse range of soundings. The high consistency of the DEM across the majority of the cases investigated suggests that, in scenarios where numerical simulations exhibit instability, the empirical model provides a more physically consistent representation of the drag trend, effectively bypassing numerical artifacts associated with near-singular resonance zones.

By capturing realistic trends such as resonance and saturation under non-linear conditions, the DEM provides a robust strategy for improving orographic drag parameterizations in weather and climate models. Its implementation can significantly reduce the computational overhead of resolving non-hydrostatic effects over steep terrain. Future research should incorporate surface friction and three-dimensional topographic effects to further extend the applicability and realism of this empirical framework.