Implementation of a roughness sublayer parameterization in the Weather Research and Forecasting model (WRF version 3.7.1) and its evaluation for regional climate simulations

The roughness sublayer (RSL) is one compartment of the surface layer (SL) where turbulence deviates from Monin–Obukhov similarity theory. As the computing power increases, model grid sizes approach the gray zone of turbulence in the energy-containing range and the lowest model layer is located within the RSL. From this perspective, the RSL has an important implication in atmospheric modeling research. However, it has not been explicitly simulated in atmospheric mesoscale models. This study incorporates the RSL model proposed by Harman and Finnigan (2007, 2008) into the Jiménez et al. (2012) SL scheme. A high-resolution simulation performed with the Weather Research and Forecasting model (WRF) illustrates the impacts of the RSL parameterization on the wind, air temperature, and rainfall simulation in the atmospheric boundary layer. As the roughness parameters vary with the atmospheric stability and vegetative phenology in the RSL model, our RSL implementation reproduces the observed surface wind, particularly over tall canopies in the winter season by reducing the root mean square error (RMSE) from 3.1 to 1.8 m s⁻¹. Moreover, the improvement is relevant to air temperature (from 2.74 to 2.67 K of RMSE) and precipitation (from 140 to 135 mm per month of RMSE). Our findings suggest that the RSL must be properly considered both for better weather and climate simulations and for the application of wind energy and atmospheric dispersion.

1 Introduction

The planetary boundary layer (PBL) is important for the proper simulation of weather, climate, wind energy application, and air pollution. Turbulence plays a critical role in the spatiotemporal variation of the PBL structure through the turbulent exchanges of momentum, energy, and water between the atmosphere and Earth's surface. Because turbulent eddies in the PBL are smaller than the typical grid size in mesoscale and global models, their impacts must be properly parameterized for atmospheric models. The surface layer (SL) occupies the lowest 10 % of the PBL where the shear-driven turbulence is dominant. In the SL, Monin–Obukhov similarity theory (MOST), which is a zero-order turbulence closure, provides the relationships between the vertical distribution of wind and scalars and the corresponding fluxes in a given stability condition (Obukhov, 1946; Monin and Obukhov, 1954). The typical numerical weather prediction (NWP) and climate models are applied for the SL parameterization based on MOST to parameterize the subgrid-scale influences of the turbulent eddies in the PBL (e.g., Sellers et al., 1986, 1996).

The SL has two parts: the inertial sublayer (ISL) and the roughness sublayer (RSL). The ISL is the upper part of the SL, where MOST is valid and vertical variation of the turbulent fluxes is negligible. The RSL is the layer near and within the surface roughness elements (e.g., trees and buildings). The turbulent transport in the RSL has a mixing layer analogy, and the atmospheric flow depends on the roughness element properties (Raupach et al., 1996). Accordingly, the flux–gradient relationships in the RSL deviate from the MOST predictions, and the eddy diffusion coefficients are larger than the values in the ISL (e.g., Shaw et al., 1988; Kaimal and Finnigan, 1994; Brunet and Irvine, 2000; Finnigan, 2000; Hong et al., 2002; Dupont and Patton, 2012; Shapkalijevski et al., 2016; Zhan et al., 2016; Basu and Lacser, 2017).

Traditionally, the RSL has not been explicitly considered in global and mesoscale models because the PBL in the model is coarsely resolved, and the lowest model layer is well above the roughness elements accordingly. As the computing power increases, the grid size of the NWP model decreases and it gets close to the grid size of gray zone for the turbulence. Nevertheless, studies on the impact of a fine vertical resolution have not been relatively performed. From this perspective, the RSL has an important implication in atmospheric modeling research. The lowest model layer is typically approximately 30 m high, and its vertical resolution continues to be better; hence, the models have more than one vertical layer in the RSL, which extends to 2–3 times the canopy height. Furthermore, model outputs are sensitive to the selection of the lowest model level height (Shin et al., 2012), but its relation to the RSL has not yet been clearly investigated. Accordingly, turbulent transport in the RSL must be incorporated particularly in the mesoscale models if the vertical model levels are inside the RSL with an increase in the vertical model resolution.

The RSL function is a popular and simple method of incorporating the effects of the RSL in the observation and model (e.g., Raupach, 1992; Physick and Garratt, 1995; Wenzel et al., 1997; Mölder et al., 1999; Harman and Finnigan, 2007, 2008; de Ridder, 2010; Arnqvist and Bergström, 2015). The RSL function is defined as the observed relationship between the vertical gradient of wind and scalar variables and their corresponding fluxes in the RSL. Accordingly, simple relationships are appropriate for the land surface model in the climate model and for the mesoscale model (Physick and Garratt, 1995; Sellers et al., 1986, 1996). Despite the importance of the RSL, the Weather and Research Forecasting (WRF) model (Skamarock et al., 2008), which is one of the widely used models in the operation and research fields, does not consider the effects of the RSL for the regional weather and climate simulations. Harman and Finnigan (2007, 2008) and Harman (2012) (hereafter, HFs) recently proposed a relatively simpler RSL function that can be used in a wide range of atmospheric models. The RSL function of the HFs is based on a theoretical background and applicable to a wide range of atmospheric stabilities by succinctly satisfying the continuity of the vertical profiles of fluxes, wind, and scalars both at the top of the RSL and at the top of a canopy. The parameterization of HFs has recently been incorporated in a one-dimensional (1-D) PBL model and a land surface model (Harman, 2012; Shapkalijevski et al., 2017; Bonan et al., 2018).

