Parameterization of a surface drag coeﬀicient in conventionally neutral planetary boundary layer

. Modern large-scale models (LSMs) rely on surface drag coefﬁcients to parameterize turbulent exchange be-tween surface and the ﬁrst computational level in the atmosphere. A classical parameterization in an Ekman boundary layer is rather simple. It is based on a robust concept of a layer of constant ﬂuxes. In such a layer (log-layer), the mean velocity proﬁle is logarithmic. It results in an universal dependence of the surface drag coefﬁcient on a single internal non-dimensional parameter, namely the ratio of a height within this layer to a surface roughness length scale. A realistic near-neutral planetary boundary layer (PBL) is usually much more shallow than the idealized Ekman layer. The reason is that the PBL is developing against a stably strat-iﬁed free atmosphere. The ambient atmospheric stratiﬁca-tion reduces the PBL depth and simultaneously the depth of the log-layer. Therefore, the ﬁrst computational level in the LSMs may be placed above the log-layer. In such a case, the classical parameterization is unjustiﬁed and inaccurate. The paper proposes several ways to improve the classical parameterization of the surface drag coefﬁcient for momentum. The discussion is focused on a conventionally neutral PBL, i.e. on the neutrally stratiﬁed PBL under the stably stratiﬁed free atmosphere. The analysis is based on large eddy simulation (LES) data. This data reveals that discrepancy between drag coefﬁcients predicted by the classical parameterization and the actual drag coefﬁcients can be very large in the shallow PBL. The improved parameterizations provide a more accurate prediction. The inaccuracy is reduced to one-tenth of the actual values of the coefﬁcients.


