Geostrophic drag law for conventionally neutral atmospheric boundary layers revisited

The geostrophic drag law (GDL), which predicts the geostrophic drag coefficient and the cross‐isobaric angle, is relevant for meteorological applications such as wind energy. For conventionally neutral atmospheric boundary layers (CNBLs) capped by an inversion, the GDL coefficients A and B are affected by the inversion strength and latitude, expressible via the ratio of the Brunt–Väisälä frequency (N) to the Coriolis parameter (f). We present large‐eddy simulations (LES) covering a wider range of N/|f| than considered previously, and show that A and B obtained from carefully performed LES collapse to a single curve when plotted against N/|f|. This verifies the GDL for CNBLs over an extended range of N/|f| within LES. Additionally, in agreement with atmospheric observations, we show that using A = 1.9 and B = 4.4 accurately predicts the geostrophic drag coefficient in the limit of weak inversion strength or high latitude ( N/|f|≲300 ). However, due to the strong dependence of B on N/|f|, corresponding predictions for the cross‐isobaric angle are less accurate. As we find significant deviations between the LES results and the original parameterization of the GDL for CNBLs, we update the corresponding model coefficients.


INTRODUCTION
The atmospheric boundary layer (ABL) is the lowest part of the atmosphere, playing a key role in weather phenomena and human activities. In simplified studies of the ABL, the atmospheric background stratification is neglected and the ABL is assumed to be neutrally stratified. In such a truly neutral atmospheric boundary layer (TNBL), the geostrophic drag coefficient u * /G and the cross-isobaric angle 0 can be determined using the classic geostrophic drag law (GDL, Rossby and Montgomery, 1935;Blackadar and Tennekes, 1968;Tennekes, 1973), where A and B are dimensionless model coefficients, which are presumably universal constants, = 0.4 is the von Kármán constant, u * is the friction velocity, G is the geostrophic wind speed, z 0 is the roughness length, and f = 2Ω sin is the Coriolis parameter. Here, Ω = 0.729 × 10 −4 rad/s is the rotation rate of the Earth and is the latitude. On the right-hand side of the second equality in Equation (1), the minus sign is used for the Northern Hemisphere and the plus sign for the Southern Hemisphere.
The TNBL height h can be estimated by using the following relation (Rossby and Montgomery, 1935;Zilitinkevich et al., 2007): where C R is an empirical dimensionless constant that depends upon the definition of the ABL height. In this work, we use one of the commonly accepted definitions for h, that is, the height at which the total momentum flux reaches 5% of the surface value. We emphasize that other definitions for the ABL height h, which are based on the vertical wind speed, heat flux, or potential temperature profiles, are also commonly used (see, e.g., Abkar and Porté-Agel, 2013;Allaerts and Meyers, 2015;Kelly et al., 2019). Based on this definition and large-eddy simulations (LES) data, Zilitinkevich et al. (2007) recommended C R = 0.6. However, ABLs are seldom truly neutral as there is at least thermal stratification in the free atmosphere. Such ABLs are known as conventionally neutral atmospheric boundary layers (CNBLs, Zilitinkevich and Esau, 2002). Therefore, the effect of free-atmosphere stratification must be taken into account to make the drag law consistent with meteorological observations and numerical simulations (Hess and Garratt, 2002a;2002b;Zilitinkevich and Esau, 2002;Hess, 2004). Zilitinkevich and Esau (2005) generalized the classic GDL to the conventionally neutral and stable ABL regimes. They showed that the GDL is still expressed by Equation (1), while the dimensionless functions A and B are not universal constants but functions of the internal stability parameter = u * ∕(|f |L s ), the external stability parameter N = N∕|f |, and the baroclinicity parameter G = (dG∕dz)∕N. Here, L s is the surface Monin-Obukhov length, N is the Brunt-Väisälä frequency, and dG/dz is the baroclinic shear. In this paper we focus on the CNBL case with zero heat flux at the surface and uniform geostrophic wind above the boundary layer such that = 0 and G = 0. As a result, N , which is named the Zilitinkevich number by Esau (2004), is the only free parameter. Based on similarity theory, alternative limits, and asymptotic matching, Zilitinkevich and Esau (2005) derived the following general form for the dimensionless coefficients: where a and b are dimensionless constants obtained from similarity theory, a 0 and b 0 are correction constants to be determined empirically, and m is the composite stratification parameter (or equivalently, the ratio of the CNBL height to the generalized turbulence length scale) originating from the alternative limits. It should be pointed out that, in the original version of the GDL for CNBLs, that is, Equation (3), Zilitinkevich and Esau (2005) assigned different composite stratification parameters for A and B, namely m A and m B , each with two additional empirical constants, that is, (C fA , C NA ) and (C fB , C NB ), corresponding to the rotational length scale L f = u * /|f | and the external static-stability length scale L N = u * /N, respectively. Based on their LES dataset, Zilitinkevich and Esau (2005) determined C fA = C fB = 1, C NA = 0.09 and C NB = 0.15. In this paper, we still assume C fA = C fB = 1 but combine C NA and C NB into a single coefficient, namely C m = 0.1 (see below). This treatment is equivalent to assuming the same turbulent length scale in both the x-and y-components of the GDL for CNBLs and thus m A = m B = m.
Following the derivation of Zilitinkevich and Esau (2005), the composite stratification parameter m can then be written as where C m is an empirical constant, and the external stability parameter or Zilitinkevich number (Esau, 2004) is Here Γ is the free-atmosphere lapse rate, g is the acceleration due to gravity, and 0 is the reference potential temperature. Similarly, accounting for the alternative limits, the CNBL height h can be estimated as (Zilitinkevich et al., 2007) where C R and C N are empirical constants. Note that, for the TNBL ( N = 0), Equation (6) reduces to Equation (2). Two parameters govern the physics in stratified boundary layers: one that characterizes the size, speed, or scale separation, and another that describes the degree of stratification. Therefore, the dynamics can be described by different parameters such as a flux or gradient Richardson number or a buoyancy Reynolds number (Coleman et al., 1992;Flores and Riley, 2011;van Hooijdonk et al., 2018;Sun et al., 2020). However, in this study, we use the Rossby number Ro = G/(z 0 |f |) and the Zilitinkevich number N = N∕|f |, as they were used in the GDL for CNBLs introduced by Zilitinkevich and Esau (2005). Previous LES studies on the GDL for CNBLs mainly focused on high-latitude cases (Zilitinkevich and Esau, 2002;, Zilitinkevich et al., 20072012). However, the applicability of these findings for a broader range of Zilitinkevich numbers is of fundamental interest and relevant for meteorological applications such as wind energy (Kelly et al., 2019). The A and B values from the GDL are mainly used to extrapolate the wind profile beyond the surface layer to obtain estimates for the available wind resources (Gryning et al., 2007;Kelly and Gryning, 2010;Kelly and Troen, 2016). Therefore, it is essential to characterize the integral properties of the CNBL over a wide range of N , that is, for Γ from 1 to 10 K⋅km −1 (Sorbjan, 1996;Allaerts and Meyers, 2015), and for from low to high latitudes. The goal of the present study is to determine a and b such that the asymptotic behavior of A and B for the relative limit N C 2 R ∕C 2 N ≫ 1 is faithfully captured over a wider Zilitinkevich number range than considered previously. Our new simulation campaign shows that there are significant differences between the original Zilitinkevich and Esau (2005) parametrization and the results from LES, especially at high Zilitinkevich numbers. We show that this finding does not depend on the employed subgrid scale (SGS) model or the used numerical resolution and propose a new parametrization based on the available LES data.
In this study, we systematically assess the effect of the latitude and the free-atmospheric lapse rate on the integral properties of the CNBL by performing LES. In Section 2, we discuss the numerical method and simulation setup. In Section 3, we determine the empirical constants in the GDL for CNBLs for a wide parameter range. Furthermore, we present a comparison of our LES results with the available atmospheric data and theoretical predictions. We conclude with a summary of the main findings in Section 4.