Based on the abovementioned background, this study incorporates the RSL parameterization based on the RSL function of the HFs into the WRF model (version 3.7.1). For this purpose, we reformulate the HFs' RSL parameterization to implement it in the SL parameterization in the WRF model and then discuss the impacts of the RSL parameterization on the regional weather and climate simulations in terms of meteorological conditions near the Earth's surface. To the best of the authors' knowledge, our study is the first extensive attempt to incorporate the RSL parameterization into the WRF model and to validate it for regional climate simulations. Section 2 presents a brief discussion of the RSL parameterization of HFs and the implementation procedures into the WRF model. Section 3 explains the experimental and observational descriptions. Section 4 presents the impacts of the RSL parameterization. Section 5 ends the study with the concluding remarks.

2 RSL theory of the HFs

The roughness sublayer parameterization by HFs is adopted herein along with an explanation of the core of HFs' model, and the relevant details on this parameterization can be found in Harman and Finnigan (2007, 2008) and Harman (2012). Appendix A lists the symbols and variables used in this study in alphabetical order.

We first define the coordinate alignment for its application to the WRF. The revised MM5 SL scheme in the WRF model defines the vertical origin by the conventional zero-plane displacement height (d₀). The same coordinate system is also applied herein. The vertical coordinates z and $\stackrel{\mathrm{̃}}{z}$ in this coordinate system are defined as the distance from d₀ and from the terrain surface, respectively; therefore, their relation is $z=\stackrel{\mathrm{̃}}{z}-d_{\mathrm{0}}$ (Fig. 1). Note that a vertical origin in the HFs is at the canopy height (h). MOST says that a variable (C), such as wind speed (u) and temperature (T), has the following logarithmic vertical profile:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{k}{C_{*}}}\left(C\left(z\right)-C_{\mathrm{0}}\right)=\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z}{L}}\right)+{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right),\end{array}$$

where k is von Kármán constant; C_* is a C scale; C₀ is C at z₀; z₀ is the roughness length; ψ_c is the integrated similarity function of C; and L is the Obukhov length. The C profile based on the RSL function of the HFs is divided into two layers depending on the relative distance between the canopy height and the redefined zero-plane displacement height in the HFs ($d_{\mathrm{t}}=h-d_{\mathrm{0}}$): the upper-canopy layer (z>d_t), where the influence of additional mixing by the canopy exists, and the lower-canopy layer (z<d_t), where the canopy is the direct source and sink for drag and heat (Fig. 1). The vertical profile in the upper-canopy layer is described as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{k}{C_{*}}}\left(C\left(z\right)-C_{\mathrm{0}}\right)& =\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z}{L}}\right)+{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right)\\ & +\underset{z}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{c}}\left(\mathrm{1}-{\widehat{\mathit{\varphi }}}_{\mathrm{c}}\right)}{z^{\prime }}}\mathrm{d}{z}^{\prime },\end{array}\end{array}$$

where ϕ_c is the similarity function of C, and ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ is an RSL function of C. In the z→∞ limit, ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ is equal to 1 and the wind speed converges to the MOST prediction. The last term on the right-hand side represents the additional mixing caused by the roughness element due to the coherent canopy turbulence and can be replaced by ${\widehat{\mathit{\psi }}}_{\mathrm{c}}$, which is an integrated RSL function of C. The vertical profile from the HFs for the RSL deviates from that of MOST because of ${\widehat{\mathit{\psi }}}_{\mathrm{c}}$, thereby adjusting the logarithmic profile. The ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ is introduced as follows:

$$\begin{array}{}\text{(3)}& {\widehat{\mathit{\varphi }}}_{\mathrm{c}}=\mathrm{1}-c_{\mathrm{1}}\exp \left[-c_{\mathrm{2}}{\displaystyle \frac{\mathit{\beta }}{l_{\mathrm{m}}}}z\right].\end{array}$$

c₁ and c₂ are then determined from the continuity of ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ at the canopy top, and the RSL function, ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$, exponentially converges to zero above the RSL. In the lower canopy layer, C has the following exponential form:

$$\begin{array}{}\text{(4)}& C\left(z\right)-C_{\mathrm{0}}=\left(C_{\mathrm{h}}-C_{\mathrm{0}}\right)\exp \left(f{\displaystyle \frac{z-d_{\mathrm{t}}}{\mathrm{2}{d}_{\mathrm{t}}}}\right).\end{array}$$