Introduction
Despite of the complexity of the Earth's surface, widely used parameterizations of the turbulent exchange in the sur-Correspondence to: I. N. Esau (igore@nersc.no)face layer generally remain rather simple.The most of the large-scale models (LSMs) participating in the atmospheric model intercomparison program (AMIP) rely on the Monin-Obukhov similarity theory (Monin and Obukhov, 1954) in their calculations for a surface drag coefficient, C D , where τ is a turbulent vertical momentum flux, and u is flow velocity at the first model level.This approach has been developed from studies of von Karman (1975) and Prandtl (1932).Within a surface layer, i.e a layer of constant turbulent fluxes, in terms of the firstorder turbulence closure (Holt and Raman, 1988)  (2) The constant in Eq. ( 2) is simply a turbulent flux at the surface, τ (z=0)=u 2 * , where u * is so-called friction velocity and |u| is the wind speed.Dimensional analysis suggests that K m =l•u s is a combination of length, l, and velocity, u s scales.von Karman proposed l=κz and u s =u * on the basis of laboratory experiments with a well-established layer of constant turbulent fluxes.Here, z is the height above the surface.The constant κ=0.41 is known as the von Karman constant.
Integration of Eq. ( 2) gives an expression for the logarithmic velocity profile in the surface layer where z 0 is surface roughness.The log-law is an essential part of C D -parameterizations.Since the atmospheric PBL is always stratified, Monin and Obukhov suggested a universal stability correction of Eq. (3) in the following form where L=−u 3 * /F bs is the Monin-Obukhov length scale and F bs is a buoyancy flux at the surface.(z/L) is an empirical function, which is not defined in the theory.Empirical essence of (z/L) has resulted in a great variety of possible forms of (z/L) (Abdella and McFarlane, 1996).However, historically first expressions proposed by Businger et al. (1971), Dyer (1974) and Webb (1970) still remain the most popular.The Monin-Obukhov similarity theory can be developed to the C D -parameterization in two ways.The first way follows from the straightforward application of Eq. ( 4) as (5) The second way is to redefine (z/L) as F −2 M =(1− (z/L)/ ln(z/z 0 )), where F M could be seen as a new empirical function.Louis (1979), Abdella and McFarlane (1996) and many others considered F M as a function of the Richardson number, Ri.It gives Here, the classical expression for a drag coefficient for momentum in the Ekman boundary layer reads (7) Lange et al. (2004) demonstrated a systematic discrepancy between C D predicted by Eq. ( 5) and measured in the PBL over the Danish Baltic Sea during the Rodsand measurement program.The discrepancy remains significant even in the near-neutral PBL at large L or small Ri.The discrepancy was found for a wide range of wind speed and at different distances from the sea shore.This work will show that it is necessary to account for an effect of stability of the atmosphere above the PBL.This effect has been overlooked in the above equations.Therefore, it is still absent in the large-scale models.Already Csanady (1974) concluded from his theoretical considerations that the effect should be significant in shallow PBLs.However, he obtained only asymptotic solutions.He only mentioned that in the absence of experimental evidence it is difficult to speculate on the intermediate variations of C D .This work and Lange et al. (2004) provide such evidences.
The stratification of the atmosphere above the PBL does not depend on the stability of the PBL.Therefore, it is reasonable to account for this effect not through F M (Ri) or (z/L) but through C Dn , as it has been done by Csanady.Such an approach will considerably alleviate the following verification.
Until very recently, the near-neutral atmospheric PBL has been considered as a truly neutral or Ekman boundary layer (Zilitinkevich and Esau, 2002).The truly neutral boundary layer was seen to be similar to the Ekman layer in laboratory experiments and numerical simulations (Hess and Garratt, 2002).The only non-dimensional governing parameter in the Ekman layer is the Rossby number, Ro=U g /(f z 0 ), where U g is a geostrophic wind speed and f is the Coriolis parameter.Using scaling analysis and the von Karman scaling, Rossby and Montgomery (1935) proposed an expression for the Ekman layer depth where C h is a constant.The friction velocity, u * , can be found from the following non-linear equation (Zilitinkevich, 1989) where C g =u * /U g is a geostrophic drag coefficient and A is a constant.DNS (Coleman, 1999) and LES (Mason and Thompson, 1987) 8) with C h =0.65.In this case, z LSM 1 would always be placed well within the log-layer of the depth h s .However, both coarse and fine resolution LSMs would have problems in the atmospheric conventionally neutral PBLs.The conventionally neutral PBL is a PBL developing against a stably stratified free atmosphere.Such PBLs are usually shallow, so the level z LSM 1 is often placed above the log-layer or even above the entire PBL.
This study addresses the problem of shallow logarithmic layers in the parameterization of the surface drag coefficients.It is worth noting that several attempts have been made to incorporate the stability effects of the free atmosphere into the C D -parameterization.For instance, Fairall et al. (1996), following Schumann (1988), proposed to use an effective wind, S=(|u| 2 +Cw 2 * ) 1/2 , in Eq. ( 1).Here, w * =(F bs H ) 1/3 is a convective velocity and C is a constant.The convective velocity is proportional to the PBL depth.This parameterization recovers Eq. ( 1) in the near-neutral case, F bs →0.Cassano and Parish (2001) proposed to use an efficient surface roughness for temperature in stably stratified PBLs.In order to account for non-local effects, z 0T z 0 was proposed.It should account for an additional surface drag due to internal waves in the PBL.Such waves become evanescent in the near-neutral PBL.Therefore, the waves cannot exert significant drag in this case.Zilitinkevich et al. (2002) suggested to add an additional term to the log-law in the stably stratified PBL.This term accounts for the stability of the free atmosphere above the PBL.The method reminds the Csanady approach and will be followed in this paper.
Contrary to the above mentioned papers, this paper is focused on the C Dn -parameterization in the conventionally neutral PBL, i.e. the PBL with (i) a negligible buoyancy flux at the surface, F bs →0, (ii) near-neutral stratification in the surface layer, Ri g (z<h s )→0, and (iii) considerable stratification of the free atmosphere above the PBL, Ri g (z>H ) =0. Advanced parameterizations for the non-neutral (convective and stably stratified) PBLs can be derived by substitution of the obtained C Dn expressions into Eqs.( 5) and (6).
Section 2 describes relevant physical processes and features of turbulence within the conventionally neutral PBL.Section 3 describes the method and problems related to calculation of C Dn from a LES database.Section 4 discusses several ways to account for the non-local effects in an improved C Dn -parameterization.Section 5 outlines the conclusions of this study.

