The effect of a stable boundary layer on orographic gravity wave drag

Numerical simulations are carried out using the Weather Research and Fore-cast (WRF) model to calculate explicitly the ratio of orographic gravity-wave drag (GWD) in the presence of a stable boundary layer (BL) to the inviscid drag in its absence, either obtained from inviscid WRF simulations or estimated using an analytical linear model. This ratio is represented as a function of three scaling variables, defined as ratios of the BL depth to the orography width, height, and stability height scale of the atmosphere. All results suggest that the GWD affected by the stable BL, D BL , is inversely proportional to the BL depth h BL , roughly following D BL ∝ h − 2BL . The scaling relations are calibrated and tested using a multilinear regression applied to data from the WRF simulations, for idealised orography and inflow atmospheric profiles derived from reanalysis, representative of Antarctica in austral winter, where GWD is expected to be especially strong. These comparisons show that the scaling relations where the drag is normalised by the analytical inviscid estimate work best. This happens because stable BL effects reduce the amplitude of the waves above the BL, making their dynamics more linear. Knowledge of the BL depth and orography parameters is sufficient to obtain a reasonable correction to the invis-cid drag without needing additional information about the wind and stability profiles. Since the drag currently available from numerical weather prediction model parametrizations comes from linear theory uncorrected for BL effects, the results reported here may be applied straightforwardly to improve those parametrizations.

Medium-Range Weather Forecasts (ECMWF) GCM by using an envelope orography, which increases the mean orography height to account for subgrid-scale variability, but the success of this approach was limited, as it is insensitive to the horizontal scale of the orography.The first orographic gravity-wave drag (GWD) parametrization was formulated by Palmer et al. (1986), to reduce a systematic westerly bias in the midlatitude circulation in the ECMWF GCM.Current drag parametrizations (e.g., Lott and Miller, 1997) still do not include the effects of the boundary layer (BL), although interaction between GWD processes and the BL has been shown to have an important impact on the representation of GWD in models (Ólafsson and Bougeault, 1997).
A number of studies have investigated the ways in which gravity waves and the BL can interact.Using gravity-wave-resolving 3D numerical simulations, Ólafsson and Bougeault (1997) found that the BL reduces GWD magnitude and delays the onset of wave breaking.Peng and Thompson (2003) hypothesised, based on wave-resolving 2D numerical simulations, that GWD is reduced because the top of the BL acts as an attenuated effective orography, forcing weaker waves, which exert less drag than would exist in inviscid conditions.
In a theoretical article, Smith et al. (2006) investigated how gravity waves are absorbed by the BL, and developed a 2D linear model of that process, where the free atmosphere is coupled to the ground using a bulk BL representation, with Rayleigh damping coefficients acting on the differences between the wind speed within the BL and both the ground and the free atmosphere.They first applied this model to understand the absorption of trapped lee waves by the BL.Subsequently, Smith (2007) used a 3D version of the same model to study BL effects on vertically propagating mountain waves, finding that the flow divergence associated with upstream shift of the wave pattern within the BL shifts the wind maximum from the downwind slope towards the mountain top and modulates the depth of the BL so that incident gravity waves are partially absorbed, reducing the GWD.
Using high-resolution numerical simulations and the theory of Smith et al. (2006), Jiang et al. (2008) evaluated in more detail the impact of the BL on gravity waves excited by flow over 2D ridges and 3D axisymmetric mountains.They found that the upstream shift of the wave pattern and weakening of the wave aloft caused by the BL lead to a reduction in GWD (by up to 60%) and momentum flux (by up to 80%).These changes were found to be governed by the nondimensional BL height ( H) and the BL shape factor (γ), defined as the ratio of the displacement thickness to the momentum thickness (Schlichting et al., 2000).Phase shift and drag reduction were found to be greater over rough surfaces, which correspond to larger H and smaller γ.BL effects were additionally seen to decrease with increasing nonlinearity; nevertheless, the inclusion of the BL was found to shrink the linear flow regime, at least within the BL itself (Jiang et al., 2008).
Based on their numerical simulations, Jiang and Doyle (2008) found a substantial variation in GWD caused by changes in the characteristics of the BL over a diurnal cycle, which was consistent with their estimates from linear theory.In the daytime, a convective BL develops, with a shallow layer of shear near the surface and a deep well-mixed layer above.Both of these features tend to reduce GWD, therefore it was found that the daytime BL can weaken mountain waves significantly and reduce momentum flux by up to 90% (Jiang and Doyle, 2008).At night-time, a shallow stable BL develops, whose alteration of GWD is governed by the Froude number (F = U/(Nh), where U is the wind speed, N is the Brunt-Väisälä frequency, and h is the mountain height).If the BL flow is supercritical (F > 1), the drag increases as F → 1 to reach a maximum at F ∼ 1, where the drag may be several times its inviscid value.If the BL flow is subcritical (F < 1) due to excessive cooling, the drag decreases with decreasing F (Jiang and Doyle, 2008).
Although guided by the theoretical framework of Smith et al. (2006) and Smith (2007), the studies of Jiang et al. (2008) and Jiang and Doyle (2008) convey a complex yet incomplete picture of BL effects on drag.Possible empirical drag scaling laws presented in table 4 of Jiang et al. (2008) apply to very specific orography types and parameters (e.g., 2D ridges or 3D axisymmetric mountains with 10 km half-width), and the drag reduction is expressed in terms of parameters that are perhaps not the most convenient or physically significant, such as the Rayleigh damping coefficient C B used in Smith's bulk BL model.
The work presented here extends Jiang et al. (2008) by investigating the effects of a stable BL on GWD over orography with an elliptical horizontal cross-section, which is the idealised form assumed in a number of GCM orographic drag parametrizations (Lott and Miller, 1997;Elvidge et al., 2019), but takes a different scaling approach.We hypothesise that the height of the stable BL is the key parameter determining surface GWD reduction.Clearly, the GWD would not be reduced if the BL was absent (i.e., had zero thickness), and must necessarily approach zero as the BL thickness increases indefinitely (although other types of drag, e.g., frictional drag, remain).The BL height is a readily available quantity (from its corresponding parametrization) in the NWP models where the GWD needs to be parametrized.The aim here is to develop a scaling for the drag reduction as a function of a nondimensional BL depth (to be defined) that could correct the inviscid drag currently used in NWP model parametrizations.
To develop this scaling, high-resolution numerical simulations are run using vertical profiles of wind and potential temperature, either from reanalysis or adapted from reanalysis above four locations around Antarctica during austral winter.This complements the study of Turner et al. (2019), which showed that at those locations and times GWD is especially strong and effects of vertical wind shear on GWD are likely to be important.
Section 2 of this article describes the numerical model setup and methods used for processing the model output, as well as a theoretical model based on Smith et al. (2006) that is used to help interpret the results.In Section 3, the results of numerical simulations with idealised and more realistic wind profiles, assuming an elliptical orography, are presented.In Section 4, simple scaling laws are inferred from the results obtained in Section 3. Section 5 presents the main conclusions of this study.