Governing equations
We use LES to simulate the CNBL flow over an infinite flat surface with homogeneous roughness. We integrate the spatially filtered Navier-Stokes equations and the filtered transport equation for the potential temperature (J. D. Albertson, 1996;Albertson and Parlange, 1999;Gadde et al., 2020): t̃+ũ ⋅ ∇̃= −∇ ⋅q.
Here, the tilde denotes spatial filtering, ⟨⋅⟩ represents horizontal averaging,ũ is the velocity,̃= ∇ ×ũ is the vorticity,p is the modified pressure departure from equilibrium, is the potential temperature, = g∕ 0 is the buoyancy parameter, G is the geostrophic wind velocity, denotes the deviatoric part of the SGS shear stressũu −ũũ (Wyngaard, 2004), and q =ũ −ũ̃represents the SGS heat flux. Viscous terms on resolved scales are neglected as the Reynolds number in the ABL flow is very high.
For closure of Equations (7) and (8), we parameterize the SGS shear stress and heat flux as where t and are the eddy viscosity and eddy diffusivity, respectively, andS is the resolved strain-rate tensor with the superscript T denoting a matrix transpose. We use the recently developed anisotropic minimum dissipation (AMD) model (Abkar and Moin, 2017) to perform the simulations. In the AMD model, the eddy viscosity and the eddy diffusivity are determined such that the energy of the subfilter scales of the LES solution does not increase with time (Abkar and Moin, 2017). We note that the AMD model has been extensively validated and has been shown to give similar predictions to the well-known Lagrangian-averaged scale-dependent SGS model (Bou-Zeid et al., 2005;Gadde et al., 2020). To confirm that the main findings do not depend on the employed SGS model, we compared the results for the lowest and highest N with results obtained using the standard Smagorinsky model (Smagorinsky, 1963;Mason and Thomson, 1992;Liu and Stevens, 2020).
As the main goal of this study is to characterize the integral properties of the CNBL case and extend the parameter range for which the GDL for CNBLs is validated, we employ the same assumptions as used in the previous studies (Zilitinkevich and Esau, 2005;Abkar and Porté-Agel, 2013;Kelly et al., 2019). We note that the f -plane approximation is one of the underlying assumptions of the GDL for CNBLs (Zilitinkevich and Esau, 2005). Therefore, the f -plane approximation (i.e., without considering the horizontal component of the rotation of the Earth) is used in all the presented simulations. This allows us to study an extensive range of latitudes ∈ [5 • , 70 • ], corresponding to 42 ≤ N ≤ 1,350 with a lapse rate of Γ = 1 ∼ 9 K⋅km −1 . It should be pointed out that previous direct numerical simulations at low Reynolds number (Coleman et al., 1990) showed that the f -plane approximation has a significant effect on u * and 0 . To assess the suitability of the f -plane approximation for this study, we performed simulations with and without this approximation. We find that, for high-Reynolds-number flows, the effect of this approximation is relatively limited; that is, at a latitude of 5 • , the effect on u * is less than 2%.