Relevant physical processes in conventionally neutral PBLs
In the Ekman boundary layer, turbulent kinetic energy (TKE) is mainly in a local balance with dissipation.The excessive TKE is generated at the top of the surface layer and then transported downward and upward.This process maintains the surface layer as the layer of constant fluxes with the logarithmic velocity profile.The planet rotation has a negligible effect on small turbulent eddies within the surface layer.However, larger eddies above the surface layer experience a significant damping effect due to the planet rotation.
The damping maintains an equilibrium PBL depth, H , given by Eq. ( 8).The equilibrium PBL depth is defined by the balance between the excessive TKE production and the excessive TKE dissipation due to the planet rotation (Tritton, 1992).In the case of the conventionally neutral PBL, the excessive TKE production is consumed both by the planet rotation and by the ambient stratification of the free atmosphere above the PBL.There are two competing mechanisms of the turbulent exchange which are sensitive to the ambient atmospheric stratification.The most obvious mechanism is gradual mixing due to the entrainment of potentially warm air from the free atmosphere.It results in a gradual deepening of the PBL.Another process is radiation of internal waves from the PBL top to the free atmosphere.The waves take out the TKE but leave temperature variations at the PBL top.Therefore, the wave radiation works against the mixing.It increases the temperature gradient across the PBL top (Zilitinkevich, 2002).These two processes result in an equilibrium PBL depth.Csanady (1974) noticed that the equilibrium depth depends on the ambient stratification of the free atmosphere.LES data allowed Zilitinkevich and Esau (2002) to derive a quantitative measure of this dependence.The stratification can be characterized by a non-dimensional Zilitinkevich number, Z n =µ N =N/|f | (Zilitinkevich and Calanca, 2000), where N is the Brunt-Väisälä frequency in the atmosphere above the PBL.The LES studies have shown that the PBL depth decreases more than 10 times from the case with Z n =1 to the case with Z n =350 (see Fig. 1).The turbulent stress at the surface also decreases (see Fig. 2).The reason is that large eddies become limited in size by the PBL depth and therefore less energetic.Zilitinkevich and Esau (2002) generalized Rossby and Montgomery (1935) expression for H , which includes Z n as an additional external governing parameter.This expression reads , where Here, C R =0.65±0.08 and C 0 =0.23±0.05are empirical constants obtained from the LES database.The ambient stability is actually the controlling factor for H in the atmospheric range of N ∈[3•10 −3 ; 3•10 −2 ] s −1 .This is true for virtually any realistic surface heat flux, i.e. |L|>10 2 .The finding has been recently confirmed in the analysis of atmospheric data by Hess (2004).
Summing up, the near-neutral atmospheric PBLs are usually much more shallow than it follows from the theoretical analysis of the idealized Ekman layer.The reason for this is considerable stratification of the free atmosphere immediately above the PBL.This stratification was neglected in the earlier theoretical works.The first computational level in the modern LSMs is usually placed outside the surface layer or even outside the entire shallow PBL.It suggests that the C Dn -parameterization given by Eq. (1) and Eq. ( 7) does not properly represent the surface drag coefficient in the majority of LSMs.The error in C Dn leads to considerable errors in C D according to Eqs. ( 5) and ( 6).