METHODOLOGY
The Weather Research and Forecast (WRF) numerical atmospheric model is used to simulate idealised flows over elliptical orography and to calculate explicitly the GWD and its dependence on the parameters of the inflow, orography, and the BL that is parametrized in the model.The numerical results are interpreted in the framework suggested by a linear wave model including the effects of a bulk BL, building on Smith (2007).These models are described in turn below.

WRF model
Version 3.8.1 of the WRF model is used to perform the idealised numerical simulations presented in this article.WRF is a fully compressible, nonhydrostatic, nonlinear model (Watson and Lane, 2012), based on discretisation using an Arakawa-C grid and an explicit time-stepping scheme (Chen and Lin, 2004).For the simulations presented here, a horizontal grid of 100 × 100 points with 200 vertical levels is used.The horizontal grid spacing is 2 km and the vertical level spacing is ∼100 m at the surface and stretched aloft to values ∼170 m.This relatively high resolution allows a good representation of gravity waves, which is especially important in the case of structured atmospheric profiles (with directional critical levels; Guarino et al., 2016), and is sufficient to resolve the stable BL in the present simulations.For example, the equilibrium thickness of the BL in the simulations using realistic atmospheric profiles from reanalysis (in the range 430-630 m) only spans 5-7 grid levels, but is consistent with corresponding values extracted directly from reanalysis (in the range 320-500 m).The top of the domain is at a height of 20 km, and a dampening layer is applied over the top 5 km of the computational domain to prevent wave reflections.Open boundary conditions are applied to the lateral boundaries.An elliptical mountain with prescribed parameters (similar to those adopted by Phillips, 1984;Teixeira and Miranda, 2006), expressed by where h 0 is the mountain height and a and b are the lengths of the main axes of the ellipse, is placed in the centre of the domain.Two simulations are run for each of the atmospheric profiles specified below using the dry dynamical core of the WRF model, one with a BL scheme and one without.The BL scheme chosen is the Mellor-Yamada-Nakanishi-Niino (MYNN) third-level turbulent kinetic energy (TKE) scheme.The reason for this selection is described below in Section 3.1.This is a TKE scheme with third-level closure, requiring additional prognostic equations for the TKE and temperature variance, and determining the fluxes of momentum and heat diagnostically via flux-gradient relationships (Nakanishi and Niino, 2004).The eddy viscosity and diffusivity are commonly expressed as K c = S c l √ e, where l is the mixing length, e is the TKE, and S c is a proportionality coefficient (Shin and Hong, 2011) dependent on the TKE, the mixing length, and vertical gradients of the mean wind and temperature (Nakanishi and Niino, 2004).TKE schemes differ in how they define these variables.In the case of MYNN, the mixing length is determined by the surface layer length, turbulent length (i.e., the typical length-scale of the turbulent eddies), and the buoyancy length (Zhang et al., 2011).Simulations are run for 48 hr to allow the flow and BL to adjust to an approximate equilibrium and the drag to stabilise.All simulations neglect rotation of the Earth (i.e., the Coriolis parameter is set to f = 0).This approximation will be shown to be reasonable below.The integration time-step is 10 s, giving a Courant number typically far below the CFL stability criterion for the flows simulated.

WRF simulations
Each WRF simulation is initialised with an input atmospheric profile and in this instance these are taken from four different points around Antarctica in the ERA-Interim reanalysis dataset (cf.Turner et al., 2019).Antarctica is chosen because GWD is strong and because the cold temperatures occurring there mean that the BL is usually shallow and stable.These conditions are ideal for studying the effect of the stable BL in the situations that are most relevant for global atmospheric GWD (Gregory et al., 1998;Scinocca and McFarlane, 2000).The locations and dates used are as follows (all dates at 1200 UTC): • (77 The dates chosen are all in austral winter to ensure stronger gravity-wave activity and hence greater GWD.On the selected dates, the GWD enhancement by wind shear we calculated in Turner et al. (2019) was expected to be relatively strong in the corresponding location at that time.Although relevant, this effect is likely to be weaker than that of the BL, on which we are focusing here.Note that location D corresponds to Davis station, which is located on Princess Elizabeth Land and is maintained by the Australian Antarctic Division.Although the coordinates listed for location D are exact, the profile data used are from the nearest grid point in the ERA-Interim one-degree dataset.Figure 1 shows the four locations on a polar stereographic map.Table 1 details the mountain height, axis half-widths and orientation of the subgrid-scale orography, and the Coriolis parameter at each location.The height of the subgrid-scale orography is defined as twice the standard deviation  of the subgrid-scale orography elevation, that is, h 0 = 2.The two axis lengths (a and b) are calculated in terms of other subgrid-scale orography parameters.They are defined as a = ∕ and b = ∕(), where  is the slope and  the aspect ratio ( = a∕b; see Lott and Miller, 1997).The orientation is defined as the angle the major axis of the ellipse makes with north.An angle of 0 • (i.e., aligned with north) denotes meridional orography, whilst an angle of ±90 • denotes zonal orography.
At each location, simulations are run with an elliptical mountain corresponding to the orography parameters specified in Table 1, with the following setups (based on ERA-Interim data): • Constant wind (U, V) and Brunt-Väisälä frequency N profiles based on the values of wind velocity and potential temperature at the first five pressure levels above the surface (with a spacing of 25 hPa, corresponding roughly to 200 m).These levels are chosen in order to cover the full depth of the stable BL in all cases.
• Real wind and Brunt-Väisälä frequency profiles.These are linearly interpolated from the ERA-Interim grid to the WRF model grid.
As it is convenient that the axes of the orography are aligned with the WRF model grid, the winds specified in the input sounding are rotated so that they are given in the frame of reference aligned with the axes of the subgrid-scale orography.Two simulations are run for each input profile, one with a BL scheme and one without.The BL scheme to be used is selected below using preliminary idealised simulations.
The surface drag is calculated from the WRF model output as the integral of the pressure perturbation induced by the wave multiplied by the terrain slope over the whole computational domain.This is defined (Teixeira, 2014) as where p is the pressure perturbation and h(x, y) is the terrain elevation.To evaluate Equation 2 numerically, we use finite differencing to approximate the derivatives of the terrain elevation, and estimate the pressure perturbation at each time-step by subtracting the mean pressure over the horizontal domain at that time-step from the pressure at each grid point.The pressure perturbation and terrain slope are then multiplied together and summed over the whole domain.The magnitude of the drag vector in Equation 2 is given as Note that this definition accounts for all possible types of drag, including frictional drag (when a BL scheme is used), but for orography with the characteristics of that addressed here (tens of km wide) GWD is expected to dominate.

Effect of the boundary layer
To test the effect of the stable BL on the drag, initially the output of simulations run using the constant profiles described above, both with and without a BL scheme, is compared.The use of constant profiles enables the effect of the BL to be isolated from other atmospheric profile effects, which would complicate the derivation of scaling laws from the drag behaviour.However, the plausibility of the scalings will then be confirmed in Section 4.2 using the real atmospheric profiles.
To examine the impact of the stable BL on GWD magnitude, the ratio of the drags obtained at the end of simulations with and without the BL is calculated.The ratios of the depth of the BL to various relevant length-scales in the problem are also evaluated, as we hypothesise that these control the drag ratio.The BL depth is diagnosed by the WRF model and has a value at each grid point.The BL depth is diagnosed based on a hybrid method that takes into account both a potential temperature threshold exceedance relative to its minimum in the BL, and a TKE threshold (especially for stable BLs; Olson et al., 2019).This latter criterion is likely to be the relevant one in the situations addressed here.A default roughness length value of z 0 = 0.1 m is adopted, which reflects a compromise between the expected smooth icy surface and the terrain roughness associated with small-scale topographic features in a mountainous area.We expect that the drag will be reduced by the stable BL, as seen by Jiang et al. (2008) for simpler (i.e., axisymmetric) orography, and that this reduction will be more significant when the depth of the BL is large compared with either the height or width of the orography, or the stability length-scale in the atmosphere.
An important aspect of this investigation is the impact that the stable BL has on the amplitude of the gravity waves.To examine this, the amplitude of each horizontal wind component over the peak of the mountain (A U and A V ) is calculated by subtracting its minimum value from its maximum value (both above the top of the BL) and calculating the magnitude, that is, The reduction in this quantity caused by the introduction of the BL is evaluated by taking the ratio between the values A BL with and A inv without the BL.The relationship between this ratio and the drag will be examined.

Theoretical model for orographic gravity-wave drag
Comparisons are also made with drag values calculated using a simple model based on the linear theory of Smith et al. (2006).The model is both linear and nonhydrostatic, and does not account for the effects of the Earth's rotation (consistent with the WRF runs).The drag associated with vertically propagating and evanescent waves is given, respectively, by where  0 is a reference (constant) density, (k, l) is the horizontal wavenumber vector, ĥ(k, l) is the Fourier transform of the (3D) orography and H B is the BL depth.N is the Brunt-Väisälä frequency of the incoming flow and where (U, V) is the (constant) incoming wind velocity outside the BL, and C B and C T are the Rayleigh damping coefficients (inverse time-scales) coupling the BL flow to the ground and to the free atmosphere, respectively, in Smith (2007).The total drag is the sum of that due to vertically propagating and evanescent waves, and its magnitude is To obtain the inviscid drag, it is sufficient to set the BL height (H B ) to zero, which yields For an elliptical mountain, described by Equation 1, this can be evaluated analytically to yield the expressions originally derived by Phillips (1984).In the inviscid approximation, there is no drag associated with evanescent waves.
The magnitude of the inviscid drag is given by The integrals in Equations 3, 4 and 8 are evaluated numerically using a Gauss-Legendre quadrature algorithm.
If the drag given by Equations 3 and 4 is made nondimensional by dividing D theo by 0 , it becomes a function of six independent nondimensional flow parameters (an aspect that appears not to have been noted explicitly before): the angle between the incoming wind and the main axes of the orography , the orography aspect ratio , a measure of how nonhydrostatic the flow is N(ab) 1/2 /(U 2 + V 2 ) 1/2 , a nondimensional height of the BL H B /(ab) 1/2 , the ratio of the two Rayleigh damping coefficients in the bulk BL model of Smith (2007) C B /C T , and a nondimensional Rayleigh damping coefficient C T (ab) 1/2 /(U 2 + V 2 ) 1/2 .All these parameters have been made nondimensional using the wind speed (U 2 + V 2 ) 1/2 and a weighted mountain half-width (ab) 1/2 , so that they are as general as possible for all wind directions and all orography anisotropies.

Preliminary tests of WRF simulations
We firstly consider flow over an axisymmetric bell-shaped mountain, defined by Equation 1 for the case b = a (Smith, 1980).In our test case, h 0 = 10 m and a = 10 km.The input profile is specified by a uniform wind velocity U = 10 m/s and V = 0, and a potential temperature increasing linearly from 293 K at the surface up to the model top at 20 km, with a Brunt-Väisälä frequency N = 0.01 s −1 .The flow is approximately hydrostatic, since Na/U = 10, and also highly linear, since the flow nonlinearity parameter is Nh 0 /U = 0.01.The purpose of these test simulations is to determine whether the model gives an accurate value of the drag and to select an appropriate BL scheme to use in the subsequent simulations.A large number of BL schemes exist within the WRF model and it is important to select one that behaves realistically.Hu et al. (2010) compared the performance of three WRF BL schemes and found that the differences between them were largely due to differences in the strength of the vertical mixing and entrainment of air from outside the BL.Shin and Hong (2011) compared five BL schemes.They found that in unstable conditions nonlocal BL schemes in which the entrainment flux is proportional to the surface flux performed best, whilst in stable conditions local TKE closure schemes performed better (Shin and Hong, 2011).Since under Antarctic conditions the BL is expected to be stable, all the schemes considered below are TKE closure schemes.
The simulation is first run with no BL.The drag grows initially before stabilising in the first few hours of the simulation, to attain a final value of 89,179 N (Figure 2), which is approximated reasonably by the value of 93,399 N expected from linear theory (given by (∕4) 0 NUah 2 0 : see Teixeira, 2014).This corresponds to an error of less than 5% (probably attributable to the fact that the flow is not perfectly hydrostatic, which reduces the drag), indicating that WRF simulates the behaviour of low-amplitude waves adequately (Smith, 1980).
The simulation is then repeated with several BL schemes applied.The reduction of the GWD due to the top of the BL acting as an attenuated effective orography, as expected from Peng and Thompson (2003), occurs here, as was previously seen by Jiang et al. (2008).Using the MYNN third-level TKE scheme (MYNN), the drag produced at the end of the simulation is 11,585 N.For the MYNN 2.5-level TKE scheme (MYNN2: Nakanishi and Niino, 2004), the final drag is 27,204 N.For the Bougeault-Lacarrere scheme (BouLac: Bougeault and Lacarrere, 1989), the final drag is 9,737 N.This wide range of variation (by a factor of ∼3) is presumably due to differences in the way in which the schemes represent various BL parameters, for example different definitions of the turbulence length-scale.The drag stabilises within the time period of the simulations when the MYNN and MYNN2 schemes are used, but when the BouLac scheme is used the drag appears to be relatively stable at 30 hr before its value spikes up (see Figure 2).This may be due to some numerical instability, as similar behaviour is detected in the wind field for this scheme (not shown).The substantial noise in the drag signal near the end of the simulations with other schemes is due largely to the fact that we are dealing with small drag values (because the orography has low amplitude).The MYNN scheme is used in all subsequent experiments, as the drag value appears most stable in this case.Additionally, the higher order of the closure of this scheme makes it preferable, as higher-order closures are generally more accurate (Holt and Sethu, 1988), because they incorporate the transport of turbulent second-order moments, such as the TKE.The wide range of variation between schemes is not an obstacle to obtaining meaningful scalings from one of them, as the overall predicted effect of the BL on GWD is qualitatively similar.

Representative orography
Now the simulations that utilise representative orography, with an elliptical horizontal cross-section from ERA-Interim at each selected location, are discussed.This corresponds to the subgrid-scale orography that was adopted in the GWD parametrization scheme of the ECMWF model runs used to create the ERA-Interim dataset.Information about the subgrid-scale orography is given in Table 1, whilst information about the input profiles is provided in Table 2. Relevant nondimensional flow parameters are listed in Table 3.As can be seen, the inverse Rossby numbers in all situations are relatively small (< 0.3), for which linear theory predicts a drag not lower than between 80% (for 2D flow: Smith, 1979, his Figure 1) and 88% of its nonrotating value (for 3D flow: equation A15 of Miranda and James, 1992).This justifies neglecting the Coriolis parameter in the WRF simulations and analytical drag model, for simplicity.It is especially meaningful to use constant atmospheric profiles coming from the lowest levels in reanalysis in our tests, as orographic GWD is also parametrized using formulas that assume constant atmospheric parameters extracted from the lowest model levels (Lott and Miller, 1997;2020).Figure 3 shows the temporal evolution of the drag over the course of all simulations.These include simulations at location A (Figure 3a,b), location B (Figure 3c,d), location C (Figure 3e,f) and location D (Figure 3g,h), both with and without a BL.The thick black lines in Figure 3 denote simulations with real atmospheric profiles, while the simulations with constant atmospheric profiles taken at pressure levels 1, 2, 3, 4, and 5 nearest to the surface are denoted by the solid red, dashed blue, dotted magenta, dash-dotted green, and dash-double-dotted orange lines, respectively.We can see that, in general, the drag has stabilised by the end of the simulations in cases both without and with a BL.Exceptions to this are location C without a BL (Figure 3e), where the drag for the real atmospheric profile oscillates, and especially location A with a BL (Figure 3b), which seems to be subject to some form of numerical instability from hour 26 onwards.Detailed analysis of the flow field suggests that this is due to interaction of the BL scheme with the open lateral boundary condition in this case, producing an amplifying monochromatic plane wave that appears first at the northeast corner of the domain (the inflow side-see Figure 5c) and extends inwards.However, we checked that the orographic wave structure is well-established after t = 24 hr, and differs little from that of the other cases shown in Figure 3. Additionally, the instability only starts to affect the flow over the hill (which determines the GWD) at around t = 32 hr.For this reason, the results before the drag diverges will be considered valid and used later on.
The drag increases very fast (in about 2 hr) to near its final value in most simulations without a BL.However, in simulations with a BL, while increasing equally fast initially, it then decreases much more slowly over a period of about 24 hr.This is a consequence of not only the deepening of the stable BL, but also the flow development within it (as will be seen in more detail below), over a period much longer than it typically takes for the inviscid wave field to be established.The drag given by the simulations with constant atmospheric profiles using the lowest reanalysis levels brackets the drag given by the simulations with real atmospheric profiles, with the exception of the simulation for location D without a BL (Figure 3g) and the corresponding one with a BL during the first 12 hr (Figure 3h), and of the simulation for location A with a BL after hour 34 (Figure 3b), because of the above-mentioned numerical instability.The level from the constant-profile simulations that provides the best fit to the drag for the real-profile simulations at the final time varies between locations and between simulations with and without a BL.Nevertheless, a slight predominance of levels 3 and 4 can be detected, especially in cases with a

Note:
The numbers 1-5 at each location correspond to the lowest five pressure levels from reanalysis, where the variables are taken.z agl is the height above ground level.The wind components (U, V) are defined in a frame of reference aligned with the main axes of the subgrid-scale orography. 0 is the potential temperature at the surface and N is the Brunt-Väisälä frequency.
BL. From Table 1, this corresponds to heights above the ground somewhat larger than, but comparable to, the BL heights displayed in Figure 4 below.It is plausible that the anomalous behaviour of Figure 3g is caused by some nontrivial atmospheric profile effect, such as wave reflection with constructive interference at a sharp atmospheric interface, which is known to be able to amplify the drag substantially (Teixeira and Argaín, 2020).
Figure 4 shows the drag normalised by D 00 =  0 N(U 2 + V 2 ) 1∕2 (ab) 1∕2 h 2 0 for a selection of simulations at the four locations.The solid blue and red lines denote the drags given by WRF, without and with a BL, respectively.The dashed blue and red lines are the corresponding results given by the linear model.The green dash-dotted line describes the depth of the BL as defined by WRF, h BL .In the linear model including the stable BL, the input parameters have been adjusted to provide an agreement as good as possible with the drag from WRF near the end of the simulations.It was assumed that and H B /(ab) 1/2 = h BL /(ab) 1/2 .The first two assumptions imply that the corresponding dimensionless parameters are constants.This means, in particular, that C T scales as (U 2 + V 2 ) 1/2 /(ab) 1/2 .Alternative scalings, such as (U 2 + V 2 ) 1/2 /h 0 and (U 2 + V 2 ) 1/2 /h BL were tested, showing worse performance.The additional alternative scaling C T ∼ N was not tested, as it does not make sense physically (it would mean that the effects of BL friction would increase with static stability).It was assumed that H B = h BL , which equates the BL height in the linear model exactly with the BL height averaged over the whole computational domain h BL as diagnosed from the WRF model.This last assumption means that H B /(ab) 1/2 , unlike the preceding input parameters of the linear model, is defined totally by the input data from WRF.All inviscid input parameters (, , and N(ab) 1/2 /(U 2 + V 2 ) 1/2 ) are also determined from the WRF input conditions, so no assumptions have to be made about them.This defines all six input parameters of the linear model.

𝝍(
, and angle between the wind and the minor axis of the mountain, at each location, for each selected reanalysis level In the inviscid case (blue lines), the agreement between the linear model and WRF is fairly good for location A with an atmospheric profile taken from level 5 (Figure 4a) and even better for location B and levels 1 and 4 (Figure 4b,c), with some underestimation (of less than 13%) of the GWD from WRF by the linear model.This is consistent with the fact that the flow is more nonlinear (higher value of Nh 0 /(U 2 + V 2 ) 1/2 ) for location A than for location B (see Table 3).In Figure 4d, for location B and level 5, the linear model overestimates the drag given by WRF (by about 40%).It is plausible that this is due to the slightly higher nonlinearity of the flow and to the fact that the wind direction is only deviated 26 • away from the along-ridge direction.Wells et al. (2008) showed that, while the drag is increased in nonlinear flow that is across an anisotropic mountain, it is, on the contrary, decreased for flow along the mountain.Figure 4e corresponds to location C and level 5.This is the most nonlinear situation considered (see Table 3) and, since the flow is 41 • away from being along the ridge, a high-drag state is created, which is substantially underestimated by the linear model (by about 30%).
For location D and level 2 (Figure 4f), the drag from WRF is predicted especially well by the linear model, and for location D and level 5 (Figure 4g) a situation somewhat similar to that of Figure 4d, but even more marked, is displayed, with substantial overestimation of the WRF drag by the linear model (by about 80%).The reasons that can be inferred from Table 3 for this behaviour are that the flow is moderately nonlinear and the wind direction is only 15 • away from being parallel to the long axis of the mountain.
When a stable BL is taken into account, it becomes more complicated to interpret the behaviour of the drag, because it depends on the chosen BL model parameters.The initial decrease of the drag in the linear model is caused directly by the growth of the domain-averaged h BL (given by the green dash-dotted lines).It is clear from Figure 4 that this decrease is much faster than the one that occurs in WRF, the reason being that it is not only the thickness of the BL that matters for the gravity waves that produce the drag, but also the characteristics of the flow within it, which evolve more slowly.The drag from WRF stabilises by the end of the simulations, except perhaps in location A and level 5 (Figure 4a).The linear model overestimates the final WRF drag in Figure 4a,b,e (by about 31%, 71%, and 75%, respectively), underestimates it in Figure 4g (by about 27%), and gives roughly a correct estimate in Figure 4c,d,f.The height of the BL has a somewhat erratic behaviour in some cases, keeping growing slightly towards (b,d,f,h) drag D BL from simulations with a BL.Black thick lines: using real atmospheric profiles interpolated from reanalysis; red solid, blue dashed, magenta dotted, green dash-dotted, and orange dash-double-dotted lines: uniform atmospheric profiles based, respectively, on data from pressure levels 1, 2, 3, 4, and 5 above the ground, from reanalysis (a) the end of the simulations in Figure 4e or even starting decreasing in Figure 4d,f.Overall, and with the input parameters estimated previously, the linear model does a good job of predicting the decrease of the drag due to the existence of the stable BL (despite the moderate relative errors).This will be shown more systematically in Figure 6.The drag behaviour in two cases presented in Figure 4 can be interpreted in the light of the wind and BL height maps presented in Figure 5. Figure 5a,b shows wind vectors and wind speed (shading) at the lowest model level from WRF for the simulations without a BL, for location A and level 5, and location C and level 5, respectively.Figure 5c,d shows the same, but for the WRF simulations with a BL.Finally, Figure 5e,f shows maps of the BL height for the same cases.In Figure 5a it can be seen that, for location A and level 5 (weakly nonlinear conditions), the flow is accelerated downstream of the mountain, as is typical of gravity waves, but this happens to a moderate degree.This acceleration is much more pronounced (resembling a downslope windstorm) in Figure 5b for location C and level 5, where the flow is considerably more nonlinear (see Table 3).This is consistent with the GWD amplification observed in this case.When a stable BL exists in the simulations, there is a clear wake of decelerated flow downstream of the mountain at location A and for level 5 (Figure 5c), which corresponds (via mass conservation) to a thickening of the BL (Figure 5e).For location C and level 5, apart from the fact that this wake is preceded by a substantial flow acceleration on the downstream slope of the mountain, there are incipient recirculation zones either and BL depth h BL as a function of time, for a number of selected cases using constant atmospheric profiles (see Table 2).(a) Location A, data from level 5; (b) location B, data from level 1; (c) location B, data from level 4; (d) location B, data from level 5; (e) location C, data from level 5; (f) location D, data from level 2; (g) location D, data from level 5. Solid blue lines: inviscid drag from WRF; dashed blue lines: inviscid drag from the linear model; solid red lines: drag with a BL from WRF; dashed red lines: drag with a BL from linear model; green dash-dotted lines: BL depth from WRF side of the wake (Figure 5d).Both phenomena denounce the higher nonlinearity of the flow (cf.Miranda and James, 1992).The wake is equally associated with a thicker BL (see Figure 5f).
Figure 6 compares the final value of the drag from each WRF simulation with the theoretical drag calculated using the linear model.The blue symbols denote the results without a BL and the red symbols the results with a BL.It is clear that, in the majority of cases, the drag calculated using the WRF model output is close to the theoretical value, particularly in cases including a stable BL.It is interesting that the values calculated using simulations for location C without a BL are mostly greater than the theoretical values (shown by the fact that they lie above the one-to-one line).This indicates some nonlinear drag enhancement, as the orography (and hence the nonlinear parameter) is larger at this location.Conversely, some drag values from location D are overestimated by the linear model (as discussed in the description of Figure 4).The drag values for the results including a BL are in comparatively good (perhaps even better) agreement between WRF and the linear model (apart from one data point at location C).This, and the intrinsically lower gravity-wave amplitude that accompanies these lower values of the drag (see below), suggests that the flow outside the BL (where the gravity waves propagate) is closer to linear when the BL is present, confirming the ideas of Peng and Thompson (2003).

SCALING THE DRAG REDUCTION AND ITS DEPENDENCE ON BOUNDARY-LAYER DEPTH
The purpose of this section is to try to assess in as systematic a way as possible how the ratio of the drag affected by a stable BL, D BL , to the inviscid drag, whether given by linear theory (D 0 ) or WRF simulations (D inv ), varies with measures of the impact of the BL.Even if we limit ourselves to the dependence of the GWD on quantities involving the domain-averaged BL depth h BL (given its obvious relevance), this quantity can be normalised using key flow parameters in three possible ways: h BL /(ab) 1/2 , h BL /h 0 , or h BL N/(U 2 + V 2 ) 1/2 .Next, we will develop scalings for D BL /D 0 and D BL /D inv in terms of all these dimensionless parameters.The theoretical model of Section 2.4 suggests that D BL /D 0 may also depend on flow parameters such as the wind angle , the aspect ratio of the mountain , and the nonhydrostatic N(ab) 1/2 /(U 2 + V 2 ) 1/2 parameter.In nonlinear conditions, there is a possible additional dependence on Nh 0 /(U 2 + V 2 ) 1/2 .However, since these are inviscid parameters they should affect both D BL and D 0 , or D BL and D inv , in a roughly similar way, so when the ratios D BL /D 0 and D BL /D inv are calculated the influence of those parameters should factor out, at least partially.
We have seen in Section 3.1 that the existence of a stable BL decreases the GWD very substantially, so it is plausible to assume that D BL is inversely proportional to h BL in some way.Therefore, we propose the following scaling relation: where , , , and  are adjustable coefficients, and D ref may be either D 0 or D inv .If the logarithm of Equation 9is taken, it becomes clear that these coefficients may be fitted to a given dataset where all the included variables are available using multilinear regression.This will be done next.

Atmospheric profiles with constant parameters
Because the gravity waves have smaller amplitude and therefore are more linear in the presence of the stable BL (as shown by the good agreement between the GWD predicted by the linear model and by WRF in Figure 6), it makes sense to try to make this fit first for D BL /D 0 .Figure 7 shows D BL /D 0 as a function of h BL /(ab) 1/2 (Figure 7a), as a function of h BL /h 0 (Figure 7b), as a function of h BL N/(U 2 + V 2 ) 1/2 (Figure 7c), and as a function of a combination of these three parameters obtained by the multilinear regression procedure introduced above (Figure 7d).The regression was only adjusted to the data from locations A, B, and D, since location C, because of its higher nonlinearity, was seen to be an outlier, and to affect an optimal fit of the data from the other locations detrimentally.The adjusted coefficients of the multilinear regression were found to produce the following form for the scaling relation: This corresponds to D BL ∝ h −2.069 BL , which makes sense physically, as the GWD varies in inverse proportion to h BL .In the dataset from the WRF simulations, all parameters h BL /(ab) 1/2 , h BL /h 0 , and h BL N/(U 2 + V 2 ) 1/2 vary simultaneously, which makes interpreting the data complicated and sometimes misleading.In Figure 7a,b, the scatter in the data is considerable.Although the expected decreasing trend of D BL with h BL for the data from each location is found, the behaviour between locations differs between Figure 7a and 7b, with the trends that can be seen at each location not aligning between locations and thus failing to collapse into a single trend.For example, in Figure 7a, D BL /D 0 increases as h BL /(ab) 1/2 increases between locations instead of decreasing, which is a sign that the appropriate scaling for D BL /D 0 is a combination between the two variables on the horizontal axes of Figure 7a,b.In Figure 7c, the scatter is substantially smaller, but the alignment of the data cloud is more vertical.This is a consequence of the fact that, in an atmosphere without rotation or where the effect of stability is dominant, the controlling factor for BL growth is N. Therefore, h BL would be expected to scale as (U 2 + V 2 ) 1/2 /N (corresponding to a vertical line).In practice, there is some departure from this relation, shown by the fact that, for the data from each location, there still seems to be some decreasing trend of D BL /D 0 as h BL N/(U 2 + V 2 ) 1/2 increases.Figure 7d shows D BL /D 0 as a function of the scaling variable on the right-hand side of Equation 10, with the dashed line denoting the corresponding fitting equation.We can see that the fit to data from locations A, B, and D (which were used to adjust the parameters) is very good.In particular, all three datasets now align along a single power-law trend.Data from location C (which were excluded from the linear regression procedure) are naturally not as well fitted, except for a single data point.It is noteworthy that 16 out of a total number of 20 data points are very well fitted by Equation 10.This gives us some confidence in this relationship.In Figure 7e, the multilinear regression is computed again, but assuming  = 0 in Equation 9.This corresponds to ignoring the scaling variable h BL N/(U 2 + V 2 ) 1/2 in the fitting procedure, given its weak correlation with the drag, found in Figure 7c.It yields which corresponds to D BL ∝ h −1.807

BL
. The advantage of using only h BL /(ab) 1/2 and h BL /h 0 as scaling variables is that only orography characteristics are involved, simplifying the application of Equation 11.However, in order to estimate D 0 , N and (U, V) are still necessary.As can be seen in Figure 7e, the fit provided by Equation 11 is not noticeably inferior to that provided by Equation 10, which means that both relations are probably useful.
Figure 8 shows the same as Figure 7, but for D BL /D inv , that is, normalising D BL using the inviscid drag from WRF instead of its linear estimate.One point worth noting is that, since D BL /D inv and D BL /D 0 differ by the factor D inv /D 0 , which is a purely inviscid quantity, this ratio should not depend on h BL .Therefore, in the form of the multilinear fit, Equation 9, that applies to D BL /D inv , the constraint  +  +  = 0.773 + 0.863 + 0.433 = 2.069 must be fulfilled as a consequence of the preceding results for D BL /D 0 expressed by Equation 10.With the enforcement F I G U R E 7 Variation of the drag given by WRF including a stable BL, D BL , normalised by the inviscid drag given by the linear model D 0 with (a) the ratio of the BL height to the weighted mountain half-width, h BL /(ab) 1/2 , (b) the ratio of the BL height to the mountain height h BL /h 0 , (c) h BL N/(U 2 + V 2 ) 1/2 , (d) an optimal combination of h BL /(ab) 1/2 , h BL /h 0 , and h BL N/(U 2 + V 2 ) 1/2 , and (e) an optimal combination of h BL /(ab) 1/2 and h BL /h 0 , for all simulations with constant atmospheric profiles.Black circles: location A; red squares: location B; blue triangles: location C; green diamonds: location D; dashed line: linear fit to data from locations A, B, and D Figure 8a-c shows D BL /D inv as a function of h BL /(ab) 1/2 , h BL /h 0 , and h BL N/(U 2 + V 2 ) 1/2 , respectively.The behaviour is relatively similar to that in Figure 7, except perhaps that trends for the data from each location are more difficult to discern. Figure 8d shows the variation of D BL /D inv with the scaling variable on the right-hand side of Equation 12. Again, the dashed line corresponds to the corresponding fitting equation.It can be seen that the scatter is substantially larger than for the D BL /D 0 fit, coming not only from the data for location C (which was not fitted here, for consistency) but also from data for location D, which does not follow the predicted trend at all.Reasons for this could have to do with the fact that D inv is unusually low for location D compared with the linear estimate (as exemplified in Figure 4g), due to nonlinear processes.This situation is not improved if  = 0 is imposed (i.e., if the effect of the stability height scale is ignored), yielding the scaling In this fit, and for similar reasons to those previously, the constraint  +  = 0.937 + 0.870 = 1.807 has been enforced, to ensure that D inv /D 0 does not depend on h BL (by reference to Equation 11). Figure 8e shows that the scatter increases substantially, the problem with location D (which was included in the fit) is not solved, and one point from location C (not fitted) becomes a distant outlier.
A noteworthy effect is that the fitted exponents in the power laws are different between Equations 10 and 12, and between Equations 11 and 13.This in principle reflects (given the imposed consistency constraints in each pair of equations) how the nonlinearity and nonhydrostatic effects differ between the linear model and the WRF simulations.Indeed, for example from Equations 10  and 12, (14) This scaling, which is not supposed to be accurate (actually D inv → D 0 when Nh 0 /(U 2 + V 2 ) 1/2 → 0, yet Equation 14predicts instead that D inv → 0), nevertheless gives a rough description of this dependence on the flow nonlinearity (quantified by Nh 0 /(U 2 + V 2 ) 1/2 ) and on nonhydrostatic effects (quantified by N(ab) 1/2 /(U 2 + V 2 ) 1/2 ).That description seems to be at least qualitatively correct, since it predicts the GWD to increase with increasing nonlinearity and hydrostaticity of the flow (Teixeira, 2014).What the inaccuracy of this scaling for D inv /D 0 shows is simply that the forms of the scalings given by Equations 10 and 12 are too simple to evaluate these effects accurately (although they may be reasonable for their own purposes).
Analysis of the wind profiles over the mountain in the WRF simulations with and without a stable BL (examples shown in Figure 9) reveals marked differences between the two.The signature of the wave oscillation is obvious in all cases without a BL, although it is modest for location B, at which the orography is small (this is expected, as the amplitude of the orography determines the amplitude of the wave).The wave oscillation is substantially reduced in cases with a BL.This large reduction in wave activity is consistent with the detected reduction in drag magnitude, and is quantified here by the ratio of the wave amplitudes in the presence and absence of a BL, A BL /A inv , defined previously.
Figure 10 explores the relationship between the reduction in the wave amplitude A BL /A inv and the drag D BL /D inv .It indicates a strong relationship between D BL /D inv and A BL /A inv .In other words, the reduction in wave amplitude appears to be a leading factor determining how strongly the drag is reduced by the stable BL.All of this is not surprising, as the wave amplitude has a very strong link with the magnitude of the pressure perturbation that generates the drag.More specifically, D BL /D inv is proportional to A BL /A inv because the GWD is proportional to the pressure perturbation induced by the waves, and this pressure perturbation is proportional to the corresponding horizontal velocity perturbation, quantified here by A. hold for the WRF simulations using the real atmospheric profiles from ERA-Interim as input.Obviously these profiles are more complex than those considered previously, and different effects will also be important (for example, changes in stability with height, and possibly critical levels).Of particular note is the fact that these profiles include wind shear, which has previously been seen to have important impacts on the drag (e.g., Miranda et al., 2009;Turner et al., 2019;Xu et al., 2019).While, in this case, flow parameters such as a, h 0 , h BL , D BL , and D inv are unique for each simulation, N and (U 2 + V 2 ) are functions of height, so it is not a priori obvious what values to select when these parameters are used to normalise h BL or to calculate D 0 .

Real atmospheric profiles
For maximum generality, we will next present results with these variables taken from the five lowest pressure levels adopted in the constant-profile simulations.Another aspect to keep in mind is that in this case it is impossible to take D BL , D inv , and h BL from the WRF simulations at the final time, because of the nonphysical behaviour displayed in Figure 3b.For that reason, t = 32 hr is chosen instead, which is a time simultaneously before the period when that nonphysical behaviour begins but after D BL and h BL have approximately stabilised (see Figure 3). in Equations 12 and 13, respectively.The symbols in each panel correspond to data from the various locations, but here multiple points from a single location have a different meaning.Namely, in Figure 11a,c the set of five crosses for each location along the horizontal direction correspond to the scaling variable h BL N/(U 2 + V 2 ) 1/2 calculated using values of N and (U, V) from the lowest five pressure levels in reanalysis.Conversely, in Figure 11a,b, the set of five crosses for each location along the vertical direction corresponds to the inviscid drag from linear theory D 0 calculated using values of N and (U, V) from the lowest five pressure levels.Solid symbols denote data points taken at the levels which, at each location, are closest to the domain-averaged BL depth.The dashed lines in Figure 11a-d correspond to Equations 10-13, respectively.
In Figure 11a,b, it can be seen that the data roughly follow the regression line, with both a decrease of D BL /D 0 and an increase of the scaled BL depth as one goes from location A to D to B. If N and (U, V) are taken from the model just above the BL top, a good fit with the scaling for the reduction in GWD associated with the stable BL is obtained.We notice that there is not a substantial difference in the quality of the fit between Figure 11a and b, although, for a practical application of the scaling relation, the issue remains of where to estimate N and (U, V) to obtain D 0 (and normalise h BL in Figure 11a).Although the circles are not the data points closest to the dashed line, they are not the most distant either, departing from the line by a factor of about 1.5 or lower.This suggests that h BL might be a good reference level to estimate the necessary parameters.
In Figure 11c,d, where the GWD has been normalised instead using the inviscid drag from WRF, the behaviour of the data from the various locations is substantially less well predicted by the corresponding scaling.The data from the four locations bracket the theoretical prediction, with locations A and D especially close to the dashed line and most of the points from locations B and C being overestimated and underestimated by the line by factors of about 2, respectively, although in Figure 11c the full set of data for location C actually brackets the line.Note that, while the data point closest to the dashed line from location C corresponds to the lowest pressure level, the data point closest to the dashed line from location B corresponds instead to the highest level.Again, the solid symbols in Figure 11c, corresponding to the levels closest to h BL , do not show the best agreement with the dashed line, but definitely not the worst either.Although these results show that the order of magnitude of the drag attenuation is certainly predicted correctly, the datasets in Figure 11c,d do not follow the trend suggested by the line.Therefore, it is apparent that Equations 10 and 11 are more reliable than Equations 12 or 13, at least for the limited datasets that were tested here.Additional testing of the scaling relations using independent datasets would be useful, but is beyond the scope of this study.
The facts that all the scaling relations found above satisfy approximately D BL ∝ h −2 BL , and that the exponents  and  in Equation 11 are close to 1, that is, D BL ∕D 0 ∝ {h BL ∕(ab) 1∕2 } −1 and D BL /D 0 ∝ (h BL /h 0 ) −1 , are curious and should be no coincidence.The linear model (Equation 3) suggests a dependence of the type D BL ∝ h −2 BL when h BL is small and D BL ∝ h −1 BL when h BL is large, but a dependence on h BL /h 0 naturally cannot be justified using the linear theory framework, where the orography height is assumed to be infinitesimal.

CONCLUSIONS
It is evident from the results presented here that inclusion of a simulated stable BL reduces GWD substantially over orography with an elliptical horizontal cross-section, confirming and extending the findings of Jiang et al. (2008).
The effect is greater when the BL depth is large compared with other relevant length-scales in the problem, namely the orography height and width and the stability height scale of the atmosphere.The existence of this relationship corroborates and quantifies the findings of Peng and Thompson (2003) in the stably stratified flow regime, and suggests that there is a need to apply a correction to drag parametrizations in order to include this effect.In highly idealised conditions inspired by atmospheric profiles for locations in Antarctica from reanalysis, a power-law relationship was found to exist between the ratio of the drag with a stable BL to that without one and the corresponding scaled BL depth.The drag varies in inverse proportion to the BL depth, as would intuitively be expected, approximately as D BL ∝ h −2 BL .It scales inversely (via relations found using a multilinear regression) with three measures of the BL depth, normalised using the mountain width, height, and the static stability and wind speed evaluated at a representative level in the atmosphere.This dependence implies, in particular, in agreement with the findings of Jiang et al. (2008), that BL effects decrease with increasing nonlinearity, but also as the flow becomes more hydrostatic.This behaviour seems logical, since a greater depth of the BL relative to the orography height and width necessarily causes greater attenuation of the effective orography 'felt' by the free atmosphere aloft, forcing weaker waves.
The drag reduction due to stable BL effects potentially depends on other flow parameters, and its power-law dependence on the scaled BL depth, adopted here, is necessarily only approximate.Nevertheless, this dependence appears intuitively (and also by reference to the BL model of Smith et al., 2006) to be the dominant one.The relation between the drag reduction and the scaled BL height, established herein, could in principle be used to represent this effect in parametrizations under stable BL conditions.Even if current inviscid drag parametrizations could be calibrated to account implicitly for this effect by a judicious choice of adjustable coefficients, it would not be possible to obtain a physically based scaling of the drag with BL depth by such a procedure, and it is certainly preferable to represent separate physical processes more accurately.
The effect of the stable BL is naturally felt not only in the reduction of the drag magnitude, but also in various other aspects of the flow.The components of horizontal wind velocity are decelerated throughout the BL, and the wave activity was found to be reduced significantly above the BL, as was hypothesised above.It was seen that this reduction in wave amplitude is strongly related to the ratio between the drags with and without the BL, more specifically D BL /D inv ∝ A BL /A inv .This is expected, given the linear relation between the wave amplitude and the amplitude of the pressure perturbation that causes the drag.Although we have focused on quantification of the surface GWD, we expect the same amplitude reduction in the vertical profile of parametrized GWD, since the structure of the wave response is proportional, but reduced in amplitude.This has implications for the parametrization of the vertical distribution of the GWD acting on the atmospheric circulation (Xu et al., 2019).
Comparison between the results derived from the WRF model output and calculations using a linear model based on Smith et al. (2006) revealed that the two sets of results can behave quite similarly, as long as the friction coefficients in the linear model are properly adjusted.The linear model was found to perform even better in predicting the drag from WRF in the case where a stable BL is present than in the inviscid simulations.We speculate that this is due to the fact that the amplitude of the waves is reduced above the BL by BL effects, and therefore they are described better by linear theory.Inviscid waves, in contrast, have higher amplitudes, where nonlinearity and the accompanying flow regimes-wave breaking, lee vortices-make the drag depart more from its linear estimate.While it is certainly true that the flow within the BL is more nonlinear than that in the absence of a BL (because the flow perturbation becomes a larger fraction of the mean flow, which is reduced, as noted by Jiang et al., 2008), the flow outside the BL is definitely more linear.
The consequence of this fact is that normalising the GWD affected by a stable BL, D BL , by its inviscid estimate from linear theory, D 0 , yields a more universal scaling of its dependence on BL depth than normalising it by the value from inviscid WRF numerical simulations, D inv .Firstly, the scaling for D BL /D 0 fits the data from WRF better, either from the constant-profile simulations that were used to calibrate the scaling relation or from the simulations with a realistic atmospheric profile from reanalysis.Secondly, D 0 is more readily available than D inv , since the former is a quantity that inviscid drag parametrizations already use.When using Equation 10 to estimate D BL , there is still the issue of the height at which to estimate N and (U, V) included in the definition of h BL N/(U 2 + V 2 ) 1/2 .A more comprehensive justification for this choice is left for future studies, but our results suggest that a height within the range of the lowest five pressure levels from ERA-Interim reanalysis (spanning approximately the lowest 1,000 m of the atmosphere) is appropriate.More specifically, a height that is of the order of the BL depth seems to work acceptably.An alternative, but almost equally accurate, choice is to use Equation 11 to account for the effect of the BL instead.This has the advantage of depending only on the BL depth and orography parameters in the simulations, forfeiting the need to specify any atmospheric parameters.
The scaling relations developed here neglect potentially relevant effects and make many simplifying assumptions: for example, direct effects of flow nonlinearity (through Nh 0 /(U 2 + V 2 ) 1/2 ), nonhydrostatic effects (through N(ab) 1/2 /(U 2 + V 2 ) 1/2 ), effects of orography anisotropy (through ), and effects associated with wind direction (through ).In relation to this last aspect, note also that we only scaled the magnitude of the drag, but did not focus on how its two components are affected.It is quite plausible that the drag direction could change as a result of BL effects, not only associated with the Ekman spiral but also due to more subtle influences.Finally, the present scaling was developed and tested for strongly stable BLs.This approach is likely to have limited accuracy for other types of BL, for example, convective ones, as mechanisms generating GWD may, in those cases, be quite different (Jiang and Doyle, 2008).Treatment of nonstable BLs is left for future studies.

F
Drag as a function of time from WRF simulations for (a,b) location A, (c,d) location B, (e,f) location C, and (g,h) location D. (a,c,e,g) Drag D inv from inviscid simulations without a BL; Drag from WRF simulations D num and from the linear model D theo normalised by

F
Horizontal cross-sections of wind vectors and wind speed (contours) at the lowest model level above the surface and BL depth at the end of the WRF simulations for constant atmospheric profiles taken from pressure level 5 at locations A (left column) and C (right column).(a,b) Results for inviscid simulations; (c,d) results for the simulations with a BL.The arrow at the top of each panel indicates the velocity magnitude.(Note that the scales for (a,b) and (c,d) are different.Colour scales have units of m⋅s −1 .)(e,f) The depth of the BL (the colour scales have units of m).Arrows indicate the direction of the wind vector in the input profile Comparison of the drag from WRF D num and the drag from linear theory D theo normalised by D 00 for all cases with constant atmospheric profiles (all locations, and atmospheric parameters taken from all reanalysis levels).Blue symbols: inviscid WRF simulations and linear model calculations; red symbols: WRF simulations and linear model calculations including a BL.Circles: location A; squares: location B; triangles: location C; diamonds: location D; dashed line: 1 to 1 relation Figure 11 illustrates how various scalings for the drag work in these real-profile simulations.Figure 11a,b shows D BL /D 0 as a function of the scaling variables included in Ratio between the drag with and without a BL, D BL /D inv , as a function of the ratio between the wave amplitude in the WRF simulations with a BL and the wave amplitude in the corresponding inviscid simulations, A BL /A inv .Black circles: location A; red squares: location B; blue triangles: location C; green diamonds: location D Equations 10 and 11, respectively, and Figure 11c,d shows D BL /D inv as a function of the scaling variables included

F
I G U R E 11 Drag from the WRF simulations with real atmospheric profiles as a function of scaled BL depth variables similar to those used in (a) Figure7d, (b) Figure7e, (c) Figure8d, and (d) Figure8e.Black circles: location A; red squares: location B; blue triangles: location C; green diamonds: location D. In (a,c), the set of crosses along the horizontal axis corresponds to h BL N/(U 2 + V 2 ) 1/2 using values of N and (U, V) from the five lowest pressure levels, and in (a,b) the set of crosses along the vertical axis corresponds to D/D 0 with D 0 calculated using values of N and (U, V) from the five lowest pressure levels.Solid symbols denote results for the level closest to h BL in the WRF simulations (2 for locations A and C and 3 for locations B and D)

Location h 0 (m) a (km) b (km) 𝜸 = a∕b 𝜽( • ) f (s −1 )
Input profile data for each location TA B L E 2 Similar to Figure7, but with the drag D BL normalised by its inviscid value given by WRF, D inv , instead of D 0 , and with the mixed variables used in (d) and (e) defined differently (see text for details) Having established scaling relations for D BL /D 0 and D BL /D inv , we would like to see whether these relations