Universal Wind Profile for Conventionally Neutral Atmospheric Boundary Layers
aa r X i v : . [ phy s i c s . a o - ph ] F e b Accepted by Physical Review Letters
Universal Wind Profile for Conventionally Neutral Atmospheric Boundary Layers
Luoqin Liu ID , ∗ Srinidhi N. Gadde ID , and Richard J. A. M. Stevens ID † Physics of Fluids Group, Max Planck Center Twente for Complex Fluid Dynamics,University of Twente, 7500 AE Enschede, The Netherlands (Dated: February 11, 2021)Conventionally neutral atmospheric boundary layers (CNBLs), which are characterized with zerosurface potential temperature flux and capped by an inversion of potential temperature, are fre-quently encountered in nature. Therefore, predicting the wind speed profiles of CNBLs is relevantfor weather forecasting, climate modeling, and wind energy applications. However, previous at-tempts to predict the velocity profiles in CNBLs have had limited success due to the complicatedinterplay between buoyancy, shear, and Coriolis effects. Here, we utilize ideas from the classicalMonin-Obukhov similarity theory in combination with a local scaling hypothesis to derive an ana-lytic expression for the stability correction function ψ = − c ψ ( z/L ) / , where c ψ = 4 . z is the height above ground, and L is the local Obukhov length based on potential tem-perature flux at that height, for CNBLs. An analytic expression for this flux is also derived usingdimensional analysis and a perturbation method approach. We find that the derived profile agreesexcellently with the velocity profile in the entire boundary layer obtained from high-fidelity largeeddy simulations of typical CNBLs. Introduction . For well over a century, wall-boundedturbulent flows have been studied extensively [1]. A focusarea is the derivation and characterization of the meanvelocity profile. In 1925, Prandtl [2] recognized that thevelocity profile in the inertial sub-layer is approximatelylogarithmic based on his mixing length hypothesis. Inthe 1930s, von K´arm´an [3, 4] derived the logarithmic lawof the wall analytically using dimensional analysis. In1956, Coles [5] showed using measurement data that thevelocity profile can be described more accurately via thesum of the logarithmic law and a wake function. Sincethen, the law of the wall has been the pillar of the de-scription of wall-bounded turbulence [6–8]. Recently, theuniversality of the law of the wall has been supported bytheoretical and experimental studies [9–11].The dynamics in atmospheric boundary layers, wheremost human activity and biological processes occur, aremuch more complicated as turbulence is generated byshear stress and buoyancy [12], while the Coriolis forcecreates a wind veer [13]. In 1954, Monin and Obukhov[14] introduced a stability correction function ψ to ac-count for deviations to the logarithmic wind speed profilecaused by thermal stratification. Based on the Bucking-ham Π theorem, Monin and Obukhov concluded that ψ is only a function of the atmospheric stability parameter z/L w , where z is the vertical height above the groundand L w is the Obukhov length based on the surface po-tential temperature flux. The universality of this wellknown Monin-Obukhov similarity theory (MOST) [14–16] has been established in the surface layers of stable andconvective atmospheric boundary layers in many field ex-periments [17], as well as large eddy simulations (LES)[18]. Therefore, the MOST is nowadays regarded as thestarting point of modern micro-meteorology [17].Conventionally neutral atmospheric boundary layers(CNBLs) are also frequently observed and are often con-sidered in fundamental studies [19–26]. In contrast to stable and convective atmospheric boundary layers, CN-BLs are characterized with zero surface potential temper-ature flux and capped by an inversion of potential tem-perature. However, the classical MOST is not applica-ble to CNBLs, because the surface potential temperatureflux is zero, due to which L w is no longer a relevant scale[18]. Many studies [23, 27–35] have tried to predict thevelocity profile in CNBLs, but so far with limited successsince none of them have considered the effect of poten-tial temperature flux near the capping inversion layer.Therefore, the logarithmic law of the wall without stabil-ity correction is still commonly used to predict the windspeed profile in CNBLs [36–40]. However, an analyticaldescription of the velocity profiles in CNBLs is of greatfundamental interest and relevant for meteorological ap-plications [12, 35]. Theory . In this work, we derive the potential tempera-ture flux profile using dimensional analysis and a pertur-bation method approach. To account for the deviation tothe logarithmic wind speed profile we use ideas from theMOST in combination with a local similarity hypothesis.Therefore, we introduce a stability correction function ψ that depends only on the local stability parameter z/L ,where L is the local Obukhov length based on the localpotential temperature flux. The canonical shape of ψ isdetermined by asymptotic analysis. The derived univer-sality profiles for the potential temperature flux and thewind speed profiles are confirmed by the excellent agree-ments with the results of high-fidelity LES.Based on the dimensional analysis (a derivation can befound in section III of the Supplementary Material [41]),the potential temperature flux can be written as βz qu ∗ = − Π (
Ro, Zi, ξ ) = − Ro r Zi s Π ( ξ ) . (1)Here q is the potential temperature flux, u ∗ is the fric-tion velocity, β is the buoyancy parameter, Π and Π are ccepted by Physical Review Letters Ro = u ∗ / ( | f | z ) is the Rossbynumber with f the Coriolis parameter and z the rough-ness height, Zi = N/ | f | is the Zilitinkevich number [28]with N the free-atmosphere Brunt–V¨ais¨al¨a frequency,and r and s are the power exponents for Ro and Zi , re-spectively. The dimensionless parameter ξ = z/h ′ , where h ′ = h/ (1 − . / ) and h is the boundary layer heightat which the total momentum flux reaches 5% of the sur-face value. This conventional definition underestimatesthe actual boundary layer height. The functional form of h ′ follows from the fact that the dimensionless total mo-mentum flux follows a power law with exponent 3 / Ro and Zi . The values of s and r will be determined later fromour high-fidelity LES data.We emphasize that other definitions for the boundarylayer height, which are based on, for example, the verti-cal wind speed or potential temperature profiles are alsocommonly used [25, 31, 35]. In particular, the bound-ary layer height h t is defined as the height at which thepotential temperature flux reaches its minimum value.Previous studies [25, 35, 43–45] showed that the poten-tial temperature flux in CNBLs decreases linearly fromzero at the surface to a minimum value at z = h t , andthen increases to zero for z ≥ h ′ . As explained in sec-tion IV of the Supplementary Material [41], the ratio h t /h ′ ≡ − ǫ is a function of Zi . However, as we willsee later the dependence of ǫ on Zi is limited over theparameter regime under consideration. Clearly, ǫ ≪ h ′ , where the potential temperature flux recov-ers steeply to zero. Therefore, we propose the followingordinary differential equation to model the potential tem-perature flux, − ǫ Π ′′ + Π ′ = c Π , Π (0) = Π (1) = 0 . (2)Here c Π is the slope of the dimensionless total poten-tial temperature flux − ( βz q ) / ( u ∗ Zi s Ro r ) in the surfacelayers, which can be determined from simulation or mea-surement data. The solution of Π readsΠ = c Π (cid:18) ξ − e ξ/ǫ − e /ǫ − (cid:19) , ξ ≤ , , ξ > . (3)Note that Eq. (2) and its solution Eq. (3) is reminiscentof the classical singular perturbation method [46]: Theouter solution (close to the wall) is a linear function of ξ and the inner solution (close to the capping inversionlayer) is controlled by a small parameter ǫ .In contrast to the classical MOST [14] where the nor-malized wind speed gradient is assumed to be a universal function, we introduce a stability correction function ψ to account for the deviation of the logarithmic profile.Therefore we write the wind speed profile as κU mag u ∗ = ln (cid:18) zz (cid:19) − ψ (cid:16) zL (cid:17) , (4)where κ = 0 . U mag is themean wind speed, and ψ is the stability correction func-tion that depends only on the dimensionless stability pa-rameter z/L . According to the local scaling hypothesis[42, 47], L is defined as the local Obukhov length,1 L ≡ − κβqu ∗ . (5)It is worth to point out that the dimensionless slope( κz/u ∗ )d U mag / d z is usually regarded as a universal func-tion of the stability parameter z/L in the stable and con-vective atmospheric boundary layers [48]. However, un-der the assumption of Eq. (4), this slope is no longer auniversal function of z/L .To determine the canonical shape of ψ , we assume ψ = − c ψ (cid:16) zL (cid:17) p , (6)where p is the power exponent to be determined analyti-cally below, and c ψ is an empirical constant. Recall thatvery close to the wall (see Eq. (3))Π → c Π ξ = c Π zh ′ as zh ′ → . (7)Then, from asymptotic analysis [49], we find that zL = κzz Π → c Π κh ′ z Zi s Ro r (cid:16) zh ′ (cid:17) as zh ′ → . (8)Zilitinkevich and Esau [23] showed that in the surfacelayers of stable, truly neutral, and conventionally neutralatmospheric boundary layers ψ = − C u z/L M . Here C u is a dimensionless constant and L M is the combined tur-bulent length scale, which in CNBLs can be estimated as | f | L M /u ∗ = (1 + C m Zi ) − / , where C m is an empiricalconstant ([23], see also the Supplementary Material [41]).To match with the result of Zilitinkevich and Esau [23]in the surface layer, we find that p = 1 /
2. Clearly, thedetermination of p is independent of the values of c Π , c ψ , r , s , and ǫ . Thus, the wind speed profile is given by κU mag u ∗ = ln (cid:18) zz (cid:19) + c ψ (cid:16) zL (cid:17) / , ξ ≤ ξ ,κGu ∗ , ξ > ξ . (9)Here G is the geostrophic wind speed, ξ is the highest in-tersection point of the curves described by the upper andlower expressions in Eq. (9), z/L is the dimensionless sta-bility parameter predicted by the potential temperature ccepted by Physical Review Letters TABLE I. Summary of all simulated cases, where the Zi and Ro range covers the values found in typical CNBLs at mid tohigh latitudes [58, 59]Case no. A B C D E F Zi Ro ǫ c Π τ /τ w . Filled symbols: LES data; solid line: theoret-ical curve given by τ /τ w = (1 − z/h ′ ) / . flux model (i.e. Eqs. (1) and (3)), and c ψ is the empiri-cal constant that can be determined from simulation ormeasurement data. Validation . To verify the universality of the wind speedprofile for CNBLs, we perform six high-fidelity LES. Inthe simulations, a CNBL over a flat surface with periodicconditions in horizontal directions is considered. Theflow is initialized with uniform geostrophic wind speedand a linear potential temperature profile with a constantgradient [35, 44]. The simulations are performed with anin-house code [50–55], which employs a pseudo-spectraldiscretization in the horizontal directions and a second-order finite difference method in the vertical direction.We employ the advanced anisotropic minimum dissipa-tion model to parameterize the sub-grid scale shear stressand potential temperature flux [56]. The horizontal do-main size is more than six times larger than the boundarylayer height, and the grid resolution is 288 . We ensurethat all simulations have reached the quasi-stationarystate and the statistics are averaged over one inertial pe-riod [57]. A summary of all simulated cases is presentedin Table I. The simulated Zilitinkevich number Zi andRossby number Ro range covers the values found in typi-cal CNBLs at mid to high latitudes [58, 59]. More detailsabout the numerical method and simulation setup can befound in the Supplementary Material [41].Figure 1 shows the vertical profile of the dimensionlessmean total momentum flux τ /τ w , where τ is the total ( a )( b )FIG. 2. Dimensionless mean potential temperature gradientin the surface layer versus (a) the Rossby number Ro and (b)the Zilitinkevich number Zi , where q ′ = d q/ d ξ . The valuesof the slope are (a) r = − . ≈ − s = 1 . ≈ momentum flux and τ w is its surface value. All differentcases in Table I are shown in the figure (filled symbols).Nieuwstadt [42] analytically determined that the totalmomentum flux profile in stable atmospheric boundarylayers scales as τ /τ w = (1 − z/h ) / . In Fig. 1 we showthat this expression is still valid for CNBLs when we con-sider the previously introduced boundary layer thickness h ′ = h/ (1 − . / ). The finding that the dimensionlessmomentum flux profiles obtained from all LES collapseto the theoretical curve (see Fig. 1) confirms that h ′ isthe appropriate boundary layer height scale to consider.To determine the values of the power indices r and s , we take the vertical derivative of Eq. (1), see detailsin section III of the Supplementary Material [41]. Fig-ure 2 shows the dimensionless mean potential tempera-ture gradient ln ( − βz q ′ /u ∗ ) versus (a) the Rossby num-ber ln Ro and (b) the Zilitinkevich number ln Zi in thesurface layer, where q ′ = d q/ d ξ . The slopes of the curveshown in the figure determine the values of the powerexponents r and s . In the parameter regime under con-sideration r = − . ≈ − s = 1 . ≈ c Π canalso be determined from the figure and the results arelisted in Table I.Figure 3 shows the vertical profile of the dimensionlessmean potential temperature flux q/ | q | max , which reducesto almost zero at z/h ′ ≥
1. The potential temperature ccepted by Physical Review Letters FIG. 3. Vertical profile of dimensionless mean potential tem-perature flux q/ | q | max . Filled circles: LES data of the presentwork; filled triangles: prediction based on a turbulence clo-sure given by Mauritsen et al. [43]; filled inverted triangles:DNS data of Jonker et al. [60]; filled diamonds: LES data ofPedersen et al. [44]; filled squares: LES data of Allaerts andMeyers [25]; filled stars: LES data of Berg et al. [45]; solidline: theoretical prediction given by Eq. (3) with ǫ = 0 . flux first decreases linearly from zero at the surface to aminimum value at z = h t ≡ (1 − ǫ ) h ′ , and then increasesrapidly to zero in a narrow region (1 − ǫ ≤ z/h ′ ≤
1) since ǫ ≪
1. The value of ǫ is expected to dependonly on Zi (see the Supplementary Material [41]). Thedata in Table I show that for the parameter range underconsideration the variation of ǫ is limited and thereforewe take ǫ = 0 .
12 to describe the data here. Evidently,all LES data of the present work (filled symbols) collapsevery well to the introduced theoretical model (solid line),which validates the chosen approach. For comparison,the prediction based on a turbulence closure given byMauritsen et al. [43], the direct numerical simulations(DNS) data performed by Jonker et al. [60], and the LESdata taken from Pedersen et al. [44], Allaerts and Meyers[25], and Berg et al. [45] are also shown in the figure. Theoverall agreement between the theoretical prediction andthe data from previous studies [25, 43–45, 60] is verygood, which confirms the universality of the proposedpotential temperature flux profile.Figure 4 shows the vertical profile of the dimensionlesswind speed for two typical cases, which covers the Zi and Ro number range of typical CNBLs at mid to highlatitudes [58, 59]. The filled symbols are LES data, thedashed line is the theoretical prediction given by the loga-rithmic law, the blue line is the prediction of Zilitinkevichand Esau [23], the yellow line is the prediction of Gryn-ing et al. [29], the red line is the prediction of Kelly et al. [35], and the black line is the prediction given by Eq. (9)with c ψ = 4 . c Π = 0 . ǫ = 0 . ( a )( b )FIG. 4. Vertical profile of mean wind speed for (a) case B and(b) case F. Filled symbols: LES data; dashed line: predictiongiven by the logarithmic law; blue line: prediction given byZilitinkevich and Esau [23]; yellow line: prediction given byGryning et al. [29]; red line: prediction given by Kelly et al. [35]; black line: prediction given by Eq. (9) with c ψ = 4 . (see Table I). The empirical constant c ψ is determinedsuch that it can predict the wind speed profiles of allcases in Table I with minimum discrepancies. The figureshows that the logarithmic law only accurately capturesthe wind speed in the lower 10% of the boundary layer,also known as the surface layer (shaded region). The the-ory given by Gryning et al. [29] focuses on capturing thewind speed at the top of the CNBL, but does not cap-ture the effect of the low-level jet. The theory given byKelly et al. [35] is focused on the lower part of the CNBL.The predictions by Zilitinkevich and Esau [23] agree wellwith the LES data in the lower part of the CNBL, butdo not capture the low-level jet, which is represented inour approach. In contrast, the agreement between theproposed profile (9) and the LES data is nearly perfectin the entire boundary layer and much better than allprevious approaches. This excellent agreement confirmsthe universality of our proposed wind profile (9) in theconsidered parameter range of CNBLs. Summary . We propose a universal velocity profile forCNBL derived using a local similarity hypothesis com-bined with ideas from the classical Monin-Obukhov sim- ccepted by Physical Review Letters ψ to account for the deviation of the logarithmiclaw. The canonical shape of ψ is determined theoret-ically as ψ = − c ψ ( z/L ) / , where c ψ = 4 . z isthe vertical height above the surface, and L is the lo-cal Obukhov length. An analytical expression for thepotential temperature flux profile is also derived fromdimensional analysis and perturbation method. The uni-versality of the proposed profile (9) has been confirmedby its excellent agreement with high-fidelity LES resultsfor Ro = [4 . × , . × ] and Zi ∈ [51 , Zi and Ro number range cover the range of valuesobserved in typical CNBLs at mid to high latitudes. Fur-ther work is required to assess the applicability of theapproach to other parameter regimes.We appreciate very much the valuable comments of theanonymous referees. We acknowledge Drs. K. L. Chongand Y. X. Li for insightful discussion. This work is part ofthe Shell-NWO/FOM-initiative Computational sciencesfor energy research of Shell and Chemical Sciences, Earthand Live Sciences, Physical Sciences, FOM and STW,and an STW VIDI grant (No. 14868). This work wascarried out on the national e-infrastructure of SURFsara,a subsidiary of SURF cooperation, the collaborative ICTorganization for Dutch education and research. ∗ [email protected] † [email protected][1] A. J. Smits, B. J. McKeon, and I. Marusic, High Reynoldsnumber wall turbulence, Annu. Rev. Fluid Mech. , 353(2011).[2] L. Prandtl, Bericht ¨uber Untersuchungen zur ausgebilde-ten Turbulenz, Z. Angew. Math. Mech. , 136 (1925).[3] T. von K´arm´an, Mechanische ¨Ahnlichkeit und Turbulenz,Nachr. Ges. Wiss. G¨ottingen Math. Phys. Klasse , 58(1930).[4] T. von K´arm´an, Mechanische ¨Ahnlichkeit und Turbulenz,Proc. Third Internat. Congr. Appl. Mech. Stockholm ,85 (1931).[5] D. Coles, The law of the wake in the turbulent boundarylayer, J. Fluid Mech. , 191 (1956).[6] H. Tennekes and J. L. Lumley, A First Course in Turbu-lence (The MIT Press, Cambridge, Massachusetts, 1972).[7] S. B. Pope,
Turbulent Flows (Cambridge UniversityPress, Cambridge, 2000).[8] P. Davidson,
Turbulence: An Introduction for Scientistsand Engineers (Oxford University Press, Oxford, 2004).[9] I. Marusic, J. P. Monty, M. Hultmark, and A. J. Smits,On the logarithmic region in wall turbulence, J. FluidMech. , R3 (2013).[10] P. Luchini, Universality of the turbulent velocity profile,Phys. Rev. Lett. , 224501 (2017).[11] M. Samie, I. Marusic, N. Hutchins, M. K. Fu, Y. Fan,M. Hultmark, and A. J. Smits, Fully resolved measure-ments of turbulent boundary layer flows up to Re τ = 20000, J. Fluid Mech. , 391 (2018).[12] G. G. Katul, A. G. Konings, and A. Porporato, Meanvelocity profile in a sheared and thermally stratified at-mospheric boundary layer, Phys. Rev. Lett. , 268502(2011).[13] M. F. Howland, A. S. Ghate, and S. K. Lele, Influenceof the geostrophic wind direction on the atmosphericboundary layer flow, J. Fluid Mech. , A39 (2020).[14] A. S. Monin and A. M. Obukhov, Basic laws of turbulentmixing in the surface layer of the atmosphere, Tr. Akad.Nauk SSSR Geophiz. Inst. , 163 (1954).[15] A. M. Obukhov, Turbulence in an atmosphere with inho-mogeneous temperature, Trans. Inst. Teoret. Geoz. Akad.Nauk SSSR , 95 (1946).[16] A. S. Monin and A. M. Yaglom, Statistical Fluid Mechan-ics, Vol. 1. Mechanics of Turbulence (The MIT Press,Cambridge, Massachusetts, 1971).[17] T. Foken, 50 years of the Monin–Obukhov similarity the-ory, Boundary-Layer Meteorol. , 431 (2006).[18] S. Khanna and J. G. Brasseur, Analysis of Monin-Obukhov similarity from large-eddy simulation, J. FluidMech. , 251 (1997).[19] R. A. Brost, D. H. Lenschow, and J. C. Wyngaard, Ma-rine stratocumulus layers. Part 1: Mean conditions, J.Atmos. Sci. , 800 (1982).[20] A. L. M. Grant, Observations of boundary layer structuremade during the 1981 KONTUR experiment, Q. J. R.Meteorol. Soc. , 825 (1986).[21] M. Tjernstr¨om and A.-S. Smedman, The vertical turbu-lence structure of the coastal marine atmospheric bound-ary layer, J. Geophys. Res.: Oceans , 4809 (1993).[22] S. S. Zilitinkevich and I. N. Esau, On integral mea-sures of the neutral barotropic planetary boundary layer,Boundary-Layer Meteorol. , 371 (2002).[23] S. S. Zilitinkevich and I. N. Esau, Resistance and heat-transfer laws for stable and neutral planetary boundarylayers: Old theory advanced and re-evaluated, Q. J. R.Meteorol. Soc. , 1863 (2005).[24] S. S. Zilitinkevich, I. Esau, and A. Baklanov, Furthercomments on the equilibrium height of neutral and stableplanetary boundary layers, Q. J. R. Meteorol. Soc. ,265 (2007).[25] D. Allaerts and J. Meyers, Large eddy simulation ofa large wind-turbine array in a conventionally neutralatmospheric boundary layer, Phys. Fluids , 065108(2015).[26] D. Allaerts and J. Meyers, Boundary-layer developmentand gravity waves in conventionally neutral wind farms,J. Fluid Mech. , 95 (2017).[27] S. S. Zilitinkevich, V. L. Perov, and J. C. King, Near-surface turbulent fluxes in stable stratification: Calcula-tion for use in general circulation models, Q. J. R. Mete-orol. Soc. , 1571 (2002).[28] I. N. Esau, Parameterization of a surface drag coefficientin conventionally neutral planetary boundary layer, Ann.Geophys. , 3353 (2004).[29] S.-E. Gryning, E. Batchvarova, B. Br¨ummer, H. Jor-gensen, and S. Larsen, On the extension of the wind pro-file over homogeneous terrain beyond the surface bound-ary layer, Boundary-Layer Meteorol. , 251 (2007).[30] M. Kelly and S.-E. Gryning, Long-term mean wind pro-files based on similarity theory, Boundary-Layer Meteo-rol. , 377 (2010).[31] M. Abkar and F. Port´e-Agel, The effect of free- ccepted by Physical Review Letters atmosphere stratification on boundary-layer flow andpower output from very large wind farms, Energies ,2338 (2013).[32] M. Optis, A. Monahan, and F. C. Bosveld, Movingbeyond Monin–Obukhov similarity theory in modellingwind-speed profiles in the lower atmospheric boundarylayer under stable stratification, Boundary-Layer Meteo-rol. , 497 (2014).[33] M. Kelly and I. Troen, Probabilistic stability and ‘tall’wind profiles: theory and method for use in wind resourceassessment, Wind Energy , 227 (2016).[34] Q. F. Jiang, S. P. Wang, and P. Sullivan, Large-eddy sim-ulation study of log laws in a neutral Ekman boundarylayer, J. Atmos. Sci. , 1873 (2018).[35] M. Kelly, R. A. Cersosimo, and J. Berg, A universalwind profile for the inversion-capped neutral atmosphericboundary layer, Q. J. R. Meteorol. Soc. , 982 (2019).[36] C. G. Rossby and R. B. Montgomery, The layers of fric-tional influence in wind and ocean currents, Pap. Phys.Oceanogr. Meteor. , 1 (1935).[37] A. K. Blackadar, The vertical distribution of wind andturbulent exchange in a neutral atmosphere, J. Geophys.Res. , 3095 (1962).[38] A. K. Blackadar and H. Tennekes, Asymptotic similar-ity in neutral barotropic planetary boundary layers, J.Atmos. Sci. , 1015 (1968).[39] H. Tennekes, The logarithmic wind profile, J. Atmos. Sci. , 234 (1973).[40] M. Horiguchi, T. Hayashi, A. Adachi, and S. Onogi,Large-scale turbulence structures and their contributionsto the momentum flux and turbulence in the near-neutralatmospheric boundary layer observed from a 213-m tallmeteorological tower, Boundary-Layer Meteorol. ,179 (2012).[41] See the Supplemental Material at [url], which includesadditional information on numerical simulations, deriva-tion of the potential temperature flux, formulations ofexisting wind speed models, and Refs. [1,2,5,7,8], .[42] F. T. M. Nieuwstadt, The turbulent structure of the sta-ble, nocturnal boundary layer, J. Atmos. Sci. , 2202(1984).[43] T. Mauritsen, G. Svensson, S. S. Zilitinkevich, I. N. Esau,L. Enger, and B. Grisogono, A total turbulent energy clo-sure model for neutrally and stably stratified atmosphericboundary layers, J. Atmos. Sci. , 4113 (2007).[44] J. G. Pedersen, S.-E. Gryning, and M. Kelly, On thestructure and adjustment of inversion-capped neutral at-mospheric boundary-layer flows: Large-eddy simulationstudy, Boundary-Layer Meteorol. , 43 (2014).[45] J. Berg, E. G. Patton, and P. P. Sullivan, Large-eddysimulation of conditionally neutral boundary layers: A mesh resolution sensitivity study, J. Atmos. Sci. , 1969(2020).[46] M. V. Dyke, Perturbation Methods in Fluid Mechanics (Parabolic Press, Stanford, California, 1975).[47] Z. Sorbjan, On similarity in the atmospheric boundarylayer, Boundary-Layer Meteorol. , 377 (1986).[48] J. A. Businger, J. C. Wyngaard, Y. Izumi, and E. F.Bradley, Flux-profile relationships in the atmosphericsurface layer, J. Atmos. Sci. , 181 (1971).[49] C. C. Lin and L. A. Segel, Mathematics Applied to Deter-ministic Problems in the Natural Sciences (Society for In-dustrial and Applied Mathematics, Philadelphia, 1988).[50] E. Bou-Zeid, C. Meneveau, and M. B. Parlange, A scale-dependent Lagrangian dynamic model for large eddysimulation of complex turbulent flows, Phys. Fluids ,025105 (2005).[51] R. J. A. M. Stevens, M. Wilczek, and C. Meneveau, Largeeddy simulation study of the logarithmic law for high-order moments in turbulent boundary layers, J. FluidMech. , 888 (2014).[52] S. N. Gadde and R. J. A. M. Stevens, Effect of Coriolisforce on a wind farm wake, J. Phys. Conf. Ser. ,012026 (2019).[53] S. N. Gadde, A. Stieren, and R. J. A. M. Stevens, Largeeddy simulations of stratified atmospheric boundary lay-ers: comparison of different sub-grid models, Boundary-Layer Meteorol. (2020).[54] L. Liu and R. J. A. M. Stevens, Effects of two-dimensionalsteep hills on the performance of wind turbines and windfarms, Boundary-Layer Meteorol. , 61 (2020).[55] L. Liu, S. N. Gadde, and R. J. A. M. Stevens, Geostrophicdrag law for conventionally neutral atmospheric bound-ary layers revisited, Q. J. R. Meteorol. Soc. , 1 (2020).[56] M. Abkar and P. Moin, Large eddy simulation of ther-mally stratified atmospheric boundary layer flow using aminimum dissipation model, Boundary-Layer Meteorol. , 405 (2017).[57] G. N. Coleman, J. H. Ferziger, and P. R. Spalart, Directsimulation of the stably stratified turbulent Ekman layer,J. Fluid Mech. , 677 (1992).[58] G. D. Hess and J. R. Garratt, Evaluating models of theneutral, barotropic planetary boundary layer using inte-gral measures: Part I. Overview, Boundary-Layer Mete-orol. , 333 (2002).[59] S. S. Zilitinkevich, S. A. Tyuryakov, Y. I. Troitskaya, andE. A. Mareev, Theoretical models of the height of the at-mospheric boundary layer and turbulent entrainment atits upper boundary, Izvestiya, Atmospheric and OceanicPhysics , 133 (2012).[60] H. J. J. Jonker, M. van Reeuwijk, P. P. Sullivan, andE. G. Patton, On the scaling of shear-driven entrainment:a DNS study, J. Fluid Mech. , 150 (2013). r X i v : . [ phy s i c s . a o - ph ] F e b Accepted by Physical Review Letters
Supplementary Material
Luoqin Liu ID , ∗ Srinidhi N. Gadde ID , and Richard J. A. M. Stevens ID † Physics of Fluids Group, Max Planck Center Twente for Complex Fluid Dynamics,University of Twente, 7500 AE Enschede, The Netherlands (Dated: February 11, 2021)
In sections I and II, we provide additional information on the used numerical method and the simulations’ setupand validation, respectively. Sections III and IV provide a derivation of Eq. (1) of the manuscript and ǫ dependenceon the Zilitinkevich number, respectively. A summary of previous models for CNBLs is given in section V. I. NUMERICAL METHOD
We perform large eddy simulations (LES) of conventionally neutral atmospheric boundary layers (CNBLs) over aflat surface with homogeneous roughness. We integrate the spatially-filtered Navier-Stokes equations and the filteredtransport equation for the potential temperature [1–3]: ∂ t e u + e ω × e u = f e z × ( G − e u ) + β ( e θ − h e θ i ) e z − ∇ e p − ∇ · τ , ∇ · e u = 0 , (1) ∂ t e θ + e u · ∇ e θ = −∇ · e q . (2)Here, the tilde denotes spatial filtering at a scale of ∆ = (∆ x ∆ y ∆ z ) / , h·i represents horizontal averaging, e u isthe velocity, e ω = ∇ × e u is the vorticity, e p is the modified pressure, e θ is the potential temperature, f is the Coriolisparameter, β is the buoyancy parameter, G is the geostrophic wind velocity, τ denotes the deviatoric part of thesub-grid scale (SGS) shear stress, and q represents the SGS potential temperature flux. Molecular viscosity anddiffusivity are neglected as the Reynolds number in the CNBL flow is very high. The SGS shear stress and potentialtemperature flux are parameterized using the recently developed anisotropic minimum dissipation model [4].Our code is an updated version of the one used by Albertson and Parlange [2]. The grid points are uniformlydistributed, 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 andpotential temperature grid planes are located at half a grid distance away from the ground. We use a second-order finite difference method in the vertical direction, while we use a pseudo-spectral discretization with periodicboundary 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 gradient, zero vertical velocity, and zero shearstress boundary condition. At the bottom boundary, we enforce a zero potential temperature flux and the classicalwall stress formulation [5, 6], τ xz = − (cid:18) κ ln z /z (cid:19) (cid:16)e u + e v (cid:17) / e u, τ yz = − (cid:18) κ ln z /z (cid:19) (cid:16)e u + e v (cid:17) / e v, (3)where the overline denotes filtering at a scale of 2∆, z is the half vertical grid distance, and τ xz and τ yz are the shearstresses in the streamwise and spanwise directions, respectively. The 2∆ filtering was suggested by Bou-Zeid et al. [6]to reduce the log-layer mismatch [7]. The reasoning for this procedure is that the filtering preserves the importantlarge scale variations. At the same time, the use of filtered velocities results in an average stress that is very close tothe stress predicted by the average similarity formulation for homogeneous surfaces [6]. The vertical derivative at thefirst horizontal plane is calculated using the Monin-Obukhov similarity theory without stability correction [5], whichis implemented following the method proposed by Albertson [1]. II. COMPUTATIONAL SETUP AND SIMULATION RESULTS
The computational domain size is 2 π km × π km × ∗ [email protected] † [email protected] ccepted by Physical Review Letters ( a ) ( b )FIG. 1. Vertical profiles of the mean dimensionless wind speed for case E with (a) different grid resolution on a 2 π km × π km × x, ∆ y, ∆ z ).( a ) ( b )FIG. 2. Vertical profiles of the mean (a) dimensionless wind speed and (b) potential temperature for cases B-E. streamwise flow structures are captured appropriately [8]. The vertical domain size is more than two times higherthan the boundary layer height. The simulations are performed on a 288 grid. The initial potential temperatureprofile is θ ( z ) = θ + Γ z , where θ = 300 K is the reference potential temperature and Γ is the free-atmospherelapse rate. The initial velocity profile is u = G e x , where G = 12 m/s is the geostrophic wind speed. Small randomperturbations are added to the initial fields of u and θ near the surface ( z ≤
100 m) to spin-up turbulence. Statisticsare collected over one inertial period, namely ∆( f t ) = 2 π , after the boundary layer has reached a quasi-stationarystate [9]. A summary of these simulations is given in Table I of the manuscript.To confirm that the main findings do not depend on the employed grid resolution and the domain size, we performedsimulations for case E (see Table I of the manuscript) using different grid resolutions (72 , 144 , and 288 ) and fordifferent domain sizes, i.e. 2 π km × π km × × ×
288 grid and π km × π km × × × III. DERIVATION OF THE POTENTIAL TEMPERATURE FLUX
Below we give a detailed explanation of Eq. (1) of the manuscript and the method we used to determine the valuesof s and r .For the boundary layers in quasi-equilibrium state, the mean profile is only a function of z , which is controlled by ccepted by Physical Review Letters q = q ( z ; f, G, β, Γ , z ) . (4)Therefore the general form of the dimensionless potential temperature can be written as qβz G = q (cid:18) zz ; Ro G = G | f | z , Zi = 1 | f | p β Γ (cid:19) . (5)Here and below q denotes the functional form of q . Using the geostrophic draw law for CNBLs [10, 11], (cid:18) κ Gu ∗ (cid:19) = h ln (cid:16) Ro G u ∗ G (cid:17) − A ( Zi ) i + B ( Zi ) , (6)where A = A ( Zi ) and B = B ( Zi ) are dimensionless coefficients, Eq. (5) can be written in terms of u ∗ as qβz u ∗ = q (cid:18) zz ; Ro = u ∗ | f | z , Zi (cid:19) . (7)In addition, the boundary-layer height for CNBLs can be estimated as [12] (cid:16) Ro z h ′ (cid:17) = 1 C R + 1 C N Zi, (8)where C R and C N are empirical constants. Thus, Eq. (7) can be written in terms of h ′ as follows qβz u ∗ = q (cid:16) ξ = zh ′ ; Ro, Zi (cid:17) = − Π( ξ ; Ro, Zi ) . (9)We assume that there is a power law dependence between Π and the two control parameters Ro and Zi , qβz u ∗ = − Π( ξ ; Ro, Zi ) = − Ro r Zi s Π ( ξ ) . (10)To determine the values of the power indices r and s in Eq. (10), we take the derivative of Eq. (10), which givesln (cid:18) − q ′ βz u ∗ (cid:19) = ln c Π + r ln Ro + s ln Zi, q ′ = d q d ξ . (11)This indicates that r and s can be determined from the slope of the curve ln( − q ′ βz /u ∗ ) vs ln Ro and ln( − q ′ βz /u ∗ )vs ln Zi , respectively, in the surface layer. Using a least-squares fitting, we find that (see figure 2 of the manuscript) r = − . ≈ − , s = 1 . ≈ . (12)Note that the value of c Π can also be determined from Eq. (11) once r and s have been determined. The correspondingresults are shown in Table I of the manuscript. IV. THE DEPENDENCE OF SMALL PARAMETER ON THE CONTROL PARAMETER
We remark that ǫ in Eq. (2) of the manuscript depends only on the control parameter Zi . To prove it, we recallthat the boundary-layer height for CNBLs can also be defined as the height at which the potential temperature fluxreaches its minimum value. We denote this height as h t , which can be estimated as [12, 13] (cid:18) Ro z h t (cid:19) = 1 C ′ R + 1 C ′ N Zi, (13)where C ′ R and C ′ N are empirical constants. Combining equations (8) and (13) gives2 ǫ ≡ − h t h ′ = 1 − C R + C N Zi C ′ R + C ′ N Zi ! / . (14)Note that the variance of ǫ with Zi is very small (see Table I of the manuscript) and therefore we take it as a constantin this work. ccepted by Physical Review Letters V. ADOPTED FORMULATIONS OF THE EXISTING WIND SPEED MODELS
Here we describe the previous models for the wind speed profiles in CNBLs, i.e. Zilitinkevich and Esau [10], Gryning et al. [14], and Kelly et al. [15], which are shown in figure 4 of the manuscript.
A. Zilitinkevich and Esau (Z2005 in figure 4 of the manuscript)
Zilitinkevich and Esau [10] showed that the wind shear in the surface layers can be written asd U mag d z = √ τκ (cid:18) z + C u L M (cid:19) , (15)where C u = 2 . L M is the combined turbulent length scale. In CNBLs, L f L M = 1 + ( C m Zi ) , L f = u ∗ | f | , (16)where C m = 0 . et al. [14], √ τ is approximated as √ τ = u ∗ (cid:16) − zh (cid:17) , (17)where h is the boundary layer height at which the total momentum flux reaches 5% of its surface value. By substitutingEq. (17) into Eq. (15), and after vertical integration we obtain κU mag u ∗ = ln (cid:18) zz (cid:19) + (cid:18) C u L M − h (cid:19) z − C u L M h z . (18) B. Gryning et al. (G2007 in figure 4 of the manuscript)
Gryning et al. [14] showed that the wind shear in the surface layer can be written asd U mag d z = √ τκ l , (19)where l is the local length scale that can be modelled by a combination of the length scale in the surface layer L SL ,in the middle of the boundary layer L MBL , and in the upper part of the boundary layer L UBL , i.e.1 l = 1 L SL + 1 L MBL + 1 L UBL . (20)In CNBLs, L SL = z, L UBL = h − z. (21)By substituting Eqs. (17), (20) and (21) into Eq. (19) and after vertical integration we obtain κU mag u ∗ = ln (cid:18) zz (cid:19) + zL MBL (cid:16) − z h (cid:17) . (22)The length scale L MBL is parameterized, for details see Gryning et al. [14] and Pedersen et al. [16], such that the windspeed at the top of the boundary layer conforms to the geostrophic wind speed.
C. Kelly et al. (K2019 in figure 4 of the manuscript)
Kelly et al. [15] proposed the wind speed profile as follows κU mag u ∗ = ln (cid:18) zz (cid:19) + a (cid:18) zL TD (cid:19) , (23) ccepted by Physical Review Letters a = 4 . L = L f Π θ Zi Ro . h θ , L f = u ∗ | f | . (24)Here Π θ = − . Ro h θ = L f /h θ is the Rossby number based on the boundary-layer height h θ , which is definedby the height at which the gradient of the potential temperature reaches its maximum value. [1] J. D. Albertson, Large Eddy Simulation of Land-Atmosphere Interaction , Ph.D. thesis, University of California (1996).[2] J. D. Albertson and M. B. Parlange, Surface length-scales and shear stress: implications for land-atmosphere interactionover complex terrain, Water Resour. Res. , 2121 (1999).[3] S. N. Gadde, A. Stieren, and R. J. A. M. Stevens, Large eddy simulations of stratified atmospheric boundary layers:comparison of different sub-grid models, Boundary-Layer Meteorol. (2020).[4] M. Abkar and P. Moin, Large eddy simulation of thermally stratified atmospheric boundary layer flow using a minimumdissipation model, Boundary-Layer Meteorol. , 405 (2017).[5] C.-H. Moeng, A large-eddy simulation model for the study of planetary boundary-layer turbulence, J. Atmos. Sci. , 2052(1984).[6] E. Bou-Zeid, C. Meneveau, and M. B. Parlange, A scale-dependent Lagrangian dynamic model for large eddy simulationof complex turbulent flows, Phys. Fluids , 025105 (2005).[7] J. G. Brasseur and T. Wei, Designing large-eddy simulation of the turbulent boundary layer to capture law-of-the-wallscaling, Phys. Fluids , 021303 (2010).[8] R. Foster, Signature of large aspect ratio roll vortices in synthetic aperture radar images of tropical cyclones, Oceanography , 58–67 (2013).[9] G. N. Coleman, J. H. Ferziger, and P. R. Spalart, Direct simulation of the stably stratified turbulent Ekman layer, J. FluidMech. , 677 (1992).[10] S. S. Zilitinkevich and I. N. Esau, Resistance and heat-transfer laws for stable and neutral planetary boundary layers: Oldtheory advanced and re-evaluated, Q. J. R. Meteorol. Soc. , 1863 (2005).[11] L. Liu, S. N. Gadde, and R. J. A. M. Stevens, Geostrophic drag law for conventionally neutral atmospheric boundarylayers revisited, Q. J. R. Meteorol. Soc. , 1 (2020).[12] S. S. Zilitinkevich, I. Esau, and A. Baklanov, Further comments on the equilibrium height of neutral and stable planetaryboundary layers, Q. J. R. Meteorol. Soc. , 265 (2007).[13] M. Abkar and F. Port´e-Agel, The effect of free-atmosphere stratification on boundary-layer flow and power output fromvery large wind farms, Energies , 2338 (2013).[14] S.-E. Gryning, E. Batchvarova, B. Br¨ummer, H. Jorgensen, and S. Larsen, On the extension of the wind profile overhomogeneous terrain beyond the surface boundary layer, Boundary-Layer Meteorol. , 251 (2007).[15] M. Kelly, R. A. Cersosimo, and J. Berg, A universal wind profile for the inversion-capped neutral atmospheric boundarylayer, Q. J. R. Meteorol. Soc. , 982 (2019).[16] J. G. Pedersen, S.-E. Gryning, and M. Kelly, On the structure and adjustment of inversion-capped neutral atmosphericboundary-layer flows: Large-eddy simulation study, Boundary-Layer Meteorol.153