The RSL functions vary with atmospheric stability through β:

$$\begin{array}{}\text{(5)}& \mathit{\beta }=\left\{\begin{array}{lc}\frac{{\mathit{\beta }}_{\mathrm{N}}}{{\mathit{\varphi }}_{\mathrm{m}}\left(z=d_{\mathrm{t}}\right)}& \frac{L_{\mathrm{c}}}{L}>-\mathrm{0.15}\\ \frac{k}{\mathrm{2}{\mathit{\varphi }}_{\mathrm{m}}\left(z=d_{\mathrm{t}}\right)}+\frac{\frac{{\mathit{\beta }}_{\mathrm{N}}}{{\mathit{\varphi }}_{\mathrm{m}}\left(z=d_{\mathrm{t}}\right)}-\frac{k}{\mathrm{2}{\mathit{\varphi }}_{\mathrm{m}}\left(z=d_{\mathrm{t}}\right)}}{\mathrm{1}+\mathrm{2}{\left|\frac{L_{\mathrm{c}}}{L}+\mathrm{0.15}\right|}^{\mathrm{1.5}}}& \frac{L_{\mathrm{c}}}{L}\le -\mathrm{0.15}\end{array}\right.,\end{array}$$

where L_c is a canopy penetration depth defined as

$$\begin{array}{}\text{(6)}& L_{\mathrm{c}}={\left(c_{\mathrm{d}}a\right)}^{-\mathrm{1}}={\displaystyle \frac{\mathrm{4}h}{\mathrm{LAI}}},\end{array}$$

where c_d is a drag coefficient at the leaf scale and a is the leaf area density. The parameters d_t and z₀ also depend on the stability because of their dependence on β:

$$\begin{array}{}\text{(7)}& {\displaystyle }{d}_{\mathrm{t}}=h-d_{\mathrm{0}}={\displaystyle \frac{l_{\mathrm{m}}}{\mathrm{2}\mathit{\beta }}}={\mathit{\beta }}^{\mathrm{2}}{L}_{\mathrm{c}},\text{(8)}& {\displaystyle }\begin{array}{rl}{z}_{\mathrm{0}}=& d_{\mathrm{t}}\exp \left[-{\displaystyle \frac{k}{\mathit{\beta }}}\right]\exp \left[-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{d_{\mathrm{t}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right)\right]\\ & \exp \left[\underset{d_{\mathrm{t}}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{m}}\left(\mathrm{1}-{\widehat{\mathit{\varphi }}}_{\mathrm{m}}\right)}{z^{\prime }}}\mathrm{d}{z}^{\prime }\right].\end{array}\end{array}$$

Figure 1 The schematic diagram to describe the vertical coordinate systems used in this study. (a) The vertical coordinate of the Yonsei University surface layer (YSL) schemes (z) is defined as the distance from the conventional zero-plane displacement height (d₀) and (b) the convectional coordinate ($\stackrel{\mathrm{̃}}{z}$) is defined as the distance from terrain surface. Here, z₀, d_t, h, and z_r are roughness length, the redefined zero-plane displacement height in Harman and Finnigan (2007), canopy height, and the lowest model layer in the conventional coordinate, respectively.

3 Incorporation of the roughness sublayer parameterization into the WRF model

The RSL parameterization of the HFs described above is implemented in the Jiménez et al. (2012) revised MM5 surface layer scheme and Noah land surface model in the WRF (hereafter called the Yonsei University surface layer (YSL) scheme) because of theoretical consistency between the HFs and PBL parameterization. To incorporate the RSL parameterization, it is necessary to modify the SL scheme as follows (Fig. 2): the first step is to compute the bulk Richardson number at the lowest model layer, Bi_b, by the original equation of Jiménez et al. (2012) (Eq. 9 in their study):

$$\begin{array}{}\text{(9)}& B{i}_{\mathrm{b}}={\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{a}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{va}}-{\mathit{\theta }}_{\mathrm{vg}}}{{\left[u\left(z_{\mathrm{r}}\right)\right]}^{\mathrm{2}}}}z.\end{array}$$

The second step is to iteratively calculate the atmospheric stability (z_r∕L) as follows with an accuracy of 0.01:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{z_{\mathrm{r}}}{L}}=B{i}_{\mathrm{b}}{\displaystyle \frac{{\left[\ln \left(\frac{z_{\mathrm{r}}}{z_{\mathrm{0}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{r}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}}}{L}\right)+{\widehat{\mathit{\psi }}}_{\mathrm{m}}\right]}^{\mathrm{2}}}{\left[\ln \left(\frac{\mathit{\rho }{c}_{\mathrm{p}}k{u}_{*}^{n-\mathrm{1}}{z}_{\mathrm{r}}}{c_{\mathrm{s}}}+\frac{z_{\mathrm{r}}}{z_l}\right)-{\mathit{\psi }}_{\mathrm{h}}\left(\frac{z_{\mathrm{r}}}{L}\right)+{\mathit{\psi }}_{\mathrm{h}}\left(\frac{z_l}{L}\right)+{\widehat{\mathit{\psi }}}_{\mathrm{h}}\right]}},\end{array}$$