Surface drag coefficient in large eddy simulations
LES is a powerful technique to study turbulent exchange in the PBL.The detailed description of the LES code and the LES database is out of scope of this paper.It can be found in Esau (2004).
However, several important aspects of the LES technique should be highlighted here.The LES technique resolves only relatively large scales of motions.The author's LES code, LESNIC, maintains the statistically correct amount of the TKE down to the first computational level.The price is a noticeable overamplification of the shortest resolved scales at the first 2-3 computational levels.An estimated sub-grid TKE at the first computational level is about 30% of the resolved TKE.The sub-grid part of the TKE decreases rapidly with height.Distinct to the LSMs, the LES resolves the three-dimensional (3-D) structure of turbulence.There is a direct energy cascade in 3-D turbulence.Since the upper part of the Kolmogorov inertial subrange of scales is resolved, small-scale turbulence has only a minor effect on the resolved turbulence.This is easy to see in comparisons of LES with different sub-grid closures and mesh resolutions (Esau, 2004).
The boundary conditions in the LES code are prescribed locally in the form of the log-law, as in Eq. (3).Thus, any LES studies of the surface layer are unavoidably affected by the log-law boundary conditions in the LES code itself.The influence of the boundary conditions, however, decreases very rapidly with height.It follows from analysis of non-dimensional gradients.Application of the log-law in the form of Eq. (3) supplies incorrect boundary conditions in the case of stable and convective PBLs.In the case of the stably stratified PBL, the Monin-Obukhov similarity theory suggests the following expression for the non-dimensional velocity gradient Here, L=−u 3 * /F bs is the Monin-Obukhov length scale, F bs is a surface buoyancy flux and C u is supposed to be a constant.Figure 3 shows C u defined from the author's LES and C u defined from the atmospheric data sets.Recall that the LES code employs Eq. ( 3) but not Eq.( 11) to account for the velocity at the first computational level.The discrepancy between the Monin-Obukhov theory and LES is clearly seen at small z/L, i. e. at the very first computational levels in LES.However, LES demonstrates rather good agreement with the Monin-Obukhov already at the second and higher computational levels, where z/L>0.25.Mayor, Tripoli and Eloranta (2003) have independently assessed the LES quality for the convective PBL.They compared 3-D turbulence structures simulated by a coarse resolution LES and obtained by a volume imaging lidar.To the authors' surprise, they had to conclude that LES reproduced reasonable structures in the surface layer where "the technique is expected to perform poorly".Therefore, except the first and perhaps the second computational levels, the effect of surface layer parameterization on LES data is rather small.
The vertical resolution of the LES is much better than the resolution of the LSMs.The LES runs have the height of the first computational level, z LES 1 , between 19.5 m in the deepest truly neutral PBL (Z n =0, H =2000 m) and 1.5 m in the PBL, capped by the strongest inversion (Z n =400, H =60 m).All LES runs have about 45 levels within the PBL and, therefore, 4-5 levels within the surface layer.Here, the LES mesh was almost isotropic, as the technique of 3-D turbulence resolving modeling requires.In most of the LES runs, the horizontal grid resolution was x = y =4z LES 1 .Only very shallow PBL, where the flow is naturally anisotropic, was simulated with a grid resolution x = y =8z LES 1 .The difference in the vertical resolution causes some problems in data analysis.Obviously, the vertical resolution in the LES runs should not be coarser than the vertical resolution of the LSMs.Hence, the coarsest LES runs with small Z n cannot be involved in the evaluation of the finest LSMs.Simultaneously, the first computational level in the LSMs must be still placed within the PBL.Hence, the finest LES runs with large Z n cannot be involved in the evaluation of the coarsest LSMs.In order to avoid these problems, three specific LSM levels were selected.The drag coefficient for momentum can be calculated in the LES as where z LES  All LES runs accurately resolve the entire PBL.The PBL depth, H , can be diagnosed as a height, where τ (H )=0.05u 2 * .It is worth noticing that other definitions of H are also possible.It means that H is an ambiguous parameter.The definition of H affects values of empirical constants in the following analytical expressions.The value of the friction velocity, u * , is also defined only approximately.Indeed, advection and pressure terms have been neglected in the layer between the surface and z LES 1 .This is not justified for the coarse LES runs.This simplification and the noticeable overamplification of the resolved TKE at z LES 1 result in a slightly larger value of the von Karman constant in the LES.The value κ=0.44 is used in this study.

Methods to improved C Dn -parameterizations
Figure 5 shows C LES Dn normalized by the classical local scaling by Eq. ( 7).Recently, comparable deviations have been found in the atmospheric data (Lange et al., 2004).This data is also shown in Fig. 5. Apparently, the local scaling, z/z 0 , does not work properly in the shallow PBLs.There are several possible ways to incorporate non-local effects of the ambient stratification into the surface drag parameterization.One can modify the von Karman length scale, l, in such a way that it would account for the actual position of z LSM 1 within the PBL.Indeed, it is understandable that l must be limited.Both LES and data show saturation of the mixing length scale at the top of the surface layer (see Fig. 6).Obviously, the saturation of l has to be taken into account at z LSM 1 >h s .Blackadar (1962) proposed the most straightforward extension of the length scale.He suggested a linear interpolation between reciprocals to the von Karman, κz, and external, κH , scales.His external scale was, however, based on Rossby and Montgomery (1935) analytical expression for the depth of the idealized Ekman layer.More exactly, Blackadar (1962) used the following external scale λ=2.7•10 −4 U g /|f |.This scale gives λ≈0.01H , where H is given by Eq. ( 8).After Blackadar (1962) work, a number of modification of this scaling for stable and convective PBLs have been proposed (Holt and Raman, 1988).Mason and Thompson (1992) found that the best result in the neutral PBL is obtained by interpolation of squared reciprocals.The Blackadar's method immediately involves H .It can be further developed to an expression in terms of large-scale variables N ,f and U g .
There is also another method possible.It does not change l.It proposes a linear combination of an additional non-local term and the log-law (Zilitinkevich et al., 2002).The method will be referred to as the Zilitinkevich's method.

