Implementation of a roughness sublayer parameterization in the Weather Research and Forecasting model (WRF version 3.7.1) and its validationevaluation 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 to the gray zone of turbulence in the 10 energy containing range and the lowest model layer is located within the RSL. In 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. 15 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 -1 . Moreover, the improvement is relevant to air temperature (from 2.74 to 2.67 K of RMSE) and precipitation (from 140 to 135 mm month -1 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.


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 spatio-temporal 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 25 atmospheric models. The surface layer (SL) occupies the lowest 10% of the ABL, 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 30 (e.g., Sellers et al., 1986Sellers et al., , 1996. The SL has two parts: inertial sublayer (ISL) and 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 35 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 40 increases, the regional and global models can be simulated with a finer spatial resolution and the grid size of the NWP model moving toward the gray zone of turbulence (the scales on the order of the energy-containing range of turbulence). Nevertheless, studies on the impact of a fine vertical resolution have not been relatively performed. In 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-45 3 times of the canopy height. Furthermore, model outputs are sensitive to the selection of the lowest model level height (Shin et al., 2011), 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 50 (e.g., Raupach, 1992;Physick and Garratt, 1995;Wenzel et al., 1997;Mölder et al., 1999;Harman andFinnigan, 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., 1986Sellers et al., , 1996. Despite the importance of the RSL, the Weather and Research Forecasting (WRF) model (Skamarock et 55 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 andFinnigan (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. 60 The parameterization of HFs has recently been incorporated in a one-dimensional (1D) 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 65 weather and climate simulations in terms of meteorological conditions near the Earth 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 is 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. 70

RSL theory of the HFs
The roughness sublayer parameterization by HFs is adopted herein along with an explanation of the core of the HF model, and the relevant details on this parameterization can be found in Harman andFinnigan (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 75 defines the vertical origin by the conventional zero-plane displacement height (d0). The same coordinate system is also applied herein. The vertical coordinates z and ̃ in this coordinate system are defined as the distance from d0 and from the terrain surface, respectively; therefore, their relation is =̃− & (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 logarithmic vertical profile: where k is von Kármán constant; C* is a C scale; C0 is C at z0; z0 is the roughness length; 4 is the integrated similarity function 80 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 (h) and the redefined zero-plane displacement height in the HFs ( 7 = ℎ − & ): the upper-canopy layer (z > dt), where the influence of additional mixing by the canopy exists, and the lower-canopy layer (z < dt), where the canopy is the direct source and sink for drag and heat. The vertical profile in the upper-canopy layer is described as follows: where 4 is the similarity function of C and E 4 is an RSL function of C. In the z→∞ limit, E 4 is equal to 1 and the wind speed converges to the MOST prediction. The last term in the right-hand side represents the additional mixing caused by the roughness element due to the coherent canopy turbulence, and can be replaced by E 4 , which is an integrated RSL function of C. The vertical profile from the HFs for the RSL deviates from that of MOST because of E 4 , thereby adjusting the logarithmic profile. The E 4 is introduced as follows: 90 (3) c1 and c2 are then determined from the continuity of E 4 at the canopy top and the RSL function, E 4 exponentially converges to zero above the RSL. In the lower canopy layer, C has the following exponential form: The RSL functions vary with atmospheric stability through β, where Lc is a canopy penetration depth defined as: where cd is a drag coefficient at the leaf scale and a is the leaf area density. The parameter dt and z0 also depend on the stability 95 because of their dependence on β:

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) The second step is to iteratively calculate the atmospheric stability (zr/L) as follows with an accuracy of 0.01: where u* n-1 is the friction velocity (u*) at previous time step. Equation (10) is different from Eq. (23) of Jiménez et al. (2012) by the RSL functions (i.e., E t and E Q ). After zr/L is determined, the third step is to iteratively update dt and β using Eqs. (5)  105 and (7) with an accuracy of 0.0001 because they are inter-correlated with each other. Subsequently, z0 is iteratively achieved with an accuracy of 0.0001 using Eq. (8) at the given zr/L, β, and dt. 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 for the upper-canopy layer. Equation (4) is used for the lower-canopy layer. Finally, * is given to and the aerodynamic conductance from z=0 to z = zr (ga) in the RSL are given to Computing time increased only 8 % with 61 vertical layers in the YSL scheme despite more iterations in the YSL scheme 110 compared to the revised MM5 SL scheme.

Numerical experimental design
This study evaluated the YSL scheme by making a 1D offline and real case simulations. The 1D offline simulations were done to test the YSL scheme performance without feedback with the atmosphere. The 1D offline simulation consisted of the YSL 115 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 120 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 (Figure 3). A 9 km horizontal resolution domain produced using the National Center for Environmental Prediction Final Analysis data (1° × 1°).

Observation data for the model validationevaluation
The model performance was examined against the surface wind speed, surface temperature and precipitation observed at 46 Automated Synoptic Observing System (ASOS) sites in Korea (Fig. 2). Quality control of the data includes gap detection,  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 the revised MM5 SL scheme, particularly in a smaller z/L (i.e., neutral and unstable conditions) and a larger Lc (i.e., small LAI and/or large h). The roughness length in a stable condition showed relatively smaller changes with z/L and Lc compared to those in the unstable condition. Our findings suggest that a small LAI in the winter season makes a larger Lc because of the smaller LAI, thereby leading to relatively larger 150 differences of z0 between the YSL scheme and the default WRF scheme. On the contrary, a similar value of z0 was observed in summer because of the larger LAI. Note that the revised MM5 SL scheme does not consider dt and β.
The RSL function, E 4 , was introduced to consider the additional mixing caused by the roughness element. Accordingly, E 4 should asymptotically converge to the MOST profile (i.e., E 4 ®1) as z increases with the continuous vertical profiles of the wind and the temperature. The YSL scheme reproduced these properties of E 4 and matched with the observed profiles inside 155 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 E 4 < 1 in the offRSL experiments. Furthermore, the YSL scheme produced wind and temperature within the canopy (i.e. ̃< & + & ), thereby providing additional useful information on the atmospheric dispersion inside the canopy.
The roughness length changes in the YSL scheme eventually produced changes in the surface energy balance with the 160 atmospheric stability (Fig. 7). In the 1D offline simulations based on the conditions in Table 1, the YSL scheme produced a larger z0 in the unstable and near-neutral conditions, but a smaller z0 in z/L > 3 compared to the offCTL. The aerodynamic conductance (ga) in the YSL scheme was larger in all the stability conditions even in the stable conditions in which the YSL provided a smaller z0 because the additional term in Eq. (13), ga2, dominated over the other effects in the ga calculation.
Accordingly, H and lE in the YSL scheme were larger than those in the revised MM5 SL scheme. Our finding implies stronger 165 fluxes from the YSL scheme when the gradient of quantity is the same. However, the impact of the increased ga was asymmetrical in H and lE depending on the soil moisture content. In this case simulation, an increase in lE 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 lE (Fig. S1).
A significant increase in lE was found along with a decrease in H in the strong unstable conditions (Fig. 7) because of the wet 170 soil moisture of 0.25 m 3 m -3 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 8 shows the real case simulation of the roughness length, 10 m wind speed (u10), and 2 m air temperature (T2). We discuss herein the real cases in the winter season because of stronger effect of the roughness sublayer. The results for the 175 summer season can be found in the Supplementary Materials. 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).

Real case simulations
Overall, the YSL scheme (rRSL experiment) produced 0.2-2.0 m larger z0 than the default values in the rCTL experiment over the tall canopies, where Lc was large. In contrast, the YSL scheme produced a similar or even slightly smaller z0 over the short canopies compared to the rCTL experiment. Importantly, the changes of z0 made direct impacts on the momentum fluxes 180 and thus surface wind speed (Fig. 8b). The typical u10 in the rCTL was larger than approximately 3 m s −1 , and a much stronger wind (> 6 m s −1 ) 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, u10 in the rRSL was reduced by approximately 30%; however, the region of the increased wind speed provided 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 m s −1 to 1.0 m s −1 and from 3.1 m s −1 and 1.8 m s −1 . 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 190 (Table 2, Fig. 10).
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 ga, particularly in the forest canopies and mountain regions (Fig. 11a).
This larger ga in the YSL scheme led to the increases of the latent heat fluxes by approximately 20 W m −2 , with an eventual reduction of the soil water content (Fig. 12a). The sensible heat fluxes in the rCTL experiments were generally approximately 195 80 W m −2 , except over the snow-covered region where H was approximately 40 W m −2 . 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 + lE) increased in the YSL scheme, and a larger lE in the rRSL led to a cooler temperature than that in the rCTL experiment (Fig. 8c).
During the winter simulation period, precipitation was observed over an extensive area in the domain, and snow was 200 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 lE, 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 205 default WRF simulation because our implemented RSL parameterization started to converge to the default WRF in a smaller Lc (i.e., larger LAI and/or smaller h) and strong synoptic influences by the summer heavy rainy period (Table S1, Figs. S2-S6).

Summary and conclusion remark
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 210 Earth's surface. MOST can be only applicable in the inertial layer and turbulence deviates from MOST in the roughness sublayer. Importantly, 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 215 simulations. For these purposes, we designed a series of offline simulations with an idealized boundary condition and a real case simulation to evaluate the performance of the YSL scheme against the observation data.
The 1D offline simulation revealed that the YSL scheme successfully reproduced the features reported in various canopies.
The RSL function, E 4 , asymptotically increased to 1, and the vertical gradients of the wind speed and the temperature decreased 220 in the RSL as z increased, thereby converging to the MOST prediction. Notably, unlike the typical assignment of the roughness parameters (i.e., z0, dt, and β) as a constant, the roughness parameters are functions of the atmospheric stability (z/L) and Lc.
The roughness parameters had a maximum in the weakly unstable condition and in larger Lc (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. 225 Meanwhile, the real case simulation showed that the RSL-incorporated WRF produced a larger z0 than the default WRF.
This increase in z0 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 z0 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 z0 also made increases in the available energy. This increased 230 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 235 time (e.g., Hu et al., 2010Hu et al., , 2013Shimada et al., 2012;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.
Code and data availability. Effective heat transfer coefficient for nonturbulent processes (Carlson and Boland, 1978;Jiménez et al., 2012) Variable at , such as and

Scale of
Conventionally defined zero-plane displacement height Redefined zero-plane displacement height in Harman and Finnigan (2007) Parameter related the depth scale of the scalar profile (º = L (¯1 + 4 4 ³ − 1))