where $u_{*}^{n-\mathrm{1}}$ is the friction velocity (u_*) at the previous time step. Equation (10) is different from Eq. (23) of Jiménez et al. (2012) by the RSL functions (i.e., ${\widehat{\mathit{\psi }}}_{\mathrm{m}}$ and ${\widehat{\mathit{\psi }}}_{\mathrm{h}}$). After z_r∕L is determined, the third step is to iteratively update d_t and β using Eqs. (5) and (7) with an accuracy of 0.0001 because they are intercorrelated with each other. Subsequently, z₀ is iteratively achieved with an accuracy of 0.0001 using Eq. (8) at the given z_r∕L, β, and d_t. The u profile is determined using Eqs. (2) and (4). Following Jiménez et al. (2012), the profile of a scalar, such as T, is determined by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{k}{C_{*}}}\left(C\left(z\right)-C_{\mathrm{0}}\right)& =\ln \left({\displaystyle \frac{\mathit{\rho }{c}_{\mathrm{p}}k{u}_{*}^{n-\mathrm{1}}{z}_{\mathrm{r}}}{c_{\mathrm{s}}}}+{\displaystyle \frac{z_{\mathrm{r}}}{z_l}}\right)-{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z}{L}}\right)\\ & +{\mathit{\psi }}_{\mathrm{c}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right)+\underset{z}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{c}}\left(\mathrm{1}-{\widehat{\mathit{\varphi }}}_{\mathrm{c}}\right)}{z^{\prime }}}\mathrm{d}{z}^{\prime },\end{array}\end{array}$$

for the upper-canopy layer. Equation (4) is used for the lower-canopy layer. Finally, u_* is given to

$$\begin{array}{}\text{(12)}& u_{*}={\displaystyle \frac{ku\left(z_{\mathrm{r}}\right)}{\left[\ln \left(\frac{z_{\mathrm{r}}}{z_{\mathrm{0}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{r}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}}}{L}\right)+{\widehat{\mathit{\psi }}}_{\mathrm{m}}\right]}},\end{array}$$

and the aerodynamic conductance from z=0 to z=z_r(g_a) in the RSL is given to

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{1}/g_{\mathrm{a}}=\underset{\mathrm{0}}{\overset{z_{\mathrm{r}}}{\int }}{\displaystyle \frac{\mathrm{d}z}{K_{\mathrm{h}}\left(z\right)}}\equiv \mathrm{1}/g_{{\mathrm{a}}_{\mathrm{1}}}+\mathrm{1}/g_{{\mathrm{a}}_{\mathrm{2}}}+\mathrm{1}/g_{{\mathrm{a}}_{\mathrm{3}}}\text{(14)}& {\displaystyle }\begin{array}{rl}\mathrm{1}/g_{{\mathrm{a}}_{\mathrm{1}}}& \equiv \underset{\mathrm{0}}{\overset{z_l}{\int }}{\displaystyle \frac{\mathrm{d}z}{K_{\mathrm{h}}\left(z\right)}}=\underset{\mathrm{0}}{\overset{z_l}{\int }}{\displaystyle \frac{\mathrm{d}z}{k{u}_{*}z+\frac{c_{\mathrm{s}}}{\mathit{\rho }{c}_{\mathrm{p}}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{k{u}_{*}}}\ln \left[{\displaystyle \frac{k{u}_{*}{z}_l}{c_{\mathrm{s}}/\mathit{\rho }{c}_{\mathrm{p}}}}+\mathrm{1}\right]\end{array}\text{(15)}& {\displaystyle }\mathrm{1}/g_{{\mathrm{a}}_{\mathrm{2}}}=\underset{z_l}{\overset{d_{\mathrm{t}}}{\int }}{\displaystyle \frac{\mathrm{d}z}{K_{\mathrm{h}}\left(z\right)}}={\displaystyle \frac{S_{\mathrm{c}}}{{\mathit{\beta }}^{\mathrm{2}}{u}_{\mathrm{h}}}}\left[\exp \left(\mathit{\beta }{\displaystyle \frac{d_{\mathrm{t}}-z_l}{l_{\mathrm{m}}}}\right)-\mathrm{1}\right],\text{(16)}& {\displaystyle }\begin{array}{rl}\mathrm{1}/g_{{\mathrm{a}}_{\mathrm{3}}}& =\underset{d_{\mathrm{t}}}{\overset{z_{\mathrm{r}}}{\int }}{\displaystyle \frac{\mathrm{d}z}{K_{\mathrm{h}}\left(z\right)}}=\underset{d_{\mathrm{t}}}{\overset{z_{\mathrm{r}}}{\int }}{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{h}}}{k{u}_{*}z}}\mathrm{d}z={\displaystyle \frac{\mathrm{1}}{k{u}_{*}}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{r}}}{d_{\mathrm{t}}}}\right)\right.\\ & \left.-{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{r}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{d_{\mathrm{t}}}{L}}\right)+{\widehat{\mathit{\psi }}}_{\mathrm{h}}\left(\mathit{\beta }{\displaystyle \frac{z_{\mathrm{r}}}{l_{\mathrm{m}}}}\right)-{\widehat{\mathit{\psi }}}_{\mathrm{h}}\left(\mathit{\beta }{\displaystyle \frac{d_{\mathrm{t}}}{l_{\mathrm{m}}}}\right)\right].\end{array}\end{array}$$