Empirical relation Figure 5 reveals that C LES
Dn systematically deviates from the value predicted by Eq. ( 7).One can assume that l/H is a logarithmic function of z/H .The author's attempts to find an empirical form for this function have resulted in an expression with at least three empirical constants.This expression gave an integral in a slowly converging series of (ln(H /λ 0 )) n , where n=0, −1, . ... This is apparently not a constructive approach.However, these attempts suggested a possible form of an empirical relation between C Dn and H , namely, Dn is given by Eq. ( 15).Symbols and lines are the same as in Fig. 5.
where C emp =0.5 is an empirical constant and λ 0 =1 m is a normalization length scale.This empirical relation is easy to apply in practical LSMs but it does not have physical argumentation behind.
4.2 Blackadar's method Blackadar (1962) incorporated the PBL depth into the loglaw by interpolating reciprocals where C 1 is a dimensionless constant, which accounts for the relative significance of the external scaling.The original Blackadar's estimation gave too large C 1 ≈50.It has been consequently reduced in Holt and Raman (1988), C 1 ≈33, and further in Mason and Thompson (1992), C 1 ≈7.In this work, C 1 is reduced even more.The best fit gives C 1 =3.
The same value was specified in Ballard, Colding and Smith (1991).Substitution of Eq. ( 14) into Eq.( 2) and integration gives Figure 7 shows C LES Dn normalized by the non-local scaling in Eq. ( 15).There is considerably smaller systematic bias of data in Fig. 7 than in Fig. 5. Equation ( 15) predicts the values of C LES Dn with an accuracy of ±20%.Equation ( 15) can be seen as a direct implementation of the Monin-Obukhov theory with a new universal function nl (z/H )= − C 1 z/H .Hence, Eq. ( 5) becomes Dn is given by Eq. ( 18).Symbols and lines are the same as in Fig. 5.
This parameterization works as well for the deep PBL capped with a weak inversion as it does for the shallow and very shallow PBLs capped with a strong inversion.Equation (15) recovers Eq. ( 7) only when H →∞. This is physically inconsistent as it should recover the classical C Dnparameterization at H given by Eq. ( 8).
Direct accounting for the equilibrium depth of the Ekman layer is problematic.However, it suggests that the Mason and Thompson (1992) length scale will probably work better.This mixing length scale reads The best fit for the empirical constant C 2 is 2.4.It is perhaps also possible to adopt C 1 =C 2 .Substitution of Eq. ( 17) into Eq.( 2) with integration gives .
Figure 8 shows C LES Dn normalized by the non-local scaling in Eq. ( 18).There is noticeable improvement in the C Dn prediction, especially for large H .The overall accuracy of the prediction is still about ±20%.
Unfortunately, one cannot use Eq. ( 15) and Eq. ( 18) in the LSMs directly.These expressions contain the unknown parameter H .To express H in terms of resolved variables and boundary conditions, one can involve Eq. ( 10).In this equation, u * is also unknown parameter.One can assume that u 2 * is the surface turbulent stress, which is induced by the mean is given by Eq. ( 19).Symbols and lines are the same as in Fig. 5.
wind in a deep PBL.Implicitly, it assumes that the effect of the large eddies is entirely accounted for through Z n .This is a strong assumption but it allows one to use the classical log-law from Eq. ( 7) for elimination of u * .This assumption also demands that one changes the values of constants C 0 , C 1 and C 2 .After the substitution, Eq. ( 15) and Eq. ( 18) read The empirical constants are Dn .The data scatter is noticeably reduced in Fig. 9.The overall accuracy of the C G1 Dn expression is about ±10%.This is within the accuracy of calculations of C LES Dn .Moreover, it is generally better than the accuracy of atmospheric measurements over complex surfaces.The C G2 Dn is given by Eq. ( 20).Symbols and lines are the same as in Fig. 5.
expression exhibits larger data scatter.Its overall accuracy is

±20%.
The major problem of the Blackadar's method is that it does not recover the classical C Dn -parameterization in the limit Z n →0 exactly.The inconsistency is small.Nevertheless, a more consistent method is still desirable.

Zilitinkevich's method
The above-mentioned inconsistency does not appear in Zilitinkevich's method.The method is based on the idea that the increasing temperature flux modifies the velocity profile above the surface layer.This modification is similar to the modification of the velocity profile in the Monin-Obukhov similarity theory.The only difference is that the temperature flux is not due to surface cooling but due to PBL top warming.Mathematically, one can write (Zilitinkevich et al. , 2002) where a u =0.35 is an empirical constant.Integration of Eq. ( 21) gives a modified expression for the turbulent friction velocity Substitution of Eq. ( 22) into Eq.( 7) results in is given by Eq. ( 23).Symbols and lines are the same as in Fig. 5.
This expression recovers the classical C Dparameterization in the limit Z n →0.
Figure 11 shows C LES Dn normalized by Eq. ( 23).Zilitinkevich's method is rather accurate; its accuracy is about ±10% within a typical range of the atmospheric H .
It is also possible to rewrite Eq. ( 23) in terms of the geostrophic wind speed, U g .The expression reads Figure 12 shows C LES Dn normalized by Eq. ( 25), an equation that Eq. ( 25) works slightly worse than Eq. ( 23); however, the accuracy is still within ±20%.

Significance of improvements
A number of field measurements (Tjernström and Smedman, 1993) reveals that the near-neutral atmospheric PBL is usually quite shallow.The typical range of the atmospheric H is from 200 m to 600 m.The near-neutral PBL can be even more shallow in the case of a slow geostrophic wind.It is worth mentioning that this study involves LES runs with U g ranging from 1 m s −1 to 15 m s −1 .The PBL depth was just about 40 m at U g =1 m s −1 and Z n =100.Figure 5 shows the deviation of the actual C LES Dn from the predicted C Dn in the classical parameterization.The deviation becomes significant at H <600 m for the coarse resolution LSMs and at H <200 for the fine resolution LSMs.This fact generally corroborates the common assumption that h s ≈0.1H .As a matter of fact, the surface layer in LES is deeper than the atmospheric surface layer.This is due to a smaller effective Reynolds number of LES.One can argue that the discrepancy between actual and predicted C Dn is even larger in the atmospheric PBL.
Figure 13 shows the ratio of the improved parameterizations of the surface drag coefficients, C G1 Dn , C G2 Dn , C E Dn and is given by Eq. ( 25).Symbols and lines are the same as in Fig. 5.
C Z Dn to the classical parameterization by Eq. ( 7).The ratio is close to unity in the deep PBLs (small Z n ) but it rapidly decreases in the shallow PBLs (large Z n ).The improved parameterizations predict significantly smaller turbulent exchange in the LSMs.For instance, the turbulent exchange can be reduced by 20%-30% under typical atmospheric conditions in the LSM.

Conclusions
Modern large-scale meteorological models parameterize the turbulent exchange between the first computational level in the atmosphere and the underlying surface using a concept of the surface drag coefficients.This concept assumes the existence of a layer of constant fluxes.In this layer, the mean velocity profile is logarithmic.Assuming the first computational level in the LSMs is placed within this layer, one can derive a very simple classical parameterization of the surface drag coefficients given by Eq. ( 7).
Direct atmospheric measurements and LES data revealed that the conventionally neutral PBL is usually much more shallow than the idealized Ekman layer.The equilibrium depth of the conventionally neutral PBL is strongly reduced by the ambient stratification of the free atmosphere above the PBL.This has been noticed by Csanady and explained in Zilitinkevich's theory of the non-local turbulence.The theory has also gained strong support from LES.The modern LSMs usually have a rather coarse vertical resolution.Hence, the LSMs may have the first computational level well above the log-layer in the realistic near-neutral PBL.The discrepancy between the parameterized and actual C Dn and C D can be significant.The classical parameterization systematically overestimates the turbulent exchange in the modern LSMs by 20%-30%.
It is possible to improve the parameterization of the surface drag coefficient.The improved parameterization must account for the non-local effect of the ambient atmospheric stratification.The effect can be quantified using the Zilitinkevich number, Z n .
To do this, two methods have been used.Blackadar's method is based on the redefinition of the mixing length scale, l, in the log-law.Zilitinkevich's method is based on the redefinition of the log-law itself.Both methods resulted in considerable improvements in the prediction of the surface drag coefficient, especially in the very shallow PBLs.However, the new parameterizations are not perfect.They still show irregular scatter of the actual values of C Dn around its predicted values.The best accuracy gives Eq. ( 19).In this case, the scatter is only ±10%.Further improvements require an LES of significantly better resolution in the surface layer.Equation ( 25) in Zilitinkevich's method gives comparably accurate parameterization.This equation is also more mathematically consistent.It recovers the classical C D -parameterization in both limits Z n →0 and z→0.

Fig. 1 .
Fig. 1.The dependence of the normalized depth, C h =H |f |/u * , of the conventionally neutral PBL on the stratification of the free atmosphere above the PBL.The strength of the stratification is expressed through the Zilitinkevich number, Z n .Data is taken from the author's LES database.Squares, circles and diamonds denote individual LES runs with the log-layer more shallow than 20 m, 50 m and 100 m.

Fig. 2 .
Fig. 2. The dependence of the turbulent surface stress, τ s =u 2 * on the PBL depth H . Data and notation are as in Fig. 1.

Fig. 3 .
Fig. 3.Values of the constant C u of the Monin-Obukhov similarity theory for the non-dimensional velocity gradient given by Eq. (11).The values of C u are plotted against non-dimensional parameter z/L, accounting both for the stability of the surface layer and the height of the computational level in LES.Squares represent LES data; other symbols represent atmospheric data.

Fig. 4 .
Fig. 4. The dependence of the surface drag coefficient for momentum, C LES Dn , on the internal dimensionless parameter z/z 0 .C LES Dn is calculated from the author's LES database at z LSM 1 =20 m (squares), z LSM 1 =60 m (circles) and z LSM 1 =150 m (diamonds).
C Dn at z LSM 1 =20 m represents the drag coefficient in the LSMs with a relatively fine vertical resolution.LES data at this level is denoted by squares in the following text.C Dn at z LSM 1 =60 m represents the drag coefficient in the LSMs with a typical vertical resolution.LES data at this level is denoted by circles.C Dn at z LSM 1 =150 m represents the drag coefficient in the LSMs with a relatively coarse vertical resolution.LES data at this level is denoted by diamonds.The analysis involves only those LES runs which satisfy z drag coefficient is shown in Fig.4as a function of the local parameter z LSM 1 /z 0 .All variables in C Dn computation are averaged over the horizontal plain and additionally over a 30-min interval of time.The level z LES 1 is always placed well within the log-layer.It ensures that the parameterization in Eq. (1) and Eq.(7) provides an accurate approximation of the turbulence exchange in the LES.The non-local effects, which are the subject of the study, are hidden in the well resolved ratio |u(z LES 1

Fig. 5 .
Fig. 5.The dependence of the normalized surface drag coefficient C LES Dn /C Dn on the PBL depth H . C Dn is the classical parameterization by Eq. (7).The solid line shows the ideal scaling.The dashed line shows an empirical dependence by Eq. (13).Crosses and errorbars represent atmospheric data for C atm D /C Charnok D from Lange et al. (2004), where H has been recalculated according to Eq. (10) in the present study from Eq. (17) in Lange et al. work.Other symbols are the same as in Fig. 4.

Fig. 6 .
Fig. 6.The dependence of the non-dimensional mixing length scale, l/H , from the non-dimensional height, z/H .The mixing length scale is calculated as l=τ (z)/∇ z |u(z)|u −1 * from the author's LES database.

Fig. 7 .
Fig. 7.The dependence of C LES Dn /C H 1 Dn on the PBL depth H . C H 1

Fig. 8 .
Fig. 8.The dependence of C LES Dn /C H 2 Dn on the PBL depth H . C H 2

Fig. 9 .
Fig. 9.The dependence of C LES Dn /C G1 Dn on the PBL depth H . C G1 Dn

Fig. 10 .
Fig. 10.The dependence of C LES Dn /C G2 Dn on the PBL depth H . C G2 Dn

Fig. 11 .
Fig. 11.The dependence of C LES Dn /C Z Dn on the PBL depth H . C Z Dn

Fig. 12 .
Fig. 12.The dependence of C LES Dn /C ZG Dn on the PBL depth H . C Z Dn