Numerical method
Our code is an updated version of the one used by Albertson and Parlange (1999). The grid points are uniformly distributed, and the computational planes for horizontal and vertical velocities are staggered in the vertical direction. The first vertical velocity grid plane is located at the ground, and the first streamwise and spanwise velocities and potential temperature grid planes are located at half a grid distance away from the ground. In the vertical direction, we use a second-order finite difference method, while we use a pseudospectral discretization with periodic boundary conditions in the horizontal directions. Time integration is performed using the second-order Adams-Bashforth method. The projection method is used to ensure the divergence-free condition of the velocity field.
At the top boundary, we enforce a constant potential temperature lapse rate, zero vertical velocity, and zero shear stress boundary condition. At the bottom boundary, we employ the classical wall stress formulation based on the Monin-Obukhov similarity theory for the ABL (Moeng, 1984;Bou-Zeid et al., 2005), where the overline denotes filtering at a scale of 2Δ with Δ as the filter scale of Equations (7), z 1 = Δz∕2 is the half vertical grid distance, and xz and yz are the shear stresses in the streamwise and spanwise directions, respectively. No stability correction functions are needed in the wall model because the surface heat flux is zero; that is, the boundary layer is neutrally stratified. The 2Δ filtering was suggested by Bou-Zeid et al. (2005) to reduce the log-layer mismatch (Brasseur and Wei, 2010). The rationale for this procedure is that the filtering preserves the important large-scale variations. At the same time, the use of filtered velocities results in an average stress that is very close to the stress predicted by the average similarity formulation for homogeneous surfaces (Bou-Zeid et al., 2005). The vertical derivative at the first horizontal plane is calculated using the Monin-Obukhov similarity theory (Moeng, 1984), which is implemented following the method proposed by J. D. Albertson (1996).