Computing time increased only 8 % with 61 vertical layers in the YSL scheme despite more iterations in the YSL scheme compared to the revised MM5 SL scheme.

Figure 2 Flow diagram of the RSL parameterization. The gray boxes indicate the iteration module.

4 Numerical experimental design

This study evaluated the YSL scheme by making a 1-D offline and real case simulations. The 1-D offline simulations were done to test the YSL scheme performance without feedback with the atmosphere. The 1-D offline simulation consisted of the YSL and the revised MM5 SL schemes coupled with the Noah land surface model (i.e., offCTL and offRSL experiments) (Table 1). In the offline simulations, offCTL and offRSL indicate numerical experiments with and without the RSL parameterization, respectively, and were driven by the idealized forcing data described in Table 1. The real case simulation consisted of two experiments: reproducing 1 month of winter (January 2016) with the revised MM5 SL scheme and the YSL scheme (hereafter referred to as the rCTL and rRSL experiments). The rCTL and the rRSL employed the same physics package, except for the SL scheme and the land surface model (Lee and Hong, 2016 and references therein). One-way nesting was applied herein in a single-nested domain with a Lambert conformal map projection to east Asia (Fig. 3). A 9 km horizontal resolution domain 2 was then embedded in the 27 km resolution domain 1 with 61 vertical layers. The initial and boundary conditions were produced using the National Center for Environmental Prediction Final Analysis data (1^∘ × 1^∘).

Table 1 Idealized boundary condition for the one-dimensional offline simulation. Constant z₀ is used only in the CTL experiment.

Figure 3 Domains and land-use category (USGS) of the real case simulation. Black circles denote the automatic synoptic observing system in South Korea used for the model evaluation.

5 Observation data for the model evaluation

The model performance was examined against the surface wind speed, surface temperature and precipitation observed at 46 Automated Synoptic Observing System (ASOS) sites in South Korea (Fig. 3). Quality control of the data includes gap detection, a limit test, and a step test based on the standard of the World Meteorological Administration and Korea Meteorological Administration (KMA) (Zahumenský, 2004; Hong et al., 2019). For the model evaluation of the real case simulation, the three different measures of the correlation coefficients, centered root mean square differences (RSMDs), and standard deviations of the model (σ_m) normalized by that of the observation (σ_o) are together shown in a Taylor diagram (Taylor, 2001). In the Taylor diagram, a point nearer the observation at a reference point (OBS) can be considered to give a better agreement with the observation. We also provide the root mean square error (RMSE) and the mean bias (MB) with the pattern correlation for the rainfall simulation evaluation.

6 Results

6.1 Offline simulations

Figure 4 shows the roughness parameters (i.e., z₀, d_t, and β) as a function of the normalized atmospheric stability (L_c∕L) from the 1-D offline simulation of the YSL scheme (offRSL). The offline YSL simulations reproduced the results of HFs. The roughness parameters varied with the atmospheric stability, L_c∕L, and had peaks at weakly unstable conditions. These dependencies of the roughness parameters on the atmospheric stability are distinct from typical manner of dealing with the roughness parameters as a constant in atmospheric models. The roughness length is indeed constant based on the land cover in all the SL schemes in the WRF.

Figure 4 Roughness length (a), redefined displacement height (b) and β (c) at a given normalized stability (L_c∕L) from the offRSL simulation with the YSL scheme. Roughness length and redefined displacement height are normalized by their values in a neutral condition (z_0N and d_tN), respectively.

Figure 5 indicates that the impacts of the RSL are also a function of L_c, which is a function of LAI and h (Eq. 6), thus leading to both diurnal and seasonal variation of canopy roughness. Consequently, the roughness parameters showed daily and seasonal variations. Overall, the roughness length in the YSL was larger than that in the revised MM5 SL scheme, particularly in a smaller z∕L (i.e., neutral and unstable conditions) and a larger L_c (i.e., small LAI and/or large h). The roughness length in a stable condition showed relatively smaller changes with z∕L and L_c compared to those in the unstable condition. Our findings suggest that a small LAI in the winter season makes a larger L_c because of the smaller LAI, thereby leading to relatively larger differences of z₀ between the YSL scheme and the default WRF scheme. On the contrary, a similar value of z₀ was observed in summer because of the larger LAI. Note that the revised MM5 SL scheme does not consider d_t and β.

Figure 5 Roughness length difference (m) between offCTL and offRSL (offRSL–offCTL) at given atmospheric stability (z∕L) and penetration depth (L_c).

The RSL function, ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$, was introduced to consider the additional mixing caused by the roughness element. Accordingly, ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ should asymptotically converge to the MOST profile (i.e., ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}\to \mathrm{1}$) as z increases with the continuous vertical profiles of the wind and the temperature. The YSL scheme reproduced these properties of ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$ and matched with the observed profiles inside canopies: the YSL scheme showed exponential profiles under the canopy top and logarithmic profiles above the canopy top (Fig. 6). The wind speed and the air temperature above the canopy top were smaller than predicted by MOST because ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}<\mathrm{1}$ in the offRSL experiments. Furthermore, the YSL scheme produced wind and temperature within the canopy (i.e., $\stackrel{\mathrm{̃}}{z}<z_{\mathrm{0}}+d_{\mathrm{0}}$), thereby providing additional useful information on the atmospheric dispersion inside the canopy.

Figure 6 (a) Profiles of the RSL function for momentum (${\widehat{\mathit{\varphi }}}_{\mathrm{m}}$, solid line) and heat (${\widehat{\mathit{\varphi }}}_{\mathrm{h}}$, dashed line), (b) wind speed (m s⁻¹), and (c) temperature (K) at a neutral condition from offCTL (black) and offRSL (gray). The height of conventional coordinate system ($\stackrel{\mathrm{̃}}{z}$) is normalized by the canopy height (h).

The roughness length changes in the YSL scheme eventually produced changes in the surface energy balance with the atmospheric stability (Fig. 7). In the 1-D offline simulations based on the conditions in Table 1, the YSL scheme produced a larger z₀ in the unstable and near-neutral conditions but a smaller z₀ in $z/L>\mathrm{3}$ compared to the offCTL. The aerodynamic conductance (g_a) in the YSL scheme was larger in all the stability conditions, even in the stable conditions in which the YSL provided a smaller z₀ because the additional term in Eq. (13), $g_{{\mathrm{a}}_{\mathrm{2}}}$, dominated over the other effects in the g_a calculation. Accordingly, H and λE in the YSL scheme were larger than those in the revised MM5 SL scheme. Our finding implies stronger fluxes from the YSL scheme when the gradient of quantity is the same. However, the impact of the increased g_a was asymmetrical in H and λE depending on the soil moisture content. In this case simulation, an increase in λE was dominant because the wet condition made more partitioning of the available energy into the latent heat flux first in the model. However, in the dry condition (i.e., less soil water content), the YSL produced a larger H without a substantial increase of λE (Fig. S1 in the Supplement). A significant increase in λE was found along with a decrease in H in the strong unstable conditions (Fig. 7) because of the wet soil moisture of 0.25 m³ m⁻³ in the offRSL simulation in Table 1. The slight increase in the net radiation was mainly associated with the reduced outgoing longwave radiation caused by the smaller surface temperature in the offRSL.

Figure 7 (a) Roughness length (m), (b) aerodynamic conductance (m s⁻¹), (c) sensible heat flux (W m⁻²), (d) latent heat flux (W m⁻²), and (e) net radiation (W m⁻²) at a given atmospheric stability (z∕L). The black lines denote offCTL, while the gray lines denote offRSL.

6.2 Real case simulations

Figure 8 shows the real case simulation of the roughness length, 10 m wind speed (u₁₀), and 2 m air temperature (T₂). We discuss herein the real cases in the winter season because of stronger effect of the roughness sublayer. The results for the summer season can be found in the Supplement. The roughness length in the rCTL experiment was prescribed from the vegetation data table (i.e., VEGPARM table in the WRF model) and modified by the vegetation fraction (Fig. 8a).

Figure 8 (a) Roughness length (m), (b) 10 m wind speed (m s⁻¹), and (c) daytime 2 m temperature (K) of the (left) rCTL experiment and (right) the difference (rRSL–rCTL). The results are averaged over a period of 1 month and masked out over the ocean.

Overall, the YSL scheme (rRSL experiment) produced 0.2–2.0 m larger z₀ than the default values in the rCTL experiment over the tall canopies, where L_c was large. In contrast, the YSL scheme produced a similar or even slightly smaller z₀ over the short canopies compared to the rCTL experiment. Importantly, the changes of z₀ made direct impacts on the momentum fluxes and thus surface wind speed (Fig. 8b). The typical u₁₀ in the rCTL was larger than approximately 3 m s⁻¹, and a much stronger wind (>6 m s⁻¹) was observed along the mountains, making a positive bias against the observation. Overall, the YSL scheme reproduced the better observed diurnal variation by reducing the positive bias of the wind speed (Table 2, Fig. 9). Over the tall forest canopies, u₁₀ in the rRSL was reduced by approximately 30 %; however, the region of the increased wind speed corresponded to the short canopies, where the roughness length increased (Fig. 8a and b). The YSL scheme particularly provided a better RMSD and correlation coefficient but less diurnal variability of wind speed because of a relatively larger reduction of the daytime wind speed (Fig. 9). MB and RMSE decreased from 2.4 to 1.0 m ⁻¹ and from 3.1 and 1.8 m s⁻¹. The Taylor diagram shows that the overall performance of the YSL scheme is better than the default WRF simulation at all the 46 sites. In the Taylor diagram, the statistics moved toward the observation, except for one site, indicating an overall improvement of 2 m air temperature in the YSL scheme; however, the impact of the RSL was not as large as the wind speed (Table 2, Fig. 10).

Figure 9 (a) A 1-month mean diurnal variation of 10 m wind speed and (b) the Taylor diagram showing the correlation coefficient, normalized centered root mean square differences (RMSDs), and standard deviations of the models (σ_m) normalized by that of observation (σ_o) from observation (gray), rCTL experiment (black), and rRSL experiment (red). The vectors indicate the changes of the statistics from rCTL to rRSL. The arrows indicate those from rCTL to rRSL. Every vector shows the movement toward the observation, thereby suggesting the model improvement.

Figure 10 Same as in Fig. 9 but for the 2 m temperature.

Table 2 Statistics of the 10 m wind speed, 2 m temperature, and rain rate. The top statistics are presented in bold.

Similar to the increases of the aerodynamic conductance in the offline simulations, the YSL scheme in the real case simulation (i.e., the rRSL simulation) simulated a larger g_a, particularly in the forest canopies and mountain regions (Fig. 11a). This larger g_a in the YSL scheme led to the increases of the latent heat fluxes by approximately 20 W m⁻², with an eventual reduction of the soil water content (Fig. 12a). The sensible heat fluxes in the rCTL experiments were generally approximately 80 W m⁻², except over the snow-covered region where H was approximately 40 W m⁻². As described in the offline simulation, the changing sign of H in the rRSL depended on the soil moisture content because evapotranspiration is limited in dry soils at given available energy (Figs. 11b and 12b). Consequently, the available energy ($=H+\mathit{\lambda }E$) increased in the YSL scheme, and a larger λE in the rRSL led to a cooler temperature than that in the rCTL experiment (Fig. 8c).

Figure 11 (a) Aerodynamic conductance (m s⁻¹), (b) daytime sensible heat flux (W m⁻²), and (c) daytime latent heat flux (W m⁻²) of the (left) rCTL experiment and (right) the difference (rRSL–rCTL). The results are averaged over a period of 1 month and masked out over the ocean.

Figure 12 (a) Difference of the soil moisture (m³ m⁻³) (rRSL–rCTL) and (b) snow cover (%) of rCTL. The results are averaged over a period of 1 month and masked out over the ocean.

During the winter simulation period, precipitation was observed over an extensive area in the domain, and snow was dominant over the northeastern side of the domain (Figs. 12b and 13). The overall total precipitation in the YSL scheme increased, and the skill score indicated a better simulation of the total amount of precipitation (Table 2, Fig. 13). The pattern correlation of precipitation also increased from 0.972 to 0.978 in the YSL scheme based on 656 rain gauge stations, indicating a better match of the precipitation bands. Despite the increase in λE, precipitation decreased in several regions (Figs. 11b and 13b). The differences were not significant in the summer season, and the skill scores in the YSL scheme were similar to the default WRF simulation because our implemented RSL parameterization started to converge to the default WRF in a smaller L_c (i.e., larger LAI and/or smaller h) and strong synoptic influences by the summer heavy rainy period (Table S1 in the Supplement, Figs. S2–S6).

Figure 13 (a) The 1-month accumulated precipitation of the rCTL experiment (mm) and (b) difference (rRSL–rCTL). (c) Taylor diagram showing the correlation coefficient, normalized centered RMSD, and the standard deviations of models (σ_m) normalized by that of the observation (σ_o) and from the rain rate (mm h⁻¹) of the rCTL experiment (black) and the rRSL experiment (red) during 1 month at 656 rain gauges.

7 Summary and concluding remarks

Turbulent fluxes regulate the planetary boundary layer; thus, they are a crucial process for weather, climate, and air pollution simulations. Most of the NWP and climate models are commonly applied for MOST to compute the turbulent fluxes near the Earth's surface. MOST can be only applicable in the inertial layer and turbulence deviates from MOST in the roughness sublayer. Importantly, the roughness sublayer, the important compartment of the SL, has not been properly parameterized in the model. Increasing the computing power enables us to use more vertical layers in the atmospheric models. Accordingly, the RSL must be incorporated into the model properly to simulate the atmospheric processes in the gray zone. This study proposed the YSL scheme, which incorporated the RSL into the WRF model, based on the RSL model proposed by Harman and Finnigan (2007, 2008) and Harman (2012). We also investigated the impacts of the RSL parameterization on the weather and climate simulations. For these purposes, we designed a series of offline simulations with an idealized boundary condition and real case simulations to evaluate the performance of the YSL scheme against the observation data.

The 1-D offline simulation revealed that the YSL scheme successfully reproduced the features reported in various canopies. The RSL function, ${\widehat{\mathit{\varphi }}}_{\mathrm{c}}$, asymptotically increased to 1, and the vertical gradients of the wind speed and the temperature decreased in the RSL as z increased, thereby converging to the MOST prediction. Notably, unlike the typical assignment of the roughness parameters (i.e., z₀, d_t, and β) as a constant, the roughness parameters are functions of the atmospheric stability (z∕L) and L_c. The roughness parameters had a maximum in the weakly unstable condition and in larger L_c (i.e., large h or small LAI). In most conditions, the YSL scheme provided a larger roughness length, thereby producing a wind speed slower than that of the revised MM5 SL scheme. The YSL scheme simulated a colder surface temperature in the unstable conditions.

Meanwhile, the real case simulation showed that the RSL-incorporated WRF produced a larger z₀ than the default WRF. This increase in z₀ and its change with atmospheric stability eventually made substantial impacts on wind, and temperature near the surface, momentum transfer, surface energy balance, and precipitation. First, an increase of z₀ produced larger momentum fluxes and a smaller 10 m wind speed when the YSL scheme was applied, leading to the mitigation of positive bias in the wind speed in the revised MM5 SL scheme. The larger z₀ also made increases in the available energy. This increased available energy is related to the surface cooling caused by the increases in the latent heat fluxes in the wet surface conditions when the RSL parameterization is applied. As a result, these changes in the climate near the surface and the surface energy balance resulted in more precipitation, thereby giving a better simulation of the amount of precipitation and its spatial pattern.

Our results indicate that the RSL parameterization can be a promising option for resolving the typical overestimation of the surface wind speed of the WRF model, particularly in the tall vegetation and low LAI, with slight increase of computing time (e.g., Hu et al., 2010, 2013; Shimada and Ohsawa, 2011; Shimada et al., 2011; Wyszogrodzki et al., 2013; Lee and Hong, 2016). The improvement caused by the RSL parameterization is useful in air quality modeling and wind energy estimation by better weather and climate in the planetary boundary layer. A further study is necessary to evaluate the characteristics of the YSL scheme in various cases particularly at gray-zone resolutions.

Appendix A: List of symbols and definitions

Symbols	Definitions
a	Leaf area density
Bib	Bulk Richardson number, at the lowest model layer
cd	Drag coefficient at the leaf level
cp	Specific heat for air
cs	Effective heat transfer coefficient for nonturbulent processes (Carlson and Boland, 1978; Jiménez et al., 2012)
C	Variable at z, such as u and T
C0	C at z=z0
Ch	C at h
C*	Scale of C
d0	Conventionally defined zero-plane displacement height
dt	Redefined zero-plane displacement height in Harman and Finnigan (2007)
f	Parameter related the depth scale of the scalar profile ($\equiv \frac{\mathrm{1}}{\mathrm{2}}(\sqrt{\mathrm{1}+\mathrm{4}{r}_{\mathrm{c}}{S}_{\mathrm{c}}}-\mathrm{1}))$
g	Gravitational acceleration
ga	Aerodynamic conductance
h	Canopy height
k	von Kármán constant
lm	Mixing length for momentum
L	Obukhov length
Lc	Canopy penetration depth ($\equiv (c_{\mathrm{d}}a{)}^{-\mathrm{1}}=\frac{\mathrm{4}h}{\mathrm{LAI}}$)
LAI	Leaf area index
p	Pressure at z
q	Water vapor mixing ratio at z
rc	Canopy Stanton number
SW	Downward shortwave radiation
Sc	Turbulent Schmidt number at canopy top
Sm	Soil moisture
T	Air temperature at z
T2	Air temperature at 2 m
Tsk	Skin temperature
u	Wind speed at z
u10	Wind speed at 10 m
uh	Wind speed at h
u*	Friction velocity
$u_{*}^{n-\mathrm{1}}$	Previous time step value of u*
z	Height from d0
$\stackrel{\mathrm{̃}}{z}$	Height from terrain surface
z0	Roughness length
zl	Viscous sublayer depth = 0.001 (Carlson and Boland, 1978; Jiménez et al., 2012)
zr	Height of the lowest model layer
zr∕L	Atmospheric stability
β	$u_{*}/u_{\mathrm{h}}$
βN	β at neutral condition (=0.374)
θa	Potential temperature of the air at zr
θva	Virtual potential temperature of the air at zr
θvg	Virtual potential temperature of the air at ground
ρ	Density of air
ϕC	Similarity function of C
${\widehat{\mathit{\varphi }}}_{\mathrm{C}}$	RSL function of C
ψC	Integrated similarity function of C
ψh	Integrated similarity function of heat
ψm	Integrated similarity function of momentum
${\widehat{\mathit{\psi }}}_{\mathrm{C}}$	Integrated RSL function of C
${\widehat{\mathit{\psi }}}_{\mathrm{h}}$	Integrated RSL function of heat
${\widehat{\mathit{\psi }}}_{\mathrm{m}}$	Integrated RSL function of momentum