The resistance law for stably stratified atmospheric planetary boundary layers

The resistance law for stably and neutrally stratified atmospheric planetary boundary layers (PBL) entered textbooks on boundary‐layer meteorology but, until now, remains practically unused in modelling applications. This is not surprising. The law has been formulated and validated only for idealised cases, such as truly neutral PBL – implying neutral stratification across the entire atmosphere, nocturnal stable PBL – stably stratified near the surface but developed against the neutrally stratified free flow, and (more recently) conventionally neutral PBLs – neutrally stratified near the surface but developed against stable stratification in the free flow. We derive and validate the general formulation of the resistance law accounting for the integral effect on PBL of stable stratifications at the surface and in the free atmosphere. Such long‐lived stable PBLs, typical of wintertime at high latitudes, were until recently overlooked in boundary‐layer meteorology, not to mention weather and climate models. The proposed general formulation of the resistance law covers long‐lived stable PBLs and opens up prospects for their improved modelling.


INTRODUCTION
The classical resistance law for the steady-state atmospheric planetary boundary layer (PBL) links the geostrophic drag coefficient C g (the ratio between the friction velocity u * and the geostrophic wind |U g |, C g = u * /|U g |) and the cross-isobaric angle with the PBL governing parameters: Here, k ≈ 0.4 is the von Kármán constant, Ro = |U g |/(z 0 f ) is the surface Rossby number, z 0 is the surface roughness length, f is the Coriolis parameter, and A and B are dimensionless coefficients, originally supposed to be universal constants but later proven to be functions of stability parameters. The sign on the right-hand side of the second equation is plus in northern and minus in southern hemispheres. For brevity, we limit below to the Northern Hemisphere.
For the first time, this law has been derived by Rossby and Montgomery (1935) from the analytical solution to the steady-state momentum equations for horizontal wind velocity by setting wind velocity to zero at the surface and to geostrophic wind at the infinite height, accounting for the Earth's rotation and using the following vertical profile of eddy viscosity: Here, K * is eddy viscosity in the PBL core, z is the height over the surface, z * is the height of the near-surface layer with a logarithmic velocity profile that was taken proportional to the Ekman rotational height-scale, √ K * ∕f , dictated by the momentum equations (Ekman, 1905). Notably, such derivation yields the resistance law, Equation (1), for any vertical profile of K M provided that close to the surface it behaves as K M = ku * z. In this formulation, the shape of the profile at z > z * affects only the constants A and B.
A quarter-century later, Kazanski and Monin (1961) have derived the resistance law without solving momentum equations -by matching the near-surface logarithmic law, U = (u * /k)ln(z/z 0 ), with the velocity defect law, [U g − U(z)]/u * = f U (fz/u * ), where f U is the universal function of dimensionless height fz/u * . In this derivation the concrete shape of the velocity defect function was not needed, whereas coefficients A and B appeared as unknown dimensionless constants to be determined empirically.
It follows that the principal formulation of the resistance law, Equation (1), with the accuracy of undefined dimensionless coefficients A and B, can be derived from the solution to momentum equations using an obviously simplified approximation of eddy viscosity satisfying only two requirements: • close to the surface, it behaves as K M = ku * z, yielding the logarithmic velocity profile; • its typical value in the PBL is defined as with no concern for the shape of its vertical profile beyond the near-surface layer. The above derivations do not take into account stratification of density, which is why coefficients A and B in Equation (1) appear as dimensionless universal constants. This is proper only for Truly Neutral (TN) PBL characterised by zero buoyancy flux at the surface, F bs = 0, and neutral stratification in the free atmosphere aloft. The formulation accounting for non-zero buoyancy flux at the surface has been derived by Zilitinkevich (1967), Zilitinkevich et al. (1967) and Zilitinkevich and Chalikov (1968). They revealed theoretically and verified against observational data that Equation 1 holds true for Nocturnal Stable (NS) PBL with negative buoyancy flux at the surface, F bs < 0, and neutral stratification in the airflow aloft, whereas coefficients A and B became functions of the dimensionless parameter: It should be stressed that here and below the derivation of Equation (1) for stably stratified PBLs involves the scaling analysis and additional assumptions, thus the result may not seem as rigorous as the near-surface logarithmic law or the velocity defect law used to derive the resistance law for TN PBL. Therefore, the "law" term for stably stratified PBLs might be misleading. Nevertheless, the authors kept using the name of "resistance law", so we follow them in this regard for the sake of coherence. Subsequently, Zilitinkevich and Deardorff (1974) have: extended the resistance law analogous to Equation (1) to non-steady PBL regimes -with PBL height, h = h(t), dependent on time, t; employed the ratio h/z 0 instead of the surface Rossby number Ro; and defined the coefficients A and B as functions of h/L, where L is the Obukhov (1946) length-scale: Much later, Zilitinkevich and Esau (2002; have revealed theoretically and verified against both observational data and large-eddy simulation (LES) that PBLs with zero buoyancy flux at the surface, traditionally identified as merely neutral, are in fact strongly affected by the static stability in the free atmosphere, whereas Equation (1) holds true, but A and B become functions of another dimensionless parameter: where N is the Brunt-Väisälä frequency in the free atmosphere. By this means, two principally different types of the presumably neutral PBL were distinguished: TN PBL (F bs = 0 and N = 0), uncommon in the atmosphere but realisable in laboratory experiments, LES and direct numerical simulation (DNS); and conventionally neutral (CN) PBL (F bs = 0, N > 0), in fact, affected by strongly stable stratification in the free atmosphere aloft. These results were supported by further analyses (e.g. Hess and Garratt, 2002a;2002b;Hess, 2004) and are now widely recognised.
We recall that traditionally the effect of static stability inherent to the free atmosphere was fully overlooked, which is why the attributes "stable" and "nocturnal" with respect to PBL were used as synonyms. The above-quoted papers account for this effect and distinguish between three principally different types of stable PBL: • NS PBL with negative buoyancy flux at the surface (F bs < 0) and neutral stratification (N = 0) in the residual layer aloft, typical only over continents at mid-and low latitudes; • CN PBL with zero buoyancy flux at the surface (F bs = 0) affected by the free-flow stability (N > 0) aloft, typical over the open ocean; • Long-lived Stable (LS) PBL (F bs < 0, N > 0) typical of wintertime at high latitudes. Zilitinkevich and Esau (2005) have demonstrated the relevance of the resistance law, Equation (1), to any stable PBLs, and defined dimensionless coefficients A and B as universal functions of either for NS PBL or N for CN PBL. In the general case of LS PBL, where A and B are functions of two arguments, and N , the resistance law remained undefined, which strongly restricted its practical applicability. This article compensates for this drawback. We develop a theoretical model specifying A = F A ( , N ) and B = F B ( , N ), validate the model against LES and use it to build an algorithm linking the surface stress (the friction velocity u * and the cross-isobaric angle ) with the geostrophic wind (U g , V g ). This could be of use, in particular, for retrieving the atmospheric pressure field from satellite observations of sea waves. Such observations yield data about the surface stress, and the resistance law recalculates them into the geostrophic wind and, hence, horizontal pressure gradients (Brown and Levy, 1986;Velden et al., 2006;Monzikova et al., 2016). Until recently this method was based on the resistance law for neutrally stratified PBLs which strongly limited its accuracy.

THEORETICAL MODEL
We derive the resistance law for steady-state stably stratified PBL employing the Ekman equations where U and V are components of the wind velocity U = (U, V) and U g and V g are components of the geostrophic wind U g = (U g , V g ).
To solve these equations we propose the following three-layer approximation of eddy viscosity: Here, h * is the PBL height-scale defined by the Ekman (1905) formula: and K * is the eddy viscosity scale for PBL core defined by conventional formulation: where u T and l T are turbulent velocity and length scales. Following , we take u T equal to friction velocity for any neutral or stable PBL: and take l T proportional to the turbulent length-scale inherent to the concrete type of PBL: (Rossby and Montgomery, 1935) for TN PBL, (Zilitinkevich and Esau, 2002) for CN PBL, (Obukhov, 1946) for NS PBL.
Equations 7, 9-11 define the eddy viscosity scale K * and the log-layer height-scale z * = l T /k for the three types of PBL with the accuracy of dimensionless constants to be determined empirically. Following Zilitinkevich and Esau (2005) and , we determine h * ∼ √ K * ∕f and z * = l T /k over the whole range of neutral and stable PBLs (TN, NS, CN and LS) through the following interpolation giving priority to stronger turbulence-restriction mechanism: ( ( In view of different sensitivities of z * and h * to static stability, we do not assume proportionality between the two sets of dimensionless empirical constants: C *TN , C *CN , C *NS -defining z * , and C TN , C CN , C NS -defining h * , and determine all these constants empirically via LES of various types of PBL. It should be noted that we define z * and h * scales as the matching heights of different model layers.
While explicitly depending on PBL stability and scaling as actual log-layer and boundary-layer heights, z * and h * do not necessarily correspond to them.
In the Appendix, we solve analytically Equation (6) in the PBL core (z * < z < h * ); thus determining vertical profiles of wind velocity (U, V) and Reynolds stress ( x = − K M U/ z, y = − K M V/ z). Then, we match them at the height z = z * with conventional formulations for the logarithmic layer: and satisfy the boundary conditions: In such derivation, the combination of boundary and matching conditions just yields the resistance law, Equation (1), with the following formulations of A and B: Here the integrals ∫ĥ 0̂x dẑ and ∫ĥ 0̂y dẑ are expressed as functions of the dimensionless height: Thus A and B are expressed as functions of u * /(fz * ) and h * /z * which, in turn, are specified via Equations 12 and 13 as functions of N and . By this means A and B are defined as functions of N and : Details of the derivations and Equations A8-A11 specifying the functions F A ( N , ) and F B ( N , ) are given in the Appendix.
The system of algebraic Equations 1, 3, 5, A8-A11 allows calculating our unknowns, u * and , given the geostrophic wind speed, U g , roughness length, z 0 , Coriolis parameter, f , Brunt-Väisälä frequency in the free atmosphere, N, and the near-surface buoyancy flux, F bs . The explicit analytical solution to this nonlinear system is generally impossible because stratification parameter involves friction velocity, u * . However, the algorithm linking the surface stress (u * and ) with governing parameters (including U g , N and F bs ) is quite simple and does not cause any difficulties.

VERIFICATION AND CALIBRATION OF THE MODEL
To calibrate the above model, it remains to determine dimensionless constants C *TN , C *CN , C *NS specifying z * in Equation (12) and C TN , C CN , C NS specifying h * in Equation (13). Because of the practical impossibility of controlling all factors affecting the PBL in field experiments (baroclinicity, large-scale vertical velocity, heterogeneity of the flow, etc.), we examine Equations 12 and 13 against data from the LES database (DATABASE64: Esau, 2009) that provides u * , and all the required parameters including F bs and N. This allows determining empirical values of coefficients A and B in Equation (1) and, with the aid of Equations 16 and 17, empirical estimates of virtual parameters of our model: z * and h * .
Available LES data allow us to consider separately the three primitive types of neutral and stable PBL: TN (F bs = 0, N = 0), CN (F bs = 0, N > 0) and NS (F bs < 0, N = 0). For TN PBL, Equations 12 and 13 reduce to.
LES data shown in Figure 1 confirm these formulations and yield C *TN = 0.10, C TN = 0.53.
Next, we consider separately CN PBL for different N and NS PBL for different , so that Equations 12 F I G U R E 1 LES validation of (a) Equation (12) for log-layer height-scale, z * , and (b) Equation (13) for PBL height-scale, h * , for truly neutral PBL when these equations reduce to Equation (19). The lines are plotted after Equation (19) with C *TN = 0.10 and C TN = 0.53. Each black dot represents a single LES experiment of DATABASE64 (Esau, 2009) (12) and (c,d) for nocturnal stable PBL when these equations reduce to Equation (21).

F I G U R E 3 LES validation of (a,b) Equation
Lines are plotted after Equation (20) with C *NS = 0.076 and C NS = 0.97; black dots are LES data from the same sources as in Figure 1 and 13 become: ( LES data shown in Figures 2 and 3 confirm these formulations and yield C *CN = 6.4, C CN = 5.9, C *NS = 0.076, C NS = 0.97. Figures 2 and 3 were checked for hidden correlation due to common parameters in the variables plotted. This revealed there was no significant effect on the accuracy of the calibration. Thus all dimensionless constants of the model are determined. We attribute a bit lower accuracy of our estimates of C *TN and C TN based on LES of TN PBL (Figure 1) to the fact that the simulation run times happened to be insufficient for attaining the steady state. Indeed, in the case of neutral stratification across the entire domain, turbulence is not suppressed at all, which essentially extends the time-scale of PBL relaxation. Moreover, large eddies originated within PBL easily propagate into the free flow and live there quite long, which creates the illusion of levelling off the PBL upper boundary. It is conceivable that longer runs are needed to approach the steady-state and thus refine our estimates of empirical constants C *TN and C TN in Equations 12 and 13. Luckily, for meteorological applications, precise values of these constants are unimportant. Indeed, TN PBLs are never observed in the atmosphere, whereas the terms containing N or in Equations 12 and 13 are at least an order of magnitude larger than those involving C −1 * TN and C −1 TN , so that the latter are practically negligible.
Notably, the above calibration of empirical constants does not use data on LS PBLs allowing us to utilise it for independent validation. But using only one specific LES code for both calibration and validation might raise a question of the general applicability of the model so we performed additional LES experiments of the stably stratified sheared boundary layer with the aid of MSU-INM DNS/LES code developed at the Moscow State University and the Institute of Numerical Mathematics (Mortikov, 2016;Glazunov et al., 2019;Mortikov et al., 2019). In contrast to DATABASE64, the MSU-INM LES code uses the dynamic Smagorinsky eddy viscosity model employing the dynamic approach supplemented with the Lagrangian averaging technique (Bou-Zeid et al., 2005) for both the evaluation of the subgrid momentum and heat fluxes. The baseline experiment corresponded to the GABLS1 case (Beare et al., 2006) while additional experiments had the same set-up with one of the parameters varied. It includes LES runs with varying geostrophic wind magnitude (from 4 to 16 m ⋅ s −1 ), surface cooling rate (from 0.15 to 0.55 K ⋅ hr −1 ), temperature lapse rate above the boundary layer (up to 0.02 K ⋅ m −1 ) and halved Coriolis parameter. All simulations were carried out for 9 hr with a grid resolution of at least ∼3 and ∼1.5 m for cases of stronger surface cooling or weak wind conditions. The last hour of the simulation was used to calculate time-averaged flow statistics following the set-up of the GABLS1 experiment. The MSU-INM LES run of the original GABLS1 experiment shows good agreement with the available estimates of the boundary-layer height, friction velocity and the cross-isobaric angle (Svensson and Holtslag, 2009), thus qualifying the whole dataset as the reliable source of LS PBL data for independent validation. Figure 4, utilising all available data, including LS PBLs by two different LES codes, shows good correspondence between the friction velocity and the cross-isobaric angle calculated through the resistance law (u *mod and mod ) and obtained directly from LES (u * and ). This proves general applicability of the resistance law, in particular, to LS PBLs, that is, beyond the conditions in which our model was calibrated.
Note that the low sensitivity of the cross-isobaric angle to various parameters in the MSU-INM LES runs is consistent with a recent analysis by Grisogono (2011). The observed scatter of values for mod in this case is likely to be explained by the fixed simulation time span of 9 hr irrespective of the parameters varied; different conditions result in different time-scales of the flow adjustment to equilibrium state. In the case of halved Coriolis parameter, the larger angle is observed in LES due to the larger inertial period. The proposed algorithm successfully reproduces this case as well.

CONCLUSIONS
We have derived and empirically validated a comprehensive PBL bulk resistance law, Equations 1 and 18, including long-lived stable PBLs and, thus, covering all currently known types of neutral and stable PBLs: • Truly Neutral with neutral stratification at the surface, F bs = 0, and above PBL, N = 0.
The dimensionless coefficients A and B in Equation (1) are now defined as universal functions of the two parameters of stratification: and N , accounting for static stability in the surface layer and the free atmosphere, respectively. Universal dimensionless constants quantifying our model are determined empirically via LES. The resistance law allows calculating the surface stress (friction velocity, u * , and cross-isobaric angle, ) through parameters of mean flow: geostrophic wind speed, U g , and Brunt-Väisälä frequency in the free troposphere, N, given the surface roughness length, z 0 , and (for NS and LS PBLs) the surface buoyancy flux, F bs < 0.
In our analyses, we have employed a schematic PBL model deliberately simplified in such a way that it yields principally correct analytical formulation of the resistance law, whereas numerical coefficients appeared in the solution are determined empirically -via calibration of the resultant resistance law against wide-ranging LES data. Hence, beyond the resistance law, the model does not claim to be quantitatively correct. As an example, dimensionless empirical coefficients C TN , C CN and C NS in Equation (13) quantifying height scale, h * , for various PBL types are determined to assure the accuracy of the resistance law rather than to estimate PBL height. The latter is already done and yields realistic formulation of the LS PBL height by .
In the real atmosphere, the friction velocity, u * , and the cross-isobaric angle, , as well as PBL height, h, essentially depend on geostrophic wind shear, = | U g / z|, and large-scale vertical velocity at the PBL upper boundary, w h . This is precisely the reason why we did not use field data to calibrate the resistance law. Modifications of the resistance law accounting for and w h can be found in papers by Zilitinkevich and Esau (2003;. Instead of calculating the surface stress given U g , the resistance law can be used to calculate U g (and, thus, the horizontal gradient of atmospheric pressure) given the surface stress. Thanks to existing methods of detecting the surface stress from satellite remote sensing of sea waves, the resistance law may allow retrieving atmospheric pressure over wide areas of ocean poorly covered by in situ observations (Brown and Levy, 1986;Velden et al., 2006;Monzikova et al., 2016).
Notably, our resistance law covers long-lived stable PBLs typical at high latitudes over oceans and, in winter, also over continents. This makes the law potentially quite useful. In particular, it provides the alternative method for calculating the surface stress in very shallow PBLs, e.g. during polar nights, when PBL height could reduce to a few dozen metres , so conventional flux-calculation schemes in weather and climate models become incapable to resolve the surface layer and, hence, to realistically calculate u * and (Holtslag et al., 2013). Besides polar areas, very shallow PBLs are inherent to warm air masses over cold water (a quite typical situation in coastal areas: for example, Von Engeln and Teixeira, 2013, Davy, 2018). In all these cases, the resistance law may be used for calculation of u * and without resolving PBL. Yet, the applicability of the proposed model for such challenging problems as heterogeneous land surfaces or air-ocean interactions within the boundary layer undoubtedly requires further careful examination.  Figure 1; no fitting involved core (z * < z < h * ) are described by Ekman equations, Equation (6). We differentiate Equation (6) with respect to z, multiplying by K M and reformulate the result in terms of Reynolds stress ( x , y ):

How
Then we introduce dimensionless variables: Equation (A1) and boundary conditions (Equations 14 and 15) become (̂x,̂y)|ẑ =0 = (1, 0), (̂x,̂y)|ẑ =ĥ = (0, 0), whereĥ = (h * − z * )∕ √ ku * z * ∕f is the dimensionless height of the PBL core (0 <ẑ <ĥ). The analytical solution to Equation (A3) readŝ where the four dimensionless constants are defined by satisfying boundary conditions (Equation (A4)): Integrating x /K M and y /K M over z yields wind velocity components, U and V, in the PBL core. Then satisfying boundary conditions (Equations 14 and 15) just yields the resistance law, Equation (1), with the following formulation of coefficients A and B: or, taking the integrals ∫ĥ 0̂x dẑ and ∫ĥ 0̂y dẑ: , (A8) Thus A and B are specified as functions of u * /(fz * ) and h = (h * ∕z * − 1)∕ √ ku * ∕(f z * ). Using Equations 12 and 13, u * /(fz * ) and h * /z * are expressed through the stratification parameters and N : Equations A8-A11 define A and B as functions of and N . Figure A1 demonstrates good accuracy for A and B calculated through Equations A8-A11 and obtained directly from LES. As discussed in Section 3 a bit lower accuracy for nearly neutrally stratified cases is of lesser importance for practical applications as TN PBLs are barely observed in the atmosphere. As a result of the aforecited derivation, A and B are defined for all known types of neutral and stable PBLs.