Computational setup
The computational domain size in our simulations is 2 km × 2 km × 2 km in the streamwise, spanwise, and vertical directions, respectively. The horizontal domain size is at least six times larger than the boundary layer height to ensure that large horizontal streamwise flow structures are captured appropriately for all cases (Foster, 2013). The positive z-direction is opposite to the gravity vector, the positive x-direction is aligned with the geostrophic wind direction at the top of the domain, and the y-direction is defined such that (e x , e y , e z ) forms an orthogonal coordinate system. The initial potential temperature profile is (z) = 0 + Γz, where 0 = 300 K is the reference potential temperature and Γ = 1, 3, 9 K⋅km −1 is the free-atmosphere lapse rate. The initial velocity profile is u = Ge x , where G = 12 m/s. Small random perturbations are added to the initial fields of u and near the surface (z ≤ 100 m) to spin up turbulence. As was pointed out by Pedersen et al. (2014) and Kelly et al. (2019), the used initial profile is not particularly relevant when one is concerned with the CNBL characteristics in the quasiequilibrium state. The roughness length is z 0 = 10 −4 m, which is a typical value for the sea surface (Hess and Garratt, 2002a). Statistics are collected when the boundary layer has reached a quasistationary state. To be specific, data are collected over the nondimensional time interval ft ∈ [9, 10]. At high latitudes (e.g., = 70 • ), this is equivalent to simulating the flow for approximately 20 physical hours (Pedersen et al., 2014). However, for very low latitude (e.g., = 5 • ), it is equivalent to about nine days. The main simulations are performed using the AMD model on a 288 3 grid, such that the resolution is 21.8 m in the horizontal directions and 6.9 m in the vertical direction. A summary of these simulations is presented in Table 1. To show that the main findings do not depend on the used SGS model, we performed simulations for N = 42 and 1350 using the Smagorinsky model (Table 2). In addition, we performed simulations on different grid resolutions, that is, 72 3 , 144 3 , and 288 3 , for the same value of N (Table 3), to confirm that the main findings do not depend on the employed grid resolution. Figure 1 shows that the temporal evolution of the domain-averaged (a) friction velocity u * and (b) cross-isobaric angle 0 for N = 42 and N = 1350 obtained on different resolutions are very similar. Besides, the figure reveals that the quasiequilibrium state has been established at ft ≥ 9 when time averaging is started. Figure 2 shows the time-averaged wind speed profile as a function of height for (a) N = 42 and (b) N = 1,350. The figure shows that, regardless of N and the adopted grid resolution, the simulated velocity profiles closely follow the theoretical logarithmic TA B L E 1 Summary of AMD simulations, performed on a 288 3 grid. Note that the GDL coefficients A and B obtained from the simulations are determined using Equation (1 curve approximately up to 50 m, which confirms that the log-layer mismatch in our simulations is minimal. This shows that our results are in the high-accuracy zone prescribed by Brasseur and Wei (2010), which is required to obtain reliable predictions in wall-modeled LES of boundary layers. This is confirmed by the results for A and B obtained from the simulations shown in Figure 5, which reveals a perfect collapse of all simulation results when plotted versus N while the effect of the grid resolution and the used SGS model on the results is limited.   Figure 3 shows the vertical profiles of the time-averaged (a) wind speed and (b) potential temperature at = 70 • for different lapse rates. To minimize the effect of slowly decaying inertial oscillations in the outer layer, the profiles are averaged over one inertial period Δ(ft) = 2 (Spalart, 1989;Coleman et al., 1992;Coleman, 1999;Peña et al., 2014;Jiang et al., 2018). A distinct feature of the CNBL is that the potential temperature is almost uniform below the inversion layer. In agreement with Abkar and Porté-Agel (2013), we observe that, for the above-described uniform initial temperature profile (Pedersen et al., 2014), the CNBL depth reduces with increasing lapse rate (Figure 3, Table 1). However, the prediction of the velocity profile shown in Figure 3a is nontrivial, but of great fundamental interest and relevant for engineering applications such as wind energy (e.g., Kelly et al., 2019). Therefore, many researchers have made efforts in this direction (see, e.g., Blackadar, 1962;Gryning et al., 2007;Esau, 2004;Abkar and Porté-Agel, 2013;Jiang et al., 2018;Kelly et al., 2019). In particular, Kelly et al. (2019) derived a model that accounts for the effect of the Coriolis parameter and free-atmosphere lapse rate to obtain an analytical expression that captures the wind speed profiles from LES up to 90% of the boundary layer height with a relative error of less than 5%. Figure 4 shows the variation of the dimensionless CNBL depth |f |h/u * versus the Zilitinkevich number N = N∕|f |. We find excellent agreement between the LES data (open symbols) and the theoretical model predictions (solid line) given by Equation (6) with C R = 0.5 and C N = 1.6, confirming that the dimensionless depth |f |h/u * only depends on the Zilitinkevich number N . LES data obtained using the Smagorinsky model on 288 3 grids (filled diamonds) and using the AMD model on 72 3 and 144 3 grids (filled hexagram) agree well with the high-resolution data, which indicates that the results do not depend on the grid resolution or the employed SGS model. The prediction obtained by using Equation (6) with C R = 0.6 and C N = 1.36 ± 0.25 as suggested by Zilitinkevich et al. (2007), as well as the field measurements taken from Zilitinkevich et al. (2012), are shown for comparison. Besides, note that Zilitinkevich and Esau (2003) obtained C R = 0.5 and C N = 1.1 ∼ 2.2 and Abkar and Porté-Agel (2013) concluded C R = 0.5 and C N = 1.5, both of which are consistent with our results. Figure 5 shows the A and B values obtained from the simulations as a function of the Zilitinkevich number N . The values of A and B (open symbols) are obtained from the simulation results using Equation (1); see also Table 1. The fact that all the LES data (open symbols) collapse onto a single curve confirms that the GDL for CNBLs is suitable to describe the dynamics in the entire Zilitinkevich number range under consideration. The figure also shows the F I G U R E 4 Dimensionless CNBL height |f |h/u * versus Zilitinkevich number N = N∕|f |. Open symbols: LES data. Filled stars: field data taken from Zilitinkevich et al. (2012); solid line: theoretical curve given by the present work using Equation (6) with C R = 0.5 and C N = 1.6; dashed line and shaded region: theoretical prediction given by Zilitinkevich et al. (2007) with C R = 0.6 and C N = 1.36 ± 0.25. Filled diamonds: LES data obtained on a 288 3 grid using the Smagorinsky model. Filled hexagram: LES data obtained using the AMD model on a 72 3 and 144 3 grid LES data obtained on a 288 3 grid using the Smagorinsky model (filled diamonds) and using the AMD model on 72 3 and 144 3 grids (filled hexagram). The simulation results obtained on coarser meshes or with the Smagorinsky model agree well with the main high-resolution simulation results obtained using the AMD model. These tests confirm that our simulation results capture the integral boundary layer properties more accurately than previous studies, and for a broader range of latitudes and free-atmosphere stratification than considered previously. It is worth mentioning that the simulation challenges result from the intrinsically developing nature of the CNBL case in combination with the wall-modeled LES that needs to be adopted. Therefore, it is not possible to show the absolute convergence of the simulations results with increasing grid resolution. We note this is most challenging for high Zilitinkevich number cases, but as shown

F I G U R E 5 Dimensionless coefficients (a) A and (b) B versus
Zilitinkevich number N = N∕|f |. Open symbols: LES data using Equation (1). Solid line: theoretical curve given by the present work using the GDL for CNBLs (Equations (3)-(6)) with the empirical constants given by Equation (11). Dashed line: theoretical prediction given by Zilitinkevich and Esau (2005). Filled diamonds: LES data obtained from the Smagorinsky model with a 288 3 grid. Filled hexagram: LES data obtained with the AMD model on a 72 3 and 144 3 grid in Figure 5, the convergence of the simulation results is still good in that limit.
However, unfortunately, we find significant differences between the LES results (open symbols) and the predictions from the GDL for CNBLs, see Equations (3)-(6), using the Zilitinkevich and Esau (2005) parametrization (a = 1.4, a 0 = 1.65, b = 10, b 0 = −2, C m = 0.1; see dashed lines). Their fit was obtained from state-of-the-art LES data with resolutions up to 64 3 available at that time, which unfortunately show significant scatter. Therefore, we use the data presented in Figure 5 to reevaluate the constants used in the GDL parametrization for CNBLs. As remarked in Section 1, we assume the same composite stratification parameter m for both A and B, thus a single coefficient C m is sufficient to determine the value of m, provided that C R = 0.5 and C N = 1.6 are accurately determined ( Figure 4). We determine the values for a and b such that the asymptotic behavior of A and B in the high Zilitinkevich number limit is well captured. Subsequently, the correction constants a 0 and b 0 are set such that A and B also capture the low N limit well. The resulting constants in Equations (3) and (4)  In meteorological applications, the GDL is employed to predict the geostrophic drag coefficient u * /G and the cross-isobaric angle 0 . Typical values of A and B reported from CNBL meteorological observations that are used to predict u * /G and 0 are in the range A = 1.3 ∼ 1.9 and B = 4.4 ∼ 4.9 (see, for example, Troen and Petersen, 1989;Hess and Garratt, 2002a;Zilitinkevich and Esau, 2002;Gryning et al., 2007;Kelly and Troen, 2016). In particular, Hess and Garratt (2002a) recommended A = 1.3 and B = 4.4, Gryning et al. (2007) adopted A = 1.9 and B = 4.9, and Kelly and Troen (2016) used A = 1.8 and B = 4.5. Therefore, it is insightful to see how our LES results compare with these findings. Figure 6 compares the (a) cross-isobaric angle 0 and (b) geostrophic drag coefficient u * /G obtained from the simulations with the model predictions using Equation (1) with A = 1.9 and B = 4.4. The coefficients A and B have been obtained by a fit of Equation (1) to the LES data. This figure indeed confirms that, regardless of the Rossby number Ro, these constants provide a good approximation for geostrophic drag coefficient u * /G for N ≲ 300, which is typical for many atmospheric observations in the limit of weak inversion strength or high latitude (Troen and Petersen, 1989;Hess and Garratt, 2002a;Zilitinkevich and Esau, 2002;Gryning et al., 2007;Kelly and Troen, 2016). However, due to the strong dependence of B on N , the corresponding predictions for the cross-isobaric angle are less accurate. The figure shows that using the updated parametrization, that is, Equations (3)- (6) and (11), allows one to accurately predict both the geostrophic drag coefficient u * /G and the cross-isobaric angle 0 over the considered parameter range. Note that Figure 6 shows that the predictions for the geostrophic drag coefficient u * /G obtained using Equation (1) and the A and B values obtained using the GDL for CNBLs (Equations (3)- (6) and (11)) exhibits a maximum at around N ≈ 100 (solid lines). An intuitive conjecture for the origin of the peak is that it is caused by the −am term in Equation (3) becoming dominant over the ln(a 0 + m) term. However, this does not seem to be the case, as the peak is still visible when this term is neglected by setting −am ≡ 0. Unfortunately, the data obtained from LES are not sufficiently accurate to confirm or rule out the existence of this very small peak.

3.2
Effect of latitude Figure 7 shows the vertical profiles of the time-averaged wind speed and potential temperature for the lapse rate Γ = 1 K/km. At higher latitude, we observe a shallow CNBL with a very slightly reduced wind speed maximum ( Figure 7a) and a lower potential temperature below the inversion layer ( Figure 7b). Interestingly, Γ has a negligible effect on the friction velocity u * , while has a strong influence on u * . This indicates that the thermal stratification mainly affects the shape of the wind speed profile in the outer region of the CNBL (Figure 3a). In contrast, the latitude significantly affects the shape of the wind speed profile throughout the entire CNBL (Figure 7a).  Tables I, II, and V of Hess and Garratt (2002a) for the roughness lengths z 0 ∈ [0.05, 0.5] mm, since the roughness length adopted in our simulations is z 0 = 0.1 mm. The theoretical curves are given by Equations (3)-(6) with the empirical constants given by Equation (11) (solid lines) and Zilitinkevich and Esau (2005) (dashed lines). The atmospheric field data for A and B show the same trends as the simulation results and theoretical predictions. We note that the scatter in the field data is caused by measurement uncertainties.

CONCLUSIONS
We used large-eddy simulations (LES) to systematically study the effect of the free-atmospheric stratification and latitude on the integral measures of conventionally neutral atmospheric boundary layers (CNBLs) for an extensive range of Zilitinkevich numbers N ∈ [42,1350]. We find that the CNBL depth h reduces significantly with increasing inversion strength (i.e., lapse rate Γ) and latitude . We find that the dimensionless CNBL height |f |h/u * only depends on the Zilitinkevich number N and can be predicted well by the empirical formula (6) using the constants C R = 0.5 and C N = 1.6. We show that the dimensionless coefficients A and B calculated from the LES data using Equation (1) collapse to a single curve when plotted versus N . This extends the range for which the geostrophic drag law (GDL) for CNBLs is verified within LES, while it shows the ability of LES to accurately capture such boundary layer dynamics. We find that the original Zilitinkevich and Esau (2005) parametrization, obtained from LES with resolutions up to 64 3 , does versus . Open symbols: LES data using Equation (1). Filled symbols: field data taken from Hess and Garratt (2002a); thick lines: theoretical curves given by the present work using the GDL for CNBLs (Equations (3)-(6)) with the empirical constants given by Equation (11); thin lines: theoretical curves given by Zilitinkevich and Esau (2005) not properly describe the simulation results. Therefore, we propose an updated parameterization of the GDL for CNBLs, which accurately describes the geostrophic drag coefficient u * /G and the cross-isobaric angle 0 over the extended Zilitinkevich number range considered in this study. We further show that the combination A = 1.9 and B = 4.4, which is in agreement with results from CNBL meteorological observations (Troen and Petersen, 1989;Hess and Garratt, 2002a;Zilitinkevich and Esau, 2002;Gryning et al., 2007;Kelly and Troen, 2016), provides a good estimate for the geostrophic drag coefficient u * /G for N ≲ 300. However, due to the strong dependence of B on N , the corresponding predictions for the cross-isobaric angle are less accurate.