A New Planetary Boundary Layer Scheme Based on LES: Application to the XPIA Campaign

We present a new planetary boundary layer scheme based on large‐eddy simulations for different atmospheric stability classes. Large‐eddy simulations results are compared with the wind speed measurements from the meteorological mast at the Test Centre for Large Wind Turbines at Høvsøre, Denmark. A generic formulation for the determination of mixing length scale is proposed and incorporated with the updated closure coefficients derived under realizability constraints by Temel et al. (2018, https://doi.org/10.1016/j.jweia.2018.01.002). The new planetary boundary layer scheme is implemented into the Weather Research and Forecasting model to perform mesoscale simulations for the Høvsøre test site as an idealized case as well as for a real‐data case at the eXperimental Planetary boundary layer Instrumentation Assessment campaign. Results of the idealized case reveal that the proposed scheme, VKI01, well represents the potential temperature and wind speed characteristics. It decreases mean absolute errors for most of the stability levels despite a slight overestimation for near‐neutral stable and very stable conditions. Regarding the real‐data case, a significant improvement has been achieved by the VKI01 for both turbulence kinetic energy and its dissipation rate in comparison to sonic anemometer measurements for a 2‐day period during the eXperimental Planetary boundary layer Instrumentation Assessment campaign.


Introduction
The planetary boundary layer (PBL), which is the lowest part of the troposphere, gained much attention in many scientific fields ranging from space applications (Martinez et al., 2013), over wind energy studies (Storm et al., 2009), optical turbulence (Basu & He, 2014), to weather forecasting and climate research (Flaounas et al., 2011). The parameterization of the vertical turbulent fluxes is essential for the correct modeling of the diurnal evolution of winds, temperature, and water vapor within the boundary layer (Shin & Hong, 2011;Zhang & Zheng, 2004) as well as for the large-scale atmospheric processes such as mesoscale convective rainfall systems (Jankov et al., 2005), hurricanes (Braun & Tao, 2000), and even for the radiation budgets, for example, in the case of boundary layer clouds in polar regions (Tjernström et al., 2005).
Except for the cases, where the spatial resolution is sufficiently high to use large-eddy simulation (LES) techniques, turbulent fluxes are parameterized using PBL schemes, which are based on the eddy diffusivity concept. This approach basically assumes a relation between turbulent fluxes and mean velocity gradients where this proportionality is called the eddy diffusivity. To calculate the vertical turbulent fluxes, eddy diffusivity is either determined by algebraic closures, for example, Yonsei University scheme (Hong et al., 2006), or by turbulence kinetic energy (TKE) closures, which solve a transport equation for the TKE and empirically models its sink term, the dissipation rate of TKE (Janjić, 1994). Several different TKE closures were proposed in the literature, for example, Bougeault-Lacarrere (BouLac) scheme (Bougeault & Lacarrere, 1989), Mellor-Yamada-Janjic (MYJ) scheme (Janjić, 1994), and Mellor-Yamada-Nakanishi-Niino (MYNN) scheme (Nakanishi & Niino, 2009). The main difference between these TKE closures depends mainly on the formulation of the length scale and closure coefficients for the calculation of the eddy momentum/heat diffusivities. These TKE closures are used as the one-dimensional turbulence closures by neglecting streamwise and spanwise turbulent fluxes as well as all the horizontal derivatives.
Although these schemes are suitable for mesoscale simulations, however, in case of high horizontal grid resolutions, the microscale atmospheric turbulence must be regarded as a three-dimensional process

Description of the New PBL Scheme
The new planetary boundary scheme proposed in this study is a nonlocal TKE closure. So the one-dimensional transport equation for the TKE can be defined as The source and sink terms on the right-hand side of equation (1) stand for the shear production of turbulence, P s , the buoyancy production/destruction of turbulence, P b , and the dissipation rate of turbulence, , which can be calculated as follows: In equation (1), l is the mixing length scale, and S k is the TKE diffusion coefficient. K m is the eddy diffusivity of momentum in equation (2). K h is the eddy diffusivity of heat, = 1∕300 K is the thermal expansion where S m and S h denote the stability functions of momentum and heat. One of the main differences between mesoscale PBL schemes is based on the mixing length scale formulation, which in turn is used to compute eddy diffusivities and the dissipation as can be seen in equations (4)-(6). By making use of the Høvsøre LES results presented in section 3, we propose a new mixing length parameterization. This modeling starts with the classical Blackadar (1962) formulation that is basically the harmonic average of two length scales where l 0 and l s denote the asymptotic and surface layer length scales, respectively. By this way, Blackadar's mixing length scale switches from the surface layer throughout the well-mixed boundary layer limit. This approach constitutes the basis of different PBL schemes as well, for instance, the University of Washington (UW) PBL scheme of Bretherton and Park (2009), the scale-adaptive TKE closure of Kurowski and Teixeira (2018), and MYJ scheme (Janjić, 1994). For the latter, Janjić (1994) set the asymptotic component of Blackadar's mixing length scale, l 0 , to where p s and p t are the pressures at the bottom and top of the model atmosphere, respectively, is an empirical constant to be 0.30, and q = √ k is the square root of TKE. Later on, this length scale (equation (8)) is used in the MYNN scheme by setting = 0.23 where they named the l 0 "the length scale dependent on the PBL depth." According to our proposed model, above the surface layer, asymptotic component of Blackadar's mixing length scale is expanded to two components, that is, l 0 = (1∕l k + 1∕l ) −1 . These are TKE length scale, l k , related to TKE of the PBL and asymptotic length scale, which is used to fit the mixing length scale to our reference LES data along the mixed layer. Furthermore, the mixing length scale becomes constant above the surface layer, and it asymptotically converges to a constant value since the grid resolution of mesoscale simulations falls in the mesoscale limit, Δ meso > L, where L denotes the size of energy containing eddies (Wyngaard, 2004). This is consistent with the observations of Ito et al. (2015) where they find the dissipation length scale converges to a unique value at the mesoscale limit. Considering the study of Ito et al. (2015), the mixing length scale becomes constant above the surface layer, that is, l∕ z → 0ifz > z s . Ito et al. (2015) also showed that the vertical profiles of dissipation length scale, l = B 1 l, reach to a nearly asymptotic shape if the horizontal grid resolution approaches to the mesoscale limit. In other words, as the grid resolution increases from LES limit toward the mesoscale limit, that is, where the infinity symbol denotes this asymptote. Note that the angular bracket denotes the horizontal-averaging operator applied in the horizontal directions, that is, streamwise (x) and spanwise (y), similar to Catalano and Moeng (2010). In this context, the total contribution of our l k and l tends to be constant in the mixed layer as the altitude increases similar to Ito et al. (2015). In addition to this, the sum of l k and l reach to an asymptotic value independent from the grid resolution. As discussed above, this is because the length scale profiles coincide at the mesoscale limit in accordance with the self-similarity. Since the l brings the l k toward the mixed-layer asymptote at the mesoscale limit, we call the l "asymptotic" length scale. Another essential point is to adjust the value of the B 1 coefficient. In this manner, we directly compute the dissipation, that is, = , as well as the TKE by making use of the LES. Since B 1 l is equal to k 1.5 ∕ and considering the value of B 1 l converges to a constant value (from the equation (9)), the problem itself becomes independent from the selection of B 1 at the mesoscale limit, Δ meso . From this aspect, we where l k is the TKE length scale, l is the asymptotic length scale, and l s is the surface layer length scale. If we further introduce this parameterization, the functional form of the proposed TKE length scale in equation (11) is based on Mellor and Yamada (1974) as the ratio of the first to the zeroth moment of the TKE velocity scale, q = √ 2k. This definition can be interpreted as the TKE tendency of the length scale related to changes in the atmospheric stability and the altitude. Accordingly, the TKE tends to decay with the stability and the altitude as can be seen from the vertical profiles of TKE in Figure 2c. Both k and q are able to represent this description; however, we prefer the TKE itself for its direct implementation rather than its velocity scale. Note that the closure constant, , has been selected to be 0.10 in Mellor and Yamada (1974), 0.30 in MYJ scheme (Janjić, 1994), and 0.23 in MYNN scheme (Nakanishi & Niino, 2009). In this study, we choose = where = 0.40 denotes the von Karman constant. The reason can be explained by the physical representation of l k , such that it corresponds to the center of TKE analogous to the center of gravity of a particular object. If is 1, this gives the absolute center of TKE where z → ∞. To make the l k more dominant along the middle boundary layer, the center of TKE can be moved to z = 0.40z i . This location can be determined as the center of z = 0.80z i where the middle boundary layer is still effective, before the top boundary layer starts.
According to our asymptotic length scale parameterization, it enables the mixing length scale to fit our reference LES results above the surface layer, that is, We introduce this parameterization by postulating that the asymptotic length scale be the function of the PBL height, z i , and the nondimensional stability parameter, that is, 1 = z 1 ∕L O , at the first wall adjacent cell center, z 1 , where L O is the Obukhov length scale where u * is the surface friction velocity, g∕ 0 is the buoyancy parameter, and w ′ ′ s is the kinematic heat flux at the surface. Then, similar to the UW scheme (Bretherton & Park, 2009), which is

10.1029/2018MS001580
where they adjust = 0.085 by using LES results. Note that Ri CL is the bulk Richardson number for the convective turbulent layer. Recalling our postulate given in equation (17), a curve fitting is made by the linear regression considering the quasi-linear behavior of mixing length scale above the surface layer. Supposing that the PBL height is proportional to the asymptotic length scale with the power of 1, that is, This enables l to satisfy the unit of meter. Then, l ∝ z i can be linked to the LES results through the atmospheric stability in different functional forms (logarithmic, exponential, hyperbolic tangent, etc.), that is, l =  ( , 1 )z i . Note that here is the fitting parameter. Following the UW scheme (Bretherton & Park, 2009), the exponential function provided the best agreement for our l model in our regression analysis. Beside this, since the exponential function tended to be more aware of the stability, we choose the exponential form as Then, the asymptotic length scale is defined as a piecewise function in three different stability regimes based on the classification of Gryning et al. (2007), in equation (12). By this way, the model itself tends to be adaptive to the atmospheric stability.
If the functional form of the surface layer length scale is addressed, it can be assumed to be proportional with the altitude based on Businger et al. (1971), Dyer and Hicks (1970), and Gryning et al. (2007), l s,neutral = 2.5 z where = 0.4 is the von Karman constant. This selected value matches well with our LES results for the near-neutral condition (LES-n case). On the other hand, for nonneutral conditions, we begin with the functional form of l s from MYNN scheme as l s,unstable = z(1 − a 1 ) p 1 and l s,stable = (1 + a 2 ) p 2 . Then, by means of curve fitting to the predicted LES results, we end up with the final functional form for the nonneutral conditions given in equation (13).
As mentioned earlier, a new length scale formulation is proposed to build a new PBL scheme which later can be used as a baseline for further extension, for instance, toward a 3-D PBL scheme. After applying the closure assumption for dissipation rate (equation (4)) into equation (5), then the turbulent viscosity can be redefined independent from the mixing length scale where S m ∕B 1 is the equivalent closure coefficient of C appearing in computational fluid dynamics simulations of turbulent flows by using three-dimensional turbulence closures. Based on the realizability constraints , C is set to 0.09 (Rodi & Mansour, 1993) for industrial flow applications and to 0.033 for wind engineering applications (Alinot & Masson, 2005). In this present study, B 1 and S m are taken to be 24.0 and 0.006, which is equivalent to a smaller value of C = 0.006, ensuring the realizability of the turbulence closure for further extensions as a 3-D PBL scheme. Finally, the eddy heat diffusivity coefficient is computed by utilizing the formulation of turbulent Prandtl number, that is, Pr t = S m ∕S h , proposed by Temel et al. (2018) as where m is the stability function of Monin-Obukhov theory for velocity gradient and e is the stability function for potential temperature gradient (Businger et al., 1971). Ri stands for the gradient Richardson number Moreover, similar to the MYJ scheme, the minimum value for the TKE is set to 0.01 m 2 /s 2 . For determination of the boundary layer height, the way of BouLac scheme, which is similar to the 1.5-theta-increase method  (Gryning et al., 2007) Case (Nielsen- Gammon et al., 2008), is applied. In the original BouLac scheme (Bougeault & Lacarrere, 1989), the boundary layer height is determined at the first height where the potential temperature is increased by 0.5 K for all stability levels. We keep this value for the convective atmosphere; however, it is set to 1.0 K to avoid the possible instantaneous temperature rise at the near surface for nonconvective conditions. As noted before, the present PBL scheme is a nonlocal scheme, so that the nonlocal term is computed in the following equation as suggested by Therry and Lacarrère (1983, see section 4.2, Eq. 40).

LES
Since the proposed scheme is derived based on LESs, its details will be presented in this section. For the LES experiments, OpenFOAM-based Simulator for Offshore Wind Farm Applications ABLSolver, an open source LES code written in C++, is utilized in the present study. This code solves the resolved continuity, Navier-Stokes, and energy equations based on the finite volume method. The details of governing equations are described in Churchfield et al. (2012).
In ABLSolver, second-order centered-difference is used for the spatial discretization, while the time integration is based on the Crank-Nicolson method. The Courant-Friedrichs-Lewy number is kept lower than a critical value, that is, Co < 0.375. The pressure-velocity coupling is preformed by successively solved implicit equations using pressure-implicit with splitting of operators algorithm (Issa, 1986). Velocity and pressure fields are handled implicitly to increase the stability of the numerical scheme. However, the buoyancy, subgrid-scale (SGS) viscosity, and Coriolis terms are explicitly treated (Churchfield et al., 2010).
For the SGS parameterization, despite the universal usage of Smagorinsky model which assumes a balance between shear production and dissipation terms, it is inadequate for the buoyancy-driven and transitional PBL flows (Moeng, 1984). Therefore, a prognostic equation for the SGS kinetic energy is solved besides the filtered governing equations. This model, namely Deardorff's (1980) SGS model, is used for the convective simulations. On the other hand, the effects of SGS modeling become more important for the nonconvective conditions. This arises particularly for the stable PBL due to the excess in unresolved small eddies (Kosović & Curry, 2000). We therefore use the Lagrangian-Averaged Scale-Invariant dynamic SGS model (Meneveau et al., 1996) for the neutral and stable conditions.
To derive a new PBL scheme, seven cases are investigated based on different stability classes. Atmospheric stability conditions are as reported by Gryning et al. (2007) and listed in Table 1 The initial velocity field is set to a constant value, that is, the initial wind speed, U 0 , at a certain height, h 0 , based on the atmospheric conditions at Høvsøre test site (Table 1). Then, to drive the LES, ABLSolver forces the simulation by keeping this given velocity at this desired altitude. On the other hand, the initial potential temperature field is determined relying on the LES experiments by Pedersen et al. (2012) performed at the same test site in Høvsøre, Denmark. Reference potential temperature is kept constant as ref = 280 K up to the PBL height, z i , and then increasing with the Table 2 The Details of the Computational Domain and the Grid Resolution Depending on the Atmospheric Stability atmospheric lapse rate Γ a = 0.0034 km −1 within the free atmosphere. To restrict the growing of the PBL, a three-layer potential temperature profile is described inserting a capping inversion similar to the study of where Γ i and Δh i are set to different values related to the different atmospheric stabilities, that is, Γ i ∕Δh i is 2.5 K/100 m, 1 K/50 m, and 0 K/0 m for unstable, neutral, and stable PBL, respectively.
The dimensions and grid resolutions of the computational domain are determined depending on the atmospheric stability and in accordance with the earlier studies. So the domain size and the grid resolution are 5 × 5 × 2 km 3 and 40 × 40 × 16 m 3 for the unstable PBL (Brooks & Fowler, 2012;Cancelli et al., 2014;Degrazia et al., 2009), 3 × 3 × 1 km 3 and 24 × 24 × 10 m 3 for the neutral PBL (Khani & Porté-Agel, 2017), and 0.4 × 0.4 × 0.4 km 3 and 6.25 × 6.25 × 6.25 m 3 for the stable PBL (Basu & Porté-Agel, 2006), respectively. The details of the domain size, L, grid resolution, Δ, and the number of grid points, N, are given in Table 2. The lower boundary is ideally flat but rough with different roughness heights given in Table 1.
Periodic boundary conditions are applied at the lateral boundaries. At the top surface, a slip wall boundary condition is imposed for the velocity, while the Schumann-Grötzbach wall model (Schumann, 1975) is applied for the surface shear stress and the kinematic heat flux at the lower surface. For the potential temperature, a Neumann boundary condition and a fixed gradient equal to the capping inversion gradient are imposed at the top and lower surfaces, respectively. For the unstable and neutral simulations, the kinematic heat flux is specified on the ground, whereas the surface cooling rate is imposed for the stable simulations following the model developed by Basu et al. (2008). The cooling rate is determined based on the surface kinematic heat flux (q s,ref ) values of Høvsøre test site (Table 1) through the trial and error approach. The corresponding cooling rates are set to 0.5, 0.42, and 0.36 K/hr for the very stable, moderately stable, and near-neutral stable conditions. Besides these, in order to create the turbulence rapidly within the LESs, divergence-free perturbations similar to De Villiers (2007) are included only close to the ground.
LESs in dry-air conditions are applied to seven different stability conditions based on the Høvsøre test site (Gryning et al., 2007). Simulations reach a quasi-equilibrium state at  = 21, 000 s (∼ 26 * − 30 * ) for the unstable PBL,  = 20, 000 s (∼ 18 * ) for the neutral, and  = 20, 000 s (∼ 11 * − 14 * ) for the stable PBL. Afterward, LES simulations are continued 3 * more for time averaging where * denotes large-eddy turnover time given in Table 3. It is defined to be * = z i ∕u * for nonconvective conditions while for the convective conditions as * = z i ∕w * (Ayotte et al., 1996). Herein, u * is the friction velocity, and w * denotes the Deardorff convective velocity and reads where g is the gravitational acceleration, 0 is the horizontally averaged surface potential temperature, and q s is the surface kinematic heat flux.

Table 3 Main Planetary Boundary Layer Parameters Based on LES Results
Case  (2010) is executed in the horizontal directions to determine the turbulence statistics due to the existence of quasi-steady turbulence after a characteristic time.
Once all the computations finish, to ensure the convergence, the relative change in the friction velocity, u * , Obukhov length scale, L O , and PBL height, z i , are checked. Their relative errors are lower than 5 × 10 −5 which can be considered as a satisfying convergence level. LES results are listed in Table 3 in terms of the main PBL parameters. Here u * are calculated at the first grid center close to the ground, that is, z 1/2 = 8, 5, and 3.125 m for unstable, neutral, and stable conditions, respectively. Table 3 also includes the PBL height, z i . It is specified at the height where the vertical profile of resolved heat flux becomes minimum for the unstable conditions, contrarily, it corresponds to maximum resolved heat flux in case of stable conditions. For the neutral condition, PBL height is determined at the height when the kinematic heat flux being 0. As  (Gryning et al., 2007). vu = very unstable; mu = moderate unstable; nu = near-neutral unstable; n = near-neutral; ns = near-neutral stable; ms = moderate stable; vs = very stable; LES = large-eddy simulation.
presented in Table 3, the PBL height is highest for the very unstable case and lowest for the very stable case, consistently with the literature, for example, Moeng and Sullivan (1994) and Udina et al. (2016).
First, the wind speed predictions by LES are compared with the wind speed measurements performed at the Høvsøre test site (Gryning et al., 2007). Time-and horizontal averaged wind speed profiles normalized by the reference friction velocity, u *,ref , are plotted in the vertical direction up to 400-m altitude ( Figure 1). Herein, ⟨U⟩ = (⟨u⟩ 2 + ⟨v⟩ 2 ) 1∕2 , where the ⟨u⟩ and ⟨v⟩ denotes the x and y components of the time-and horizontal averaged mean velocity. Results indicate that the LES-based wind profiles perform well and are consistent with the Høvsøre experiment. Furthermore, the normalized wind speed gradient increases from convective conditions through neutral and stable conditions, so that the strongest gradient exists once the atmosphere becomes strongly stable, that is, very stable PBL, as in Figure 1.
In Figure 2, time-and horizontal averaged wind speed and potential temperature profiles are depicted in the vertical direction depending on various atmospheric stabilities. The altitude expressed herein is normalized by the PBL height, z i . Averaged wind speed profiles show that the highest velocity gradient appears in case of the stable cases, as a result of the suppressed turbulence due to the negative buoyancy. Unlike the stable conditions, a weak velocity gradient appears for the unstable conditions thanks to the strong buoyancy and weaker winds. This leads to slightly varying wind speed in the surface layer and nearly constant wind speed until the entrainment zone consistently with the previous studies, for example Moeng and Sullivan (1994)  be expected, in comparison to stable cases, wind speed gradients are lower for the neutral case, revealing the well-reported change of wind gradient from very stable conditions through weakly stable condition (Huang & Bou-Zeid, 2013). Furthermore, the averaged potential temperature profiles exhibit the typical behaviors related to the stability conditions, such as positive and negative potential temperature gradients on the surface for stable and unstable conditions, as well as the capping inversion at the top of PBL at z∕z i = 1.0.
One of the most important features of the atmosphere is its TKE that is the basis for mesoscale parameterizations of turbulence effects in the PBL where k is the total TKE and u ′ i is the fluctuating velocity. Herein, E ( ) denotes the energy spectrum in the wave number space, , and such that the integration of E ( ) d along the wave number space also gives the total kinetic energy. In Figure 2, the vertical variations of the total and SGS TKE are plotted for different stability conditions. It can be stated that the contribution of SGSs are negligible compared to resolved scales except near the ground. Another observation is that the TKE decays with the atmospheric stability, such that it is quite high because of the strong convective motions during the daytime. But it decays from the most unstable condition to neutral condition and toward the most stable condition. This is in fact the observation Figure 3. Time-and streamwise-averaged energy spectra versus spanwise wave number, 2 , for very unstable (top), near-neutral (middle), and very stable (bottom) conditions at two different altitudes, that is, z∕z i = 0.05 and z∕z i = 0.80. Black, blue, and red solid lines denote the averaged energy spectra of resolvedũ 1 ′u 1 ′,ũ 2 ′u 2 ′, andũ 3 ′u 3 ′ fluctuations, respectively. Gray dashed line denote the theoretical Kolmogorov −5/3 slope. that gives an idea to model the turbulence length scale later on. By this way, TKE and its length scale can be linked similar to the modeling of MYNN scheme by Nakanishi and Niino (2009).
To ensure whether the spatial resolution is sufficient for well-resolved LES, the energy spectra is investigated on the basis that if it exhibits for a −5/3 cascade, that is, E ( ) ∝ −5∕3 , within the inertial subrange in accordance with Kolmogorov (1941) −5/3 theory. Therefore, the one-dimensional energy spectral densities (or energy spectra) of resolvedũ ′ 1 u ′ 1 ,ũ ′ 2 u ′ 2 , andũ ′ 3 u ′ 3 fluctuations are computed along the spanwise y direction, that is,  Figure 3, the time-and streamwise-averaged energy spectra along the spanwise wave number, that is, 2 , are presented for the very unstable, neutral, and very stable conditions at two different vertical locations, that is, z∕z i = 0.05 and z∕z i = 0.80. The energy spectra represent the spectral energy transfer from larger scales, that is, low wave numbers, to the dissipative smaller scales, that is, high wave numbers, via the energy cascading process. In Figure 3, low wave numbers correspond to the high-level spectral energy because of the TKE production by larger scales, while it decays steeply through the high wave numbers due to the turbulence dissipation. In addition to this, all the LES computations exhibit the Kolmogorov spectral −5/3 slope within the inertial subrange, which enables us to ensure the reliability of the LES. Besides that the averaged energy spectra are much higher at the altitude z∕z i = 0.05 than the altitude at z∕z i = 0.80 for all the conditions, this is because the fact that the TKE is decaying toward the higher altitudes as depicted in Figure 2. Energy spectral density is also effected by the effect of atmospheric stability, such that the spectra within the inertial subrange decrease from very unstable condition, that is, < E ii > z∕z i =0.05 ∼ (10 1 ) and < E ii > z∕z i =0.80 ∼ (10 0 ), toward very stable condition, that is, < E ii > z∕z i =0.05 ∼ (10 0 ) and < E ii > z∕z i =0.80 ∼ (10 −1 ). And finally, the highest wave number increases from the very unstable condition, that is, 2 ∼ 8e −2 , to the very stable condition, that is, 2 ∼ 5e −1 , because of the fact that we employ higher horizontal grid resolution for the stable atmosphere, that is, Δ x,y = 6.25 m, than the convective atmosphere, that is, Δ x,y = 40.0 m.
One of our main aims is to determine the length scale depending on the atmospheric stability to derive a new PBL scheme. The length scale can be obtained from the turbulence dissipation and TKE, employing according to Kolmogorov's hypothesis, by making use of the equation (4). Then, it becomes possible to determine the length scale, that is, l = k 1.5 ∕( B 1 ), by getting TKE and dissipation fields from our LESs. Before this, we investigate the budget of total TKE in order to see the relative importance of each term in this equation and to judge the effect of the atmospheric stability. To do this, as presented in the Appendix A, the one-dimensional total TKE budget equation is derived as follows: where  denotes the advection, T T , T P , and T sgs denote the turbulence, pressure, and SGS transport, respectively. P S and P B are the shear and buoyancy productions, and stands for the dissipation. Note that the advection term is negligible in equation (28) since the horizontally averaged ⟨w⟩ = 0 as a consequence of the continuity equation. Since we aim to get the dissipation field for the determination of the length scale, we do not followed the modeled expression for the , which is based on Kolmogorov hypothesis in equation (4) and common in the literature, for example, Sullivan and Patton (2011), Ramachandran and Wyngaard (2011), and Kosović and Curry (2000). Instead, we determine as the double inner product of two second-rank tensors which are the fluctuating components of the deviatoric SGS stress tensor and the strain-rate tensor, =< D′ i S ′ i >. Following the equation (28), the TKE budget of very unstable, near-neutral, and very stable PBLs are presented in Figure 4. It shows the relative contribution to the total TKE term by term along the vertical direction. The positive and negative values of the TKE budget indicate if the TKE is produced or destructed, respectively. Figure 4 reveals that the transport terms seem to be less important than the production and dissipation terms. Its effect is only evident in case of the convective condition. This is due to the transport of strong updraughts and downdraughts different from the neutral and stable conditions. And the integral of the transport term over the altitude goes to ∼0 consistently with the earlier studies, and the overall behaviors of both pressure and turbulent transport terms are in well agreement with the LES results of Moeng and Sullivan (1994) for the convective and neutral PBL and Jiménez and Cuxart (2005) for the stable PBL. The main contribution to the TKE arises from the existence of production and dissipation terms. Regarding the production term, the buoyancy production is the most sensitive term to the atmospheric stability. To illustrate, it becomes positive because of the upward buoyancy in convective conditions. It makes a first peak around z∕z i = 0.1 and decreases until the PBL height, z∕z i = 1.0, where it makes the second peak which is relatively small and negative because of the existence of stably stratified capping inversion. This feature of the buoyancy production is a well-known feature of the unstable atmosphere (Catalano & Moeng, 2010;Moeng & Sullivan, 1994). Under the convective condition, the buoyancy production dominates the total TKE variation compared to the shear production that is only evident close to the ground due to the existence of wind shear; however, it disappears above the surface layer, that is, z∕z i = 0.1. After the sunset, convective motions dramatically weaken, and the stable atmosphere takes place. The buoyancy production is no longer a TKE production in this condition, and it pretends to be weakly destructive as can be seen in Figure 4. If we look at the neutral atmosphere, there is neither a buoyancy production nor a buoyancy destruction; instead, it is characterized by the high shear production due to the existence of strong wind shear. Besides, the final observation is that the dissipation, , behaves as the most dominant source for the destruction of TKE. Figure 4 also includes the vertical profile of the dissipation up to 1.5 km for all the stability classes. The dissipation takes its highest value close to the Earth's surface, and then it goes down with the altitude since the dissipation mainly arises close to the ground where the mechanical energy transforms into the thermal energy produced by the smallest scales. Regarding the effect of the atmospheric stability, the dissipation decreases from unstable condition to neutral and stable conditions above the surface layer; however, near the ground, this tendency become reversed. The main finding about the dissipation variation is that it decreases with the atmospheric stability. For the unstable atmosphere, the dissipation becomes less from very unstable to near-unstable condition. For the stable atmosphere, similar to the observation of Huang and Bou-Zeid (2013), an increase in stability leads to a lower dissipation value.
Finally, the results of the proposed length scale model are depicted in Figure 5 in comparison to the Høvsøre LES results. The proposed length scale matches well with these LES results until 0.9z i for all the stability conditions. However, above this altitude, there appears a large variation in LES profiles, particularly for the unstable conditions. This is because our model does not presently account for the entrainment and inversion effects that exist between the upper PBL and free atmosphere. Therefore, for the altitudes higher than 0.9z i , the mixing length scale is determined based on the MYJ scheme (Janjić, 1994), as follows: where Δz corresponds to the vertical grid spacing.

Idealized PBL Simulations
In this section, we present the results of idealized PBL simulations in comparison to LES results, which are presented in section 3. The idealized PBL simulations are performed by using an open-source numerical weather prediction code, that is, the WRF model version 3.5.1 (Skamarock et al., 2008). The model forcing is conducted by imposing the surface heat fluxes presented in Table 1. However, the momentum flux is calculated by using a surface layer scheme, MM5 similarity scheme (Beljaars, 1995). In this section, we make use of LES results to verify our proposed model. The concept of using high-fidelity simulations, such as LES, as a verification tool for PBL parameterization schemes, has been widely applied in the literature when a novel turbulence closure is initially designed (Bretherton & Park, 2009;Efstathiou & Beare, 2015;Nakanishi & Niino, 2009). Even if the LES and mesoscale models used for this comparison are different, a recent LES intercomparison study by Mirocha et al. (2017) has shown that both dynamical cores (Simulator for Offshore Wind Farm Applications and WRF) provide comparable results in canonical ABL simulations. As described in section 2, it is possible to state that our proposed PBL scheme has some similarities both with MYNN and BouLac schemes. The mixing length scale parameterization is based on three components as it is in case of MYNN scheme, and the closure coefficient to determine turbulent viscosity is constant similar to the BouLac scheme. Therefore, idealized PBL simulations are also performed by using MYNN (Nakanishi & Niino, 2009) and BouLac schemes (Bougeault & Lacarrere, 1989), enabling us to evaluate the performance of our model against these two state-of-the-art PBL schemes. As far as the computational details are concerned, the horizontal domain size is set to L x,y = 15 km × 15 km for unstable PBL, L x,y = 9 km × 9 km for neutral PBL, and L x,y = 4 km × 4 km for stable PBL. The horizontal grid resolution is taken to be 1 km in order to avoid passing the gray zone. It is worth noting that the horizontal length of the mesoscale domain does not affect the results of idealized PBL simulations due to the idealized conditions. The vertical height is set to 2 km where 80 grid points are located along this direction, by clustering near the ground. Time integration is conducted setting time step to Δt = 5 s. The total simulation time is same as our LES computations, that is,  = 23, 400 s for unstable condition and  = 24, 000 s for neutral and stable conditions. Note that periodic boundary conditions are applied to lateral boundaries and the Coriolis force is accounted consistently with our Høvsøre LES database.
As noted before, the idealized PBL simulations are performed for the same stability conditions presented in Table 3. Compared to the LES results, vertical profiles of horizontally averaged potential temperature and wind speed, predicted by VKI01, MYNN, and BouLac schemes, are presented in Figures 6 and 7. As presented in Figure 6, within the mixed layer of unstable PBLs, VKI01 and BouLac provide better predictions than the MYNN scheme. However, within the surface layer, potential temperature gradient is overestimated by VKI01 scheme compared to both schemes. Considering the mixed layer characteristics, VKI01 presents the best agreement in moderately unstable PBL where BouLac significantly overestimates the potential temperature. However, BouLac is slightly better than VKI01 for very and near-neutral unstable conditions. MYNN scheme underestimates the potential temperature in the mixed layer where it underestimates the potential temperature. In Table 4, this behavior is quantified with a performance metric, that is, mean absolute error (MAE), against to the Høvsøre LES results. VKI01 produces the lowest error of potential temperature for moderately unstable PBL, while BouLac is the best in case of very and near-unstable PBLs. For the stable PBL, although VKI01 and MYNN schemes perform very similarly, VKI01 stands slightly better than MYNN below the PBL height except for near the ground. While VKI01 improves near-neutral stable PBL with the lowest MAE, MYNN is slightly better than VKI01 for moderately stable and very stable PBLs. For the BouLac scheme, it significantly overestimates the potential temperature all along the stable boundary layer. Note that none of the PBL schemes is able to predict correctly the capping inversion, showing the intrinsic limitations of these PBL schemes for stably stratified boundary layer. Regarding the near-neutral condition, VKI01 performs superior than other two PBL schemes considering its mixing layer and PBL top behavior. Its MAE is equal to 0.142, which is much lower than that of MYNN and BouLac schemes. The latter fails predicting correctly the potential temperature both below and above the PBL for the neutral conditions (Figure 6g). Along with the potential temperature profiles, predicted wind speed profiles for three PBL schemes in comparison to LES results are also given in Figure 7. For all the stability levels, BouLac scheme fails to represent a good agreement with LES results. However, although VKI01 scheme overestimates the potential temperature gradient near the ground for convective cases, it provides better estimation for the wind speed compared to the MYNN scheme. VKI01 performs better also through the convective mixed layer. Quantitatively, VKI01 gives the lowest mean absolute errors (MAEs) for all the stability levels except for near-neutral and very stable PBLs where MYNN improves the wind speed slightly better than VKI01.

Application to the XPIA Field Campaign
In the previous section, the proposed model is verified with the LES results to some degree that despite its simpler closure coefficient formulation, the new PBL scheme is able to provide similar predictions to MYNN scheme, which solves a set of algebraic equations to determine the closure coefficients. However, the validity of the newly proposed PBL scheme for a real case study should also be tested based on a full-scale field campaign. In this section, we present the mesoscale simulations for the XPIA campaign (40.05 • N, 105 • W) conducted at the Boulder Atmospheric Observatory (BAO) from 2 March to 31 May in 2015 (Lundquist et al., 2017). The BAO tower contains a total of 12 three-dimensional sonic anemometers (Campbell CSAT3) mounted on northwesterly (NW, 334 • ) and southeasterly (SE, 154 • ) pointing booms every 50 m starting from 50 m above ground level up to 300 m. This experimental campaign involves high-frequency 20-Hz measurements of wind, temperature, and relative humidity. Recently, Muñoz-Esparza et al. (2018) derived the turbulence dissipation rate, , within the PBL using these sonic anemometer measurements. These results are available at z = 50 m, and they are classified into daytime/convective and nighttime/stable conditions. These measurements will help us to validate mesoscale simulations, performed by the new PBL scheme, and again compare with MYNN and BouLac schemes. The mesoscale domain depicted in Figure 8 consists of two one-way nested domains, the largest domain has a horizontal resolution of 27 km with a domain size of 3,375 × 2,700 km 2 . The smallest domain, on the other hand, has a horizontal resolution of 9 km with a domain size of 1,332 × 1,197 km 2 . For the vertical discretization, 53 grid points are used with the finest resolution of 2.5 m close to the ground, while the upper boundary is located at 15 km. Boundary and initial conditions are extracted from the National Center of Environmental Prediction North American Regional Reanalysis data sets. For the PBL parameterization, BouLac, MYNN, and VKI01 schemes are used. In addition to the PBL schemes, the other physical parameterizations utilized for the WRF setup are the longwave and shortwave Rapid Radiative Transfer Model (Iacono et al., 2008) for the radiation, the Thompson microphysics parameterization (Thompson et al., 2008) and MM5 Similarity Scheme (Beljaars, 1995) for the surface layer parameterization, and the Community Land Model version 4 (Lawrence et al., 2011) for the land surface coupling.
Nested mesoscale simulations are performed during the 6-9 March 2015 at the XPIA field campaign for 72 hr, but the first 24 hr is discarded considering the model spin-up, so that the results are presented between the  7 and 9 March. In Figure 9, mesoscale simulations are compared against to sonic anemometer measurements at z = 50 m. Figure 9a and 9b present the time evolution of the TKE and dissipation rate, respectively. Compared to the MYNN and BouLac schemes, VKI01 matches well with the XPIA measurements presenting the best agreement where it is able to reproduce the diurnal evolution of the boundary layer. Both MYNN and BouLac significantly underestimate the TKE and dissipation rate profiles during the first (t = 7.0-7.5 days) and second (t = 8.0-8.5 days) nighttime. Such that VKI01 and XPIA profiles are ∼ (10 −2 ) for TKE and ∼ (10 −4 ) for the dissipation rate. However, TKE for MYNN and BouLac are ∼ (10 −3 ) and ∼ (10 −4 ) where dissipation rate for MYNN and BouLac are ∼ (10 −5 ) and ∼ (10 −7 ), respectively. Although MYNN and BouLac schemes predict the TKE better during the daytime, VKI01 is much better for the dissipation rate in case of convective conditions. In Table 5, in order to evaluate the performance of simulated PBL schemes, MAE for TKE, its dissipation rate, temperature, and wind speed are presented against to XPIA measurements from Muñoz-Esparza et al. (2018). VKI01 reproduces the dissipation rate better by reducing MAE compared to MYNN and BouLac schemes. However, for the MAE of k, MYNN scheme produces slightly lower error than both PBL schemes. On the other hand, Figures 9c and 9d show the time evolution of wind speed and temperature at z = 50 m for three PBL schemes. Even though the wind speed profiles are able to capture the general behavior of the diurnal pattern at the XPIA campaign, all the PBL schemes, particularly VKI01, overpredict the wind speed during nighttime conditions. While they underestimate it during the daytime condition (i.e., t = 7.7-8.2 days). But in general, BouLac scheme predicts the wind speed much better. In Figure 9d, the time evolution of the temperature acquired by three PBL schemes are depicted. It can be stated that all the PBL schemes are in good agreement except for the first night time till the first morning transition, that is, t = 7.0-7.5 days. Eventually, if the MAEs are concerned, BouLac scheme improves both wind speed and temperature resulting in the lowest errors as MAE U = 1.14 m/s and MAE T = 1.04 • C, respectively.
In Figure 10, the time variation of the wind direction and vertical velocity are depicted compared to the sonic anemometer measurements at z = 50 m. Unlike the wind speed trend, VKI01 matches well with the XPIA measurements along the whole period except for the peak around t = 7.8 days. Moreover, both stable and convective conditions as well as the morning/evening transitions in terms of the wind direction are well represented compared to the sonic measurements. On the other hand, MYNN scheme causes large biases for the wind direction at the first nighttime and around the first noon as well as during the morning and evening transitions. Even these biases appear similarly for the BouLac scheme, during the first nighttime and the morning/evening transitions, its overall performance in terms of MAE is slightly better than VKI01 (Table 5). In addition to these, all the PBL schemes are able to reproduce the wind direction during the second morning transition. Figure 10d indicates the time evolution of the vertical wind component at z = 50 m. This result shows that all the simulated PBL schemes agree each other; however, the measured vertical speed in the XPIA campaign is much higher than these predicted results. Figure 11 compares three PBL schemes to get more insights into their PBL characteristics. Figure 11a presents the time evolution of friction velocity at the surface. Despite the general agreement among the three PBL schemes, there are some unrealistic trends for the friction velocity. To illustrate, MYNN scheme continuously oscillates on an hourly basis for the whole period, and it approaches to zero during morning/evening transitions. During the second nighttime, both MYNN and BouLac make a peak where this peak also exists in the wind speed profile, as can be seen in Figure 9c. Considering the overall behavior of u * , the evolution of VKI01 takes place falling in the between MYNN and BouLac schemes. Figure 11b shows the PBL height

10.1029/2018MS001580
variation where there is an agreement during the stable conditions. However, for the first nighttime, MYNN makes three instantaneous peaks until z i = 1, 000 m that becomes unrealistic for the stable atmosphere. Another observation can be stated regarding the value of PBL height during the convective condition. Even the results of MYNN and BouLac are consistent at the first daytime as z i ∼ 700 m, VKI01 reaches to such a higher altitude as of z i ∼ 2, 000 m. However, these high difference disappears during the second daytime. In Figure 11c, time variation of the temperature at z = 2 m is demonstrated. Similar to the temperature profiles at z = 50 m, as shown in Figure 9d, the results acquired by PBL schemes at z = 2, more generally, are in good agreement between each other. Eventually, the final result is depicted in Figure 11d, which is the time evolution of the mixing length scale at z = 300 m. VKI01 matches generally well with the MYNN scheme in comparison to the BouLac scheme. Furthermore, the average value of the l is around 10 m for the stable condition, while it rises to 50 m for the convective condition. These values are in consistency with that obtained from the Høvsøre LES results, as can be seen in Figure 5. Accordingly, it can be argued that the developed PBL scheme produces notably realistic mixing length scale profiles depending on the atmospheric stability getting more confidence to the ability of the proposed scheme. On the other hand, MYNN and BouLac models converge to the zero during the stable condition, which stands for quite unrealistic. Another unrealistic observation is that the BouLac model gives a high value of length scale corresponding to the 4 times of the VKI01 and MYNN schemes.

Conclusion
In the present study, we presented a new PBL scheme, VKI01. Our proposed scheme incorporates with a new formulation for the mixing length scale, which has been applied to idealized and real test cases, and benchmarked against two state-of-the-art PBL schemes (MYNN and BouLac PBL schemes).
For idealized cases, mesoscale simulations with a horizontal grid resolution, Δ meso = 1 km, have been performed for three PBL schemes in the WRF framework. Results have been compared against reference LESs performed over the Høvsøre test site in Denmark (Gryning et al., 2007). To account for the effect of atmospheric stability over a significant range, seven different stability levels have been employed from very unstable regime through the very stable regime based on the classification of Gryning et al. (2007).
Considering the results of idealized mesoscale simulations, the proposed scheme performs overall better in terms of predicting horizontally averaged wind speed and potential temperature profiles. For the convective atmosphere, VKI01 gives better predictions along the convective mixed layer. However, it overestimates the potential temperature inside the surface layer where the MYNN and BouLac schemes are tend to be better. On the other hand, for the wind speed predictions, VKI01 is superior for all the stability levels except for near-neutral stable and very stable PBLs where MYNN acts to be slightly better (Table 4).
Regarding the real case test, one-way nested mesoscale simulations with the horizontal grid resolutions of 27 and 9 km have been performed for the XPIA campaign at the BAO in 2015. To evaluate the skill of the proposed PBL scheme in a real application, its performance have been compared against to the sonic anemometer measurements at the altitude of z = 50 m during 7-9 March. Our proposed scheme appears to exhibit improved skill in predicting TKE and its dissipation rate except for the overestimation of TKE during convective conditions. If the MAEs are concerned, VKI01 produces the lowest error of , whereas that of k is slightly higher than both schemes due to the overprediction in the convective regime. However, statistically robust analysis needs to be performed to further confirm this interpretation. In addition to k and , time evolution of wind speed and temperature profiles have been investigated. All the performed PBL schemes are able to reproduce the general behavior of XPIA measurements for the temperature, while some large biases have been observed for the wind speed. However, MYNN and BouLac schemes bring an advantage over VKI01 estimating them better. To illustrate, BouLac scheme improves both wind speed and temperature resulting in the lowest errors. Unlike the wind speed and temperature profiles, in case the evolution of wind direction is concerned, VKI01 has a good agreement with the XPIA data together with BouLac scheme.
Regarding the general PBL characteristics, all the PBL schemes fail to predict the PBL height during the stable conditions because the predicted z i is around ≈10 m for all the schemes. This is due to the failure in the algorithm that predicts the z i . We believe that the way of the calculation for the PBL height needs to be a further investigation as a future study. Nevertheless, MYNN scheme makes three instantaneous peaks until z i = 1, 000 m that stands unrealistic for the stable condition. For the convective, MYNN and BouLac are consistent with each other having z i ∼ 700 m at the first daytime while VKI01 rising to such a higher altitude as of z i ∼ 2, 000 m. This high difference in terms of the PBL height disappears during the second daytime.
Besides the trend of PBL height, time evolution of the mixing length scale at z = 300 m reveals that the VKI01 agrees well with the MYNN scheme in general. It can also be stated that the developed VKI01 PBL scheme leads to a realistic mixing length scale profiles depending on the atmospheric stability, providing more confidence with regard to the proposed scheme. This is because MYNN and BouLac models unrealistically converge to the zero during the stable conditions. As well as the BouLac schemes result in a quite-high length scale values, corresponding to the 4 times of the VKI01 and MYNN schemes. In this context, even if the proposed PBL scheme performs well in predicting k and except for the slight overestimation during daytime, the BouLac scheme seems better in terms of wind speed and potential temperature estimations. One reason may be the use of an existing parameterization, that is, MYJ scheme, above z∕z i > 0.9 (see equation (29)), while presently, the proposed mixing length scale does not take into account the interaction between the upper PBL and the free atmosphere. The second reason for this mismatch might be caused by the parameterization of the nonlocal term given in equation (24) which is based on Therry and Lacarrère (1983). Moreover, another possibility for the mismatch between the observations and the predicted results might be resulted from either the definition of the PBL height or the physical parameterizations used by WRF, such as the microphysics, surface layer, and land-surface coupling schemes. All these possibilities will be investigated in detail as a future study.
As a result, considering both idealized and real case tests, the developed PBL scheme gives promising predictions despite some mentioned exceptions. It is possible to state that our objective has not been to propose an ultimately improved PBL scheme but to develop a new formulation that can be considered as a baseline for further studies to develop a three-dimensional PBL scheme. Our proposed PBL scheme is based on the use of constant closure coefficients to determine eddy diffusivities, as opposed to the MYNN scheme. In this way, our proposed formulation stands as a baseline for a more straightforward extension to a three-dimensional PBL parameterization. However, it must be noted that the development of a three-dimensional parameterization will require further LES studies in order to quantify the contribution of horizontal turbulent fluxes in both horizontally homogeneous PBLs as well as flow over complex terrain, including various stability conditions. Finally, these idealized PBL simulations are performed in order to verify the implementation of our new PBL scheme. And, as expected, our model performs mostly better than the other schemes thanks to the fact that its mixing length definition is derived from the LES data set given in section 3, which is also used to be compared. Therefore, the present results do not discard the applicability of the other PBL schemes. In order to perform a quantitative realistic intercomparison study of various PBL schemes, long-term mesoscale simulations must be performed and compared with near-surface meteorological observations (Hu et al., 2010), which is beyond the scope of the present study.

10.1029/2018MS001580
where ⟨K⟩ denotes averaged resolved kinetic energy that can be decomposed into the mean field resolved kinetic energy, ⟨K⟩ , and the resolved TKE,k ′ , that is, After applying horizontal averaging to the equation (A1), that is, where the brackets, that is, ⟨·⟩, denote the horizontal-averaging operator. The second term is moved to the right-hand side and expanded in two terms (I and II). After multiplying by ⟨ũ i ⟩, that is, Now we will make further arrangements starting with the first term If we continue with the III-V terms, respectively, Substituting all the derived components into equation (A7) gives the mean field resolved kinetic energy equation, Eventually, horizontally averaged resolved TKE equation forms by subtracting ⟨K⟩ (equation A12) from ⟨K⟩ (equation A4), that is, Substituting equation (A14) into equation (A13) gives the final form of the resolved TKE budget equation, that is, where  denotes the advection, T T , T P , T sgs denote the turbulence, pressure, and SGS transport, respectively. P S and P B are the shear and buoyancy productions, and stands for the dissipation.