## Abstract

This is the second part of a study about turbulence and vertical fluxes in the stable atmospheric boundary layer. Based on a suite of large-eddy simulations in Part I where the effects of stability on the turbulent structures and kinetic energy are investigated, first-order parameterization schemes are assessed and tested in the Geophysical Fluid Dynamics Laboratory (GFDL)’s single-column model. The applicability of the gradient-flux hypothesis is first examined and it is found that stable conditions are favorable for that hypothesis. However, the concept of introducing a stability correction function *f _{m}* as a multiplicative factor into the mixing length used under neutral conditions

*l*is shown to be problematic because

_{N}*f*computed a priori from large-eddy simulations tends not to be a universal function of stability. With this observation, a novel mixing-length model is proposed, which conforms to large-eddy simulation results much better under stable conditions and converges to the classic model under neutral conditions. Test cases imposing steady as well as unsteady forcings are developed to evaluate the performance of the new model. It is found that the new model exhibits robust performance as the stability strength is changed, while other models are sensitive to changes in stability. For cases with unsteady forcings, which are very rarely simulated or tested, the results of the single-column model and large-eddy simulations are also closer when the new model is used, compared to the other models. However, unsteady cases are much more challenging for the turbulence closure formulations than cases with steady surface forcing.

_{m}## 1. Introduction

The representation of the stably stratified atmospheric boundary layer (SABL) in large-scale atmospheric simulations remains an open challenge with significant impacts on the forecasting ability of weather and climate models, especially near the earth’s surface. Turbulent mixing in the SABL directly influences the thermodynamic conditions of the land surface and the magnitude of the fluxes between the atmospheric boundary layer (ABL) and the free troposphere. Large-scale weather forecasting and climate models, particularly operational ones, normally obtain a globally good performance (e.g., Derbyshire 1999a; Viterbo et al. 1999) and avoid runaway surface cooling (Mahrt 1998) by artificially enhancing turbulent mixing in the SABL. Overestimating mixing efficiency leads to a spuriously deep ABL (e.g., Seidel et al. 2012), while an underestimation of mixing efficiency results in a strong vertical temperature gradient as well as local high concentrations of contaminants. These potential errors associated with mixing efficiency then interact with other aspects of the earth system and lead to errors at large temporal and spatial scales (e.g., Bintanja et al. 2012; McNider et al. 2012). Therefore, it remains crucial to develop and test a more realistic and robust representation of the SABL with good performance under all stabilities.

Most current representations of the SABL in atmospheric models rely on the gradient-flux hypothesis, which represents a broad class of turbulent transport models. In these models, turbulent fluxes are related to the local gradients of the corresponding mean quantities by eddy viscosities or diffusivities. For the vertical turbulent transfer of momentum and heat, these relationships can be written as

where *U*, *V*, and Θ are the mean streamwise velocity, mean spanwise velocity, and mean potential temperature, respectively; *u*′, *υ*′, and *θ*′ are turbulent departures from their corresponding mean values; the angle brackets represent the operation of horizontal and/or temporal averaging at quasi-equilibrium state; and *K _{m}* and

*K*are the eddy viscosities for momentum and heat, respectively. One popular parameterization scheme in research and operational models uses first-order closure directly derived from the Monin–Obukhov similarity (MOS) theory (e.g., Garratt 1992, 49–52) to model the effects of stability on eddy diffusivities. Although the MOS theory was originally proposed for the surface layer, the closely related local similarity theory (Nieuwstadt 1984) has extended its use to the entire ABL under stable stratification. The local similarity theory assumes that the gradient Richardson number

_{h}where *g* is the gravitational acceleration, is constant at a sufficient distance from the surface in the SABL. A local *K _{m}* can then be formulated as follows:

where *l _{N}* is the mixing length under neutral conditions and

*f*is a stability correction function, which can be expressed as

_{m}While Eq. (6) generally yields good results for ABL observations under weakly stable conditions (Derbyshire 1999b), it introduces a threshold value Ri_{g} = 0.2, beyond which turbulent mixing is totally shut down. This is not in agreement with observations (e.g., Mahrt and Vickers 2005) and numerical simulations (e.g., Zilitinkevich et al. 2008) where turbulent mixing is shown to persist for Ri_{g} well above 0.2. To remedy this shortcoming, the so-called sharp form, which is an empirical relaxation based on Eq. (6), has been proposed to eliminate the upper limit of Ri_{g} (e.g., King et al. 2001):

As mentioned earlier, large-scale operational models typically require more mixing than what is suggested by observations; this is often accomplished by imposing an *f _{m}*, which decays more slowly with stability than Eqs. (6) and (7). A typical implementation is the long-tails function (Beljaars and Holtslag 1991) written as

In addition to Eqs. (6)–(8), there are many other forms of *f _{m}* used in research and operational climate models (e.g., Cuxart et al. 2006). Given the popularity of the first-order closure parameterization using

*l*and

_{N}*f*, we investigate in this paper their adequacy in representing turbulent fluxes in the SABL, with the aim of developing a general parameterization that does not have an ad hoc critical Richardson number or artificially enhanced diffusivity, and that agrees with our current understanding of the SABL under various regimes. However, we note that higher-order closures have also been developed and tested. Canuto et al. (2008) for example developed a second-order closure model with no critical Richardson number (see also references therein).

_{m}To that end, we resort to a high-resolution large-eddy simulation (LES), which has been described by Chimonas (1999) to provide the best hope for synthesizing a “realistic model” of the SABL. Beare et al. (2006) conducted an intercomparison study of LES and concluded that LES is able to reliably simulate quasi-equilibrium, moderately stable cases, but also pointed out that subgrid-scale (SGS) models and resolution do play important roles in the accuracy of the results. Using a state-of-the-art LES code with a scale-dependent Lagrangian dynamic SGS model (Meneveau et al. 1996; Porté-Agel et al. 2000; Bou-Zeid et al. 2005, 2008), Part I of the present paper (i.e., Huang and Bou-Zeid 2013) further expanded the stability range and systematically investigated the effects of stability on the bulk dynamics, turbulent structures, and TKE budgets, as well as the applicability of the local similarity theory in the SABL. In Part I we found that 1) the vertical extent of turbulent structures is reduced with increasing stability, 2) buoyant destruction of turbulent kinetic energy becomes more important than viscous dissipation under the highest stabilities, and 3) the *z*-less range of scaling in the SABL starts at lower heights than previously anticipated. As Part II of Huang and Bou-Zeid (2013), the present study mainly focuses on the representation of the SABL in large-scale atmospheric models. We use the same LES dataset as in Huang and Bou-Zeid (2013), which only includes steady-state surface forcings, to conduct a priori tests of first-order eddy-diffusion models and their dependence on stability. Based on these steady-state tests, we propose a novel mixing-length model. Subsequently, this new model is implemented in the Geophysical Fluid Dynamics Laboratory (GFDL)’s single-column model (SCM) (Anderson et al. 2004) and systematically evaluated using the steady-state test cases, as well as new cases with unsteady surface forcings.

This paper is structured as follows. In section 2, the GFDL SCM and the test cases are briefly introduced. In section 3, we examine the applicability of the gradient-flux hypothesis across varied stabilities. The new mixing-length model is then introduced in section 4. We run test cases with steady and unsteady surface forcings to assess the performance of the new model using the GFDL SCM in sections 5 and 6, respectively. Finally the paper is summarized in section 7.

## 2. Methodology

### a. The GFDL SCM

The SCM we use in this study is based on the GFDL atmospheric model, version 2.1 (AM2.1). The model solves the Reynolds-averaged equations for momentum and thermal energy conservation, including the Coriolis force. In our current simulations, we ignore radiative flux divergence and water phase changes. A comprehensive description of AM2.1 can be found in Anderson et al. (2004); here, we give the relevant details about the surface flux module and the current parameterization scheme of the SABL. The Monin–Obukhov similarity relationship is used as the wall model only at the first grid point close to the surface in the stable surface layer following

where *M* is the horizontal wind speed defined as , *κ* is the von Kármán constant, is the Obukhov length, *u*_{*} is the friction velocity, *θ*_{*} is the surface temperature scale, and Θ_{s} is the surface potential temperature. Following the recommendations of the Global Energy and Water Cycle Experiment (GEWEX) Atmospheric Boundary Layer Study (GABLS) SCM intercomparison study (Cuxart et al. 2006), the coefficients *β _{m}* and

*β*are set as 4.8 and 7.8, respectively. The stability correction function in the GFDL SCM is a modified version of the sharp form given in Eq. (6); it is largely relaxed in the ABL at high stabilities to prevent the laminarization of the flow and the decoupling of the temperature of the land surface from that of the atmosphere, and takes the following form:

_{h}In the model, Eqs. (11) and (6) are blended in the SABL to provide a smooth transition between the ABL where Eq. (11) is used to the free troposphere where Eq. (6) is applied. However, since we are only concerned with the ABL, we use only Eq. (11) for simplicity. The vertical resolution of the SCM runs is 2.5 m and the time step is 180 s.

### b. Case description

All test cases used in this study are based on the case described in the GABLS project (Beare et al. 2006; Cuxart et al. 2006), which simulates a dry clear-air SABL driven by a moderate surface cooling rate (i.e., −0.25 K h^{−1}) in the Beaufort Sea Arctic Stratus Experiment (BASE). However, in this study, we increase the stabilizing surface cooling rate significantly beyond the GABLS value. The initial mean potential temperature is 265 K up to 100 m, with an overlying inversion of strength 0.01 K m^{−1}. A constant geostrophic wind of *U _{g}* = 8 m s

^{−1}is imposed, and the Coriolis parameter

*f*is set to 1.39 × 10

_{c}^{−4}s

^{−1}, corresponding to latitude 73°. Stress-free and no-penetration conditions are imposed at the top of the computational domain—that is, ∂

_{3}

*u*

_{1,2}=

*u*

_{3}= 0, where 1, 2, and 3 (or

*x*,

*y*, and

*z*) refer to the streamwise, spanwise, and vertical directions, respectively. As described in Part I, we study six test cases with a range of steady surface cooling rates—namely, −0.25 K h

^{−1}(case A), −0.5 K h

^{−1}(case B), −1 K h

^{−1}(case C), −1.5 K h

^{−1}(case D), −2 K h

^{−1}(case E), and −2.5 K h

^{−1}(case F). In addition to cases A–F, we also introduce two new cases with temporally varying surface cooling rates. This is a step toward parameterizing more realistic SABLs and certainly presents more challenges.

As shown in Fig. 1, the two new unsteady simulation cases last 10 h with the surface cooling rates following a hat function or a step function. The resulting surface temperature of the hat function decreases gradually with an inflection point in the middle, while that of the step function is piecewise linear with three segments. Before initiating stable conditions, we let the ABL develop to a quasi-equilibrium state under neutral stability for 2 h for LES and 5 h for SCM; we define quasi-equilibrium as the state where the height of the SABL and the surface fluxes of momentum and heat change relatively slowly with time. For convenience we denote the case with the hat function as “HAT” and the other one as “STEP.” The simulations of the six cases with steady forcing are run with 80^{3} nodes for the first 6 h and then the resulting outputs are interpolated to 162 × 162 × 160 nodes and run for another 4 h. The HAT and STEP cases are simulated with 162 × 162 × 160 nodes for the entire duration.

## 3. Applicability of the gradient-flux hypothesis

Before introducing the new model, we first examine the general applicability of the gradient-flux hypothesis using the LES and the GFDL SCM. The question we ask here is the following: given the values of *K _{m}* and

*K*that we can compute from the LES, can the SCM reproduce the vertical profiles of the mean quantities and turbulent fluxes in the LES using Eqs. (1)–(3)? This is the best-case scenario one can expect from the SCM as the eddy diffusivites are known perfectly. A major issue that can challenge the gradient-flux hypothesis and the ability of the SCM to reproduce the LES profiles is that

_{h}*K*used in Eqs. (1) and (2) is a scalar that takes the same value in both equations. This presupposes that the deflection angle between the gradient vector (∂

_{m}*U*/∂

*z*, ∂

*V*/∂

*z*) and the turbulent flux vector (

*τ*,

_{uw}*τ*) is zero—that is, the deflection angle . It is a specific implication of a more general assumption of the gradient-flux model: according to this, model fluxes at a given location are produced by local eddies of limited size and can thus be modeled using only local parameters. This assumption clearly fails in the convective ABL where coherent turbulent eddies contribute significantly to fluxes and are of the scale of the entire ABL (e.g., Kaimal and Finnigan 1994). However, in the SABL, the reduction in the vertical eddy size, which was observed in Huang and Bou-Zeid (2013) and other studies, suggests that the gradient-flux hypothesis might hold much better.

_{υw}In Fig. 2, the vertical profiles of *α* for cases A–F are plotted. Note that only the resolved parts of *τ _{uw}* and

*τ*are used to calculate

_{υw}*α*, since the angle between the SGS parts of

*τ*and

_{uw}*τ*and the velocity gradients is necessarily zero because of the use of a Smagorinsky-type SGS model. The deflection angle remains close to 0° in the lower SABL. However, in the upper SABL,

_{υw}*α*tends to be consistently negative with the minimum value around −20° at the SABL top for case A. This is probably associated with the increase of the Coriolis force and the resultant Ekman turning with height, and is of minor significance because the magnitudes of

*τ*and

_{uw}*τ*decrease with height and approach zero at the SABL top. By comparison, results for the neutral ABL and the convective ABL (not shown here) generally depict much larger angles between the stress and the strain (reaching up to 30° for the convective ABL).

_{υw}An important advantage of using LES is that the nonlinear dynamics of turbulence are largely resolved and can be used to deduce the temporally and spatially averaged emerging structure of the SABL. In this study, LES allows us to assess the applicability of the gradient-flux hypothesis by computing the a priori *K _{m}* and

*K*from the LES, imposing them in the SCM, and then comparing the resulting profiles of mean quantities and turbulent fluxes between the LES and the SCM. This is not a trivial test: we are requiring the simple linear eddy diffusivity model to reproduce the average behavior of the highly nonlinear turbulent fluxes resolved in LES. In connection with Eqs. (1)–(3), the LES-deduced

_{h}*K*and

_{m}*K*are calculated as follows:

_{h}Note that the mean quantities *U*, *V*, and Θ are the resolved parts, while the turbulent fluxes *τ _{uw}*,

*τ*, and

_{υw}*q*are the total fluxes that we obtain by adding the resolved and modeled SGS fluxes. This is needed since the SGS contribution to the fluxes is significant especially near the surface, while the direct SGS contribution to the means (i.e., the means of the SGS fields) is negligible. Both the mean quantities and turbulent fluxes are averaged horizontally and temporally over the last 3 h of the simulation. The vertical profiles of and are then imposed in the SCM to obtain the comparison between LES and SCM depicted in Fig. 3. A generally good match is achieved for both case A and case F, indicating that the gradient-flux hypothesis is generally valid for the SABL and could be used as a basis for SABL parameterization. There are nevertheless discrepancies between the LES and the SCM for the profile of the wind speed

*M*around the SABL top, which may be due to the relatively large magnitudes of

*α*in the same range shown in Fig. 2.

It is interesting to note that despite differences in the magnitude of wind speed between the LES and SCM around the SABL top, the two profiles are approximately parallel, yielding vertical gradients that are not significantly different. This in turn results in a good agreement of the stress profile between LES and SCM. The same applies to the temperature and heat transport; the agreement of the fluxes is better than the agreement of the mean profiles. Compared to case A, the agreement between the LES and the SCM for the profiles of wind speed and Θ for case F is less satisfactory. This may be caused by smaller fluxes and larger vertical gradients of wind speed and Θ in case F, which together result in smaller and less accurately determined diffusivities. Despite these discrepancies, the results presented in this section provide robust quantitative support to the application of the gradient-flux theory in the SABL.

## 4. A novel mixing-length model

The previous section validated the gradient-flux hypothesis for turbulence parameterization in the SABL using LES-derived eddy diffusivities. However, to be applicable in weather and climate models, a SABL parameterization must compute these diffusivities as functions of available fields. We proceed to assess currently available parameterization schemes for *K _{m}* and

*K*in terms of their accuracy as well as sensitivity to stabilities. The classical form of the

_{h}*l*proposed by Blackadar (1962) as

_{N}is used throughout this paper, where *l*_{∞} is a constant representing a turbulence mixing length far above the land surface. Typical values of *l*_{∞} used in large-scale atmospheric models range from 40 to 200 m (e.g., Beare et al. 2006; Cuxart et al. 2006; McCabe and Brown 2007). However, it is known that such a range (both the lower and the upper limits) is too high. In Cuxart et al. (2006), *l*_{∞} was reduced from 50 to 10 m for the Japan Meteorological Agency model and a closer fit with LES results was achieved. Kim and Mahrt (1992) derived *l*_{∞} = 14.2 m by analyzing aircraft data from three different field programs. A recent field study of the 1999 Cooperative Atmosphere–Surface Exchange Study (CASES-99) by Sun (2011) also revealed that *l*_{∞} is around 6 m above *z* = 20 m. Sorbjan (2012) obtained good agreement between the SCM and the LES using *l*_{∞} = 12 m based on the Surface Heat Budget of the Arctic Ocean (SHEBA) and the CASES-99 data. Our LES results indicate that a value of *l*_{∞} = 7 m is optimal to obtain a good match between the LES-computed and the one obtained from the traditional form in Eq. (6) under weak stabilities [Ri_{g} < 0.1, where the form in Eq. (6) is known to work rather well]. In Fig. 4, we compare the sharp form [Eq. (7)], the long-tails form [Eq. (8)], and the GFDL form [Eq. (11)] of *f _{m}* with those calculated from the LES for cases A–F; that is,

The figure illustrates that turbulent mixing remains significant, even when Ri_{g} is well beyond 0.2 in the LES; this contradicts the use of a critical Ri_{g} of 0.2 as in Eq. (6). Furthermore, does not appear to be a universal function of Ri_{g} only, since at a given Ri_{g} it changes between different cases. Also for case F, the change with Ri_{g} is not monotonic and one Ri_{g} can correspond to multiple values of . This undermines the use of any form of *f _{m}* as defined in Eq. (5). However, among the three parameterizations, the sharp form fits the LES results relatively better than the other two, while the long-tails form overestimates

*f*for the entire stability range of [0, 1] and the GFDL form overestimates

_{m}*f*for strong stabilities. One obvious reason for the failure of this form of the parameterization is that, as the stability increases, the dependence on the height above the surface does not decrease; the mixing length maintains the same dependence on

_{m}*z*, but is reduced as a function of Ri

_{g}. This does not agree with a large range of observations that suggest that the SABL dynamics over most of the SABL depth obey a

*z*-less similarity (Wyngaard 1973; Nieuwstadt 1984), where the distance to the surface is not relevant.

This local behavior is further confirmed by analyses of LES simulation data in the first part of this study (Huang and Bou-Zeid 2013), which suggests the existence of a *z*-independent mixing-length-scale *λ* that is inversely proportional to stability in the *z*-less part of the SABL (see Fig. 14 in Part I). Ri_{g} is selected as the stability parameter since it is fully determined by the mean variables solved for by the SCM and other coarse atmospheric models, so that the mixing-length scale in the *z*-less SABL is given by

In this study, we will determine *λ* empirically. Ideally, one would relate *λ* to some physical length scale. However, since good candidates for such a length scale (the Ozmidov scale or some length scale based on the turbulent kinetic energy for example) all depend on higher-order quantities not directly available in a SCM, such an approach would eliminate the succinctness of the parameterization that we wish to keep.

These findings prompt us to develop a new type of parameterization that does not rely on the framework of Eq. (5). We seek a formulation where the dependence on *z* is reduced as stability increases, and eliminated for very strong stabilities. A natural and simple way to combine the neutral limit of the mixing length [Eq. (14)] and the *z*-less stable limit [Eq. (16)] is to sum them inversely:

This formulation is similar to the one proposed by (Delage 1974), but instead of the local Obukhov length used by Delage in the third term on the right, we use Ri_{g} to keep the model explicit and to use a gradient-based stability measure rather than a flux based one since the gradient-based measures were found to perform better by Bou-Zeid et al. (2010). In the neutral limit where Ri_{g} = 0, Eq. (17) reduces to Eq. (14); in the stable limit where Ri_{g} becomes large such that the first two terms on the right side of Eq. (17) are negligible, this model yields the scaling sought in Eq. (16). In Fig. 5, the modeled mixing-length using Eq. (17) with parameters *l*_{∞} = 7 m and *λ* = 0.27 m (more details about the determination of these coefficients are presented in the next section) is plotted against computed as

The figure illustrates that this hybrid model generally works very well for cases B–F; for case A, the near-neutral case, the model underestimates the mixing length for the middle range of the SABL but captures the mixing length well near the surface and near the SABL top. For this near-neutral limit of case A, the impact of the last term we added is very minimal and the weakness of the model is related to the shortcomings of Eq. (14), which could be due to the value of *l*_{∞} = 7 m. As discussed above, there is currently no consensus on the value of *l*_{∞} under neutral conditions, which for example could depend on the depth of the ABL rather than be a constant [see more discussion on modeling of *l*_{∞} in Blackadar (1962)]. The failure of Eq. (14) could also be due to nonlocal effects related to the larger eddies in the neutral ABL. One can try to improve on the near-neutral formulations, but this is beyond the scope of this study since it is not critical in the stable ABL on which we focus here.

## 5. Evaluation with steady forcing

The proposed mixing-length model still needs to be calibrated (i.e., its empirical constants determined to match LES results) before it can be evaluated through the LES–SCM comparison. For convenience of notation, this new model will be called “HBG” (Huang, Bou-Zeid, and Golaz), and it will be evaluated versus LES results and also by comparing its performance to the sharp and GFDL forms of *f _{m}*. The long-tails form is dropped mainly because it artificially enhances mixing beyond what is suggested by the LES as shown in Fig. 4, and because we found its performance in the SCM to be poor. One could determine the optimal values of the model parameters

*l*

_{∞}and

*λ*using an a priori approach to match the mixing lengths and diffusivities produced by the LES. This is in fact what we have performed in the previous section, yielding

*l*

_{∞}= 7 m and

*λ*= 0.27 m. However, given that the diffusivities computed from the LES produce some errors when imposed in the SCM as illustrated in Fig. 3, an alternative a posteriori approach that determines the coefficients that best match the stress profile of the LES might be more appropriate to correct for the effect of the nonlinear dynamics of the Reynolds-averaged Navier–Stokes equations that are missing from the gradient-flux model. For this purpose, we use case A as the reference case for parameter determination for both the HBG model and the sharp model; thus we underline that both models are calibrated independently and benefit from the LES results. The results for the vertical stress profiles with different parameters are shown in Fig. 6 for the HBG and the sharp model. The results using the GFDL model with

*l*

_{∞}= 40 m and

*l*

_{∞}= 200 m, respectively, are also shown for reference and confirm that even

*l*

_{∞}= 40 m is too high to reasonably model turbulence in the SABL (note that the GFDL default value of

*l*

_{∞}is 500 m). While the a priori tuned parameters for the HBG model gives very reasonable results, a better agreement with LES is obtained when the model coefficients are calibrated a posteriori, yielding

*l*

_{∞}= 9 m and

*λ*= 0.45 m. The sharp model’s performance is also improved significantly with this calibration yielding

*l*

_{∞}= 9 m; this value is also found to improve the GFDL model results (not shown here) and hence will be used in the rest of the analysis for all three models. This is desirable since all models asymptotically approach the same form in the neutral limit.

The vertical profiles of wind speed and stress are depicted in Fig. 7 for cases A and F, with all models calibrated as discussed above using only case A. Since case A is the reference calibration case, all three models present a relatively good match between the SCM and the LES except the GFDL model, which slightly overestimates the stress in the upper SABL and underestimates the stress in the lower SABL. However, significant differences emerge between the HBG model and the other two when the stability is increased, as illustrated in the figure for the most stable simulation, case F. While the HBG model continues to yield a good LES–SCM agreement, the performance of the other two models deteriorates significantly. The sharp model largely underestimates stress throughout the SABL, with an error at the surface of about 25%, which is due to the underestimation of *f _{m}* for case F as shown in Fig. 4. This also causes a sharper vertical gradient of wind speed in the lower SABL. As for the GFDL model, the stress is also underestimated near the surface, but overestimated in the middle and upper SABL because in this model

*f*is forced to be much larger than indicated by the LES, as shown in Fig. 4 under high stabilities. This results in a spuriously deep SABL and a much weaker low-level jet (LLJ) around the SABL top with the GFDL model. In Fig. 8, we contrast the surface stress obtained by the three mixing-length models to that given by the LES. It is clear that the HBG model achieves a very good match with the LES for all six cases, while the sharp and GFDL models highly underestimate the surface stress for cases C–F. Overall, the contrast between the performance of the three models for cases A–F confirms the insensitivity of the HBG model performance to the strength of surface cooling, and the deterioration of the sharp and GFDL models’ performance as the strength of surface cooling is increased. This is despite the fact that all models were calibrated using results for case A.

_{m}All of the SCM results depicted thus far were obtained assuming *K _{h}* to be the same as

*K*(i.e., Prandtl number Pr =

_{m}*K*/

_{m}*K*= 1). However, our LES results indicate that Pr mostly falls in the range of [0.6, 0.9] within the SABL (see Huang and Bou-Zeid 2013, their Fig. 12). An overestimated Pr leads to a lower

_{h}*K*, which in turn gives rise to a lower heat flux and a sharper vertical gradient of mean temperature. This is shown in Fig. 9 for cases A and F where the profiles of mean temperature and heat flux obtained with different values of Pr are shown; cases B–E behave similarly. The figure illustrates that a value of Pr = 0.85 yields better results than Pr = 1, and will thus be used for the rest of this paper. A Pr varying with stability was also tested but the improvement in the results was insignificant. The deviations of the SCM results from those of the LES are reduced significantly when Pr = 0.85 is used, although the heat flux is now slightly underestimated in the upper SABL for case F.

_{h}## 6. Evaluation with unsteady forcing

In the previous section, we evaluated the HBG model using six cases with various steady surface cooling rates. It was shown that the new model is relatively insensitive to the change of surface cooling and the resulting change in stability, while the sharp and GFDL models do not yield good results for cases where the surface cools down rapidly. In this section, we further apply the three models to the two cases with the unsteady surface cooling rates illustrated in Fig. 1 (i.e., HAT and STEP). We use the same model coefficients determined in the previous section, where the models were all tuned similarly using LES data for near-neutral case A at quasi equilibrium. These unsteady cases present significant challenges to the SCM and to the turbulence parameterization since they feature surface forcing variability at a time scale comparable to or smaller than the time it takes the ABL to equilibrate to a new surface forcing. The real ABL under such unsteady forcing conditions will be continuously adapting to the forcing variability and will not be in quasi equilibrium. The LES can capture such unsteady effects since it resolves most of the turbulent scales and as such can capture the memory effect causing turbulent kinetic energy and other fields to have a lag in adjusting to the forcing variability. The SCM can capture some of these transient effects since it solves the momentum and sensible heat budgets, but it cannot capture the effects of unsteadiness on the turbulent kinetic energy since, by construction, the gradient-flux model assumes that the turbulence is in equilibrium with the mean fields. As such, we should not expect the models to produce results that match the corresponding LES results exactly, but rather a good LES–SCM match would consist of a relatively quick adjustment by the SCM following a forcing change.

We run the LES and the SCM separately under neutral conditions for a warmup period before hour 0. Note that we did not directly impose the mean profiles under neutral condition obtained from the LES in the SCM at hour 0, since the SCM would then need spinup time to adapt to these mean profiles. Consequently, the SCM and the LES do not exactly match at the start of the transient simulations but are rather started from the equilibrium neutral profiles each of them produces. Thus the criteria for comparing the different SCM closures would be how fast they can recover from the initial discrepancy and how fast they can adapt to subsequent changes in surface cooling.

In Fig. 10, we plot the time history of the wind speed at *z* = 10 m *M*_{10} and the friction velocity *u*_{*} for the HAT and STEP cases. As indicated by Fig. 7, increasing surface cooling tends to lower the height of the LLJ and consequently increase wind speed at the level below the LLJ. For the HAT case LES, this causes *M*_{10} generally to increase, after an adaptation period of about 2 h following the abrupt change from neutral conditions at time 0, until around hour 5 when the strongest surface cooling occurs, and then to decrease as the surface cooling rate decreases. This trend is well captured by the HBG model except during the initial adaptation period of the first several hours. However, the sharp and GFDL models start to significantly overestimate the wind speed after hour 4 by up to 0.4 m s^{−1} for the sharp model and 0.6 m s^{−1} for the GFDL model, which is consistent with Fig. 7 where it was shown that the sharp and GFDL models largely underestimate the stress and overestimate wind speeds in the lower SABL.

In the STEP case, the performance of the three models is comparable to their performance in the HAT case: the HBG model produces a good match of *M*_{10} with the LES, while the other two models significantly overestimate it starting at hour 2.5 where the surface cooling rate is suddenly increased. Note that around hour 2.5, both *M*_{10} and *u*_{*} present a sharp spike that reflects the development of the LLJ before the change in the surface cooling rate and the breakdown of this development because of this change. This feature is not well captured by any of the three models. As far as *u*_{*} is concerned, it generally decreases as the surface cooling rate increases at quasi equilibrium (see Part I of this study). This is to be expected since stronger thermal stability tends to suppress mechanical turbulence and the downward transport of momentum. However, since it takes time for the SABL to adapt to the change of surface forcing, this relationship is not clearly reflected in unsteady results depicted in Fig. 10. Nevertheless, the three SCM models generally produce the temporal trend of *u*_{*} in LES. However, the sharp model and the GFDL model tend to underestimate *u*_{*} and yield higher biases than the HBG model, which overestimates *u*_{*} after the initial adaption period for both unsteady cases.

In Fig. 11, we plot the temporal evolution of mean potential temperature at 10 m Θ_{10}, and surface heat flux *q _{s}*. The Θ

_{10}values produced by the HBG model and the LES are almost indistinguishable, while the GFDL model performs the worst with the maximum overestimation exceeding 1 K for both cases. For the HAT case, the magnitude of

*q*increases first and then decreases, similarly to the surface cooling rate. However, the strongest

_{s}*q*occurs with a delay of about 1 h after the surface cooling rate reaches its maximum. This is not observed for the STEP case since the surface forcing rate changes abruptly and affects

_{s}*q*immediately. All three models generally perform well in predicting

_{s}*q*, though the HBG model predicts the timing of the peak negative flux for the HAT more accurately than the other models. The underestimation of the amplitude of

_{s}*q*by the sharp model becomes noticeable after hour 5 for the two cases, which is again due to the sharp model’s weak mixing under high stabilities. Together with the results shown in Fig. 10, this discrepancy again underlines the variability of the performance of the sharp and GFDL models under different thermal stabilities, while the HBG models performance is superior and more consistent.

_{s}The time history of the SABL height *h* determined by two methods for the STEP case is shown in Fig. 12. Following Beare et al. (2006) and Cuxart et al. (2006), the first method defines *h* as the height where stress falls to 5% of its surface value divided by 0.95. The second method defines *h* as the lowest maximum of the wind speed, which is the LLJ when it exists (Melgarejo and Deardorff 1974). To distinguish the two methods, we denote the SABL height determined by the first method as *h*_{1} and that determined by the second as *h*_{2}. Over the 10 h of simulations, both *h*_{1} and *h*_{2} of the LES decrease from over 200 m for the initial neutral ABL to around 100 m. The *h*_{1} of the LES responds to the change of the surface cooling rate with a delay of about 1 h: it starts to decrease around hour 1 and again around hour 3.5. However, *h*_{1} of the three SCM models responds to the change of the surface cooling rate with a delay of about 4 h. In addition, *h*_{1} is highly overestimated by the GFDL model, while the other two models gradually converge to the LES after hour 4. The overestimation by the GFDL model is related to the overestimation of the turbulent diffusivities under high stabilities (Ri_{g} > 0.2 from Fig. 4), which usually occurs in the upper SABL (cf. Fig. 7). In the right panel of the figure, it is interesting to note the association of *h*_{2} with the formation of the LLJ. In the neutral ABL, the geostrophic speed at the ABL top is the maximum. As the SABL develops, the LLJ starts to form around hour 2 for the LES and hour 4 for the three SCM models, *h*_{2} decreases sharply afterward. Subsequently, the three SCM models converge to the LES results gradually.

In addition to the analysis of time histories, we also compare the vertical profiles of various parameters in Fig. 13 for the STEP case at hour 5. Both similarities and differences emerge as we compare Fig. 13 and the results for cases with steady forcing at quasi equilibrium in Fig. 7. Since the SCM does not respond to the change of surface forcing at the same rate as the LES, the LES–SCM match for the LLJ is less satisfactory in Fig. 13 than in Fig. 7. In Fig. 13 the LES wind speed peaks around *z* = 120 m, while it is around *z* = 170 m for the SCM. This also causes the HBG model to slightly underestimate *τ* and *q* compared to the LES. However, compared to the HBG model, the sharp model and the GFDL model performance is poorer, as was observed for quasi-equilibrium results in Fig. 7. The sharp model and the GFDL model generally underestimate the fluxes *τ* and *q* more significantly than the HBG model in the lower SABL, while the GFDL model largely overestimates the fluxes in the upper SABL.

## 7. Summary

In this paper, we conducted a suite of simulations using LES and the GFDL SCM to examine the effects of stability on the applicability of first-order parameterizations of the SABL using Ri_{g}. The simulations significantly extend the range of stabilities tested beyond previous studies tackling this problem. Unsteady conditions with variable surface forcing are also examined in detail.

We first investigated the applicability of the gradient-flux hypothesis in the SABL since it is the basis of the parameterizations we assess. We found that the deflection angle between the mean velocity gradient and the stress is generally close to 0° in the lower SABL, where turbulent transport is most important. We also imposed the turbulent diffusivities calculated from the LES in the SCM, and a good agreement was obtained for mean quantities and fluxes across varying stabilities. This suggests that the performance of first-order SABL parameterizations can be reasonably accurate if good models for turbulent diffusivities are formulated.

From the LES, we were able to obtain a priori stability correction functions *f _{m}*(Ri

_{g}) for cases with different stabilities and compare them to the forms suggested by currently available models. A fundamental problem all the models share is that

*f*is not found to be a universal function of Ri

_{m}_{g}. With this observation, we proposed a novel mixing-length model (the HBG model), which continues to use Ri

_{g}to characterize stability but does not use

*f*and does not invoke a critical Ri

_{m}_{g}. Instead, in the new model the mixing length is inversely proportional to Ri

_{g}under very high stabilities and converges to the model proposed by Blackadar (1962) under neutral conditions. The match between the mixing length predicted by this model and the one deduced from LES is noticeably better than the more widely used models, particularly under high stabilities.

The new model was then incorporated in the GFDL SCM and compared with other models using *f _{m}*(Ri

_{g}). Test cases with steady forcing as well as with unsteady forcing were used to examine its performance. The results illustrate that the performance of the traditional models using stability correction functions to reduce the neutral ABL mixing length is very sensitive to the change of stability. If we calibrate the parameters using cases under low stabilities, the performance of the traditional parameterization schemes break down under high stabilities. However, the new proposed model performs very well in all the steady-state test cases, even when Ri

_{g}increases up to around 1. For the unsteady surface forcing simulations, the new model continues to perform better than the other models we tested, but these conditions are significantly more challenging to closure formulations than steady-state cases.

## Acknowledgments

This work is supported by NSF Physical and Dynamic Meteorology Program under AGS-1026636 and by the Siebel Energy Challenge of Princeton University. The simulations were performed on the supercomputing clusters of the National Center for Atmospheric Research and of Princeton University.

### APPENDIX

#### The Effects of Δ*z* and Δ*t* on the Performance of the SCM

The relationship between the performance of the SCM and the imposed vertical-resolution Δ*z* and time step Δ*t* are explored here. We conduct two experiments with the SCM for case A using the HBG model with *l*_{∞} = 9 m, *λ* = 0.45 m, and Pr = 0.85. In the first one, Δ*z* varies from 2.5 to 50 m and Δ*t* is fixed at 180 s, which is small enough to maintain numerical stability for all three Δ*z* values. In the second experiment, Δ*t* ranges from 180 to 1200 s while Δ*z* is set as 2.5 m. The resulting vertical profiles of stress are shown in Fig. A1. The results of other variables share similar features and are not shown here. It is clear that Δ*z* does not have a significant effect on the performance of the SCM as long as the corresponding Δ*t* is small enough. However, the SCM does encounter numerical instability if Δ*t* is larger than a certain threshold, which corresponds to the minimum time of turbulent diffusion over Δ*z* as shown in the right panel of the figure.

To further quantify this threshold on Δ*t*, we run the SCM with a series of combinations of Δ*t* and Δ*z*, and calculate the root-mean-square error (RMSE) of the stress from the ground up to 250 m, using the LES profiles as the “true value.” In Fig. A2, the RMSE of stress is plotted against the dimensionless parameter *u*_{*}Δ*t*/Δ*z*. The figure illustrates that *u*_{*}Δ*t*/Δ*z* < 30 is required for the SCM to perform successfully, which is equivalent to Δ*t* < 30Δ*z*/*u*_{*}. Note that the *u*_{*}Δ*t*/Δ*z* parameter is equivalent to the more conventionally used *ν _{T}*Δ

*t*/Δ

*z*

^{2}, with the turbulent viscosity

*ν*=

_{T}*u*

_{*}Δ

*z*.

## REFERENCES

*The Atmospheric Boundary Layer.*Cambridge University Press, 316 pp.

*Atmospheric Boundary Layer Flows: Their Structure and Measurement.*Oxford University Press, 289 pp.

*Workshop on Micrometeorology,*D. A. Haugen, Ed., Amer. Meteor. Soc., 101–149.

## Footnotes

Current affiliation: CSIRO Marine and Atmospheric Research, Pye Laboratory, Canberra, Australian Capital Territory, Australia.