Evaluating the performance of two surface layer schemes for the momentum and heat exchange processes during severe haze pollution in Jing-Jin-Ji in eastern China

The turbulent flux parameterization schemes in the surface layer are crucial for air pollution modeling. There have been some deficiencies in the prediction of air pollutants by atmosphere chemical models, which is closely related to the uncertainties of the momentum and sensible heat fluxes calculated in the surface layer. The differences between two surface layer schemes (Li and MM5 schemes) were discussed, and the performances of two schemes were mainly evaluated based on the observed momentum and sensible heat fluxes during a heavy haze episode in Jing-Jin-Ji in eastern China. The results showed that the aerodynamic roughness length z_0m and the thermal roughness length z_0h played major roles in the flux calculation. Compared with the Li scheme, ignoring the difference between z_0m and z_0h in the MM5 scheme induced a great error in the calculation of the sensible heat flux (e.g., the error was 54 % at Gucheng station). Besides the roughness length, the algorithm for the surface turbulent flux as well as the roughness sublayer also resulted in certain errors in the MM5 scheme. In addition, magnitudes of z_0m and z_0h have significant influence on the two schemes. The large z_0m and z_0m∕z_0h in megacities with a rough surface (e.g., Beijing) resulted in much larger differences of momentum and sensible heat fluxes between Li and MM5, compared with the small z_0m and z_0m∕z_0h in suburban areas with a smooth surface (e.g., Gucheng). The Li scheme could better characterize the evolution of atmospheric stratification than the MM5 scheme in general, especially for the transition stage from unstable to stable atmospheric stratification, corresponding to the PM_2.5 accumulation. The biases of momentum and sensible heat fluxes from Li were lower, about 38 % and 43 %, respectively, than those from MM5 during this stage. This study indicates the superiority of the Li scheme in describing regional atmospheric stratification and an improved possibility of severe haze prediction in Jing-Jin-Ji in eastern China by coupling it into atmosphere chemical models.

1 Introduction

Adequate air quality modeling relies on the accurate simulation of meteorological conditions, especially in the planetary boundary layer (PBL) (Hu et al., 2010; Cheng et al., 2012; Xie et al., 2012). The PBL is tightly coupled with the earth's surface through turbulent exchange processes. As the bottom layer of the PBL, the surface layer (SL) reflects the surface state by calculating momentum, heat, water vapor, and other fluxes, and influences the atmospheric structure through a turbulent transport process. Many studies have illustrated the important roles of meteorological factors in the SL during air pollution formation. It has been demonstrated that weak wind speed, high relative humidity (RH), and strong temperature inversion are favorable for the concentration of haze (Zhang et al., 2014; Yang et al., 2015; Liu et al., 2017; Zhong et al., 2018). Strong stable stratification and weak turbulence are mainly responsible for many haze events. The relationship between the flux and the atmospheric profile in the atmospheric surface layer is a critical factor for air pollution diffusion, especially under stable stratification conditions (Li et al., 2017). However, there are still some uncertainties in the study of the stable boundary layer due to the poor description of surface turbulent motion. The simulation study on a severe haze in eastern China by the Weather Research and Forecasting with Chemistry (WRF-Chem) model concluded that current PBL schemes have a weak ability to distinguish between haze days under stable conditions and clean days under unstable conditions (T. Li et al., 2016). Another study (Vautard et al., 2012) of mesoscale meteorological models also pointed out there was a systematic overestimation of near-surface wind speed in the stable boundary layer which should contribute to the underestimation of surface concentrations of primary pollution. In addition, atmospheric conditions in both the PBL and upper layers are highly dependent on turbulent fluxes which are computed in the SL (Ban et al., 2010). Flux parameterization in the SL plays an important role in studies of the hydrological cycle and weather prediction (Yang et al., 2001; Li, 2014). An adequate SL scheme is crucial in providing an accurate atmospheric evolution by numerical models (Jiménez et al., 2012) and hence it may introduce significant impacts on air pollution simulation.

The bulk aerodynamic formulation based on Monin–Obukhov similarity theory (hereinafter MOST; Monin and Obukhov, 1954) is usually employed to calculate surface fluxes in numerical models. Turbulent fluxes are parameterized by wind, temperature, humidity in the lowest layer in the model, and temperature and humidity at the surface. Many international scholars verified MOST using field experiments and then proposed the universal functions, the most commonly used of which is the Businger–Dyer (BD) equation (Businger, 1966; Dyer, 1967). With the development of observation technology, the coefficients in the BD equation have been further modified (Paulson, 1970; Webb, 1970; Businger et al., 1971; Dyer, 1974; Högström, 1996). In addition to the BD equation, some other schemes have been put forward and they performed better, especially for strongly stable stratification (Holtslag and De Bruin, 1988; Beljaars and Holtslag, 1991; Cheng and Brutsaert, 2005). The schemes can be divided into two types according to the computing characteristics. One type is called an iterative algorithm (Paulson, 1970; Businger et al., 1971; Dyer, 1974; Högström, 1996; Beljaars and Holtslag, 1991), and it keeps MOST completely with less approximation so that the results can be more precise. However, many more steps are necessary for it to converge, and hence the CPU time is longer, which reduces the computational efficiency of modeling (Louis, 1979; Li et al., 2014). The other one is called a non-iterative algorithm (Louis et al., 1982; Launiainen, 1995; Wang et al., 2002; Wouters et al., 2012). There is no requirement for loop iteration in the calculation due to the approximate treatment. This algorithm is much simpler and less CPU-time-consuming, but the results are based on the loss of the calculation accuracy.

A new non-iterative scheme proposed by Li et al. (2014, 2015; Li hereinafter) speeds up effectively with a higher accuracy compared with classic iterative computation. It is remarkable that this new scheme has just been theoretically evaluated and it has never been applied in any models. Haze pollution has occurred frequently in recent years in eastern China. The concentration of PM_2.5 may reach up to 1000 µg m⁻³ in the Beijing–Tianjin–Hebei (Jing-Jin-Ji) region in winter (Wang et al., 2014), while it is generally underestimated by current air quality models (Zhang et al., 2015; T. Li et al., 2016; Liu et al., 2017). The Li and another classic SL scheme (Zhang and Anthes, 1982; MM5 hereinafter) were compared in detail in this study. The observed momentum and sensible heat flux data covering one complete haze process at Gucheng station were used to evaluate the two schemes, focusing on the transition stage from unstable to stable atmospheric stratification, corresponding to the PM_2.5 accumulation. The evaluation is in view of both local and regional scales. This study may provide the prerequisite for coupling the Li scheme into atmosphere chemical models in the future.

2 Theory

The definitions of momentum and sensible heat fluxes as well as the detailed algorithms of the Li and MM5 schemes are introduced in this section.

2.1 Introduction of momentum and sensible heat fluxes

The turbulent fluxes from the ground surface are defined as follows:

$$\begin{array}{}\text{(1a)}& {\displaystyle }& {\displaystyle }\mathit{\tau }=\mathit{\rho }{u}_{\ast }^{\mathrm{2}},\text{(1b)}& {\displaystyle }& {\displaystyle }H=-\mathit{\rho }{c}_p{u}_{\ast }{\mathit{\theta }}_{\ast },\end{array}$$

where τ is the momentum flux, H is the sensible heat flux, ρ is the air density, c_p is the specific heat capacity at constant pressure; u_∗ and θ_∗ are the friction velocity and the temperature scale, respectively, and they represent the intensity of the vertical turbulent flux transport and are approximately independent of height in the SL.

Both the Li and MM5 schemes are based on bulk flux parameterization. As an important dimensionless parameter related to the stability, the bulk Richardson number Ri_B is defined as

$$\begin{array}{}\text{(2)}& {\mathit{Ri}}_{\mathrm{B}}={\displaystyle \frac{gz\left(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{g}}\right)}{\mathit{\theta }{u}^{\mathrm{2}}}},\end{array}$$

where g is the acceleration of gravity, z is the reference height which is the lowest level in models, θ is the mean potential temperature at height z, θ_g is the surface radiometric potential temperature, and u is the mean wind speed at height z. Thus, Ri_B can be computed through meteorological variables from at least two levels.

2.2 The Li scheme

This new scheme employs a non-iterative algorithm to compute the surface fluxes. Its basic idea is to parameterize the stability parameter ζ directly with Ri_B and roughness lengths (z_0m and z_0h). Specifically, bulk transfer coefficients of the momentum and sensible heat fluxes (C_M and C_H) are expressed as

$$\begin{array}{}\text{(3a)}& {\displaystyle }{C}_{\mathrm{M}}& {\displaystyle }={\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{u^{\mathrm{2}}}}={\displaystyle \frac{\mathit{\tau }}{\mathit{\rho }{u}^{\mathrm{2}}}},\text{(3b)}& {\displaystyle }{C}_{\mathrm{H}}& {\displaystyle }={\displaystyle \frac{u_{\ast }{\mathit{\theta }}_{\ast }}{u\left(\mathit{\theta }-{\mathit{\theta }}_g\right)}}=-{\displaystyle \frac{H}{\mathit{\rho }{c}_{p}u\left(\mathit{\theta }-{\mathit{\theta }}_g\right)}}.\end{array}$$

Based on MOST and considering the roughness sublayer (RSL) effect at the same time, the relationships between the bulk transfer coefficients and the profile functions corresponding to wind and potential temperature are usually expressed as

$$\begin{array}{}\text{(4a)}& {\displaystyle }{C}_{\mathrm{M}}& {\displaystyle }={\displaystyle \frac{k^{\mathrm{2}}}{{\left[\ln \frac{z}{z_{\mathrm{0}m}}-{\mathit{\psi }}_{\mathrm{M}}\left(\frac{z}{L}\right)+{\mathit{\psi }}_{\mathrm{M}}\left(\frac{z_{\mathrm{0}m}}{L}\right)+{\mathit{\psi }}_{\mathrm{M}}^{\ast }\left(\frac{z}{L},\frac{z}{z_{\ast }}\right)\right]}^{\mathrm{2}}}},\text{(4b)}& {\displaystyle }{C}_{\mathrm{H}}& {\displaystyle }={\displaystyle \frac{k^{\mathrm{2}}}{\begin{array}{c}R\left[\ln \frac{z}{z_{\mathrm{0}m}}-{\mathit{\psi }}_{\mathrm{M}}\left(\frac{z}{L}\right)+{\mathit{\psi }}_{\mathrm{M}}\left(\frac{z_{\mathrm{0}m}}{L}\right)+{\mathit{\psi }}_{\mathrm{M}}^{\ast }\left(\frac{z}{L},\frac{z}{z_{\ast }}\right)\right]\\ \left[\ln \frac{z}{z_{\mathrm{0}h}}-{\mathit{\psi }}_{\mathrm{H}}\left(\frac{z}{L}\right)+{\mathit{\psi }}_{\mathrm{H}}\left(\frac{z_{\mathrm{0}h}}{L}\right)+{\mathit{\psi }}_{\mathrm{H}}^{\ast }\left(\frac{z}{L},\frac{z}{z_{\ast }}\right)\right]\end{array}}},\end{array}$$

where k is the von Kármán constant, which is 0.4 in both schemes; R is the Prandtl number, which is 1.0 in the two schemes; z_0m and z_0h are the aerodynamic roughness length and the thermal roughness length, respectively; ψ_M and ψ_H are the integrated stability functions for momentum and sensible heat, respectively, which are also called universal functions. L is the Obukhov length ($\mathit{\zeta }=\frac{z}{L})$, ${\mathit{\psi }}_{\mathrm{M}}^{\ast }$ and ${\mathit{\psi }}_{\mathrm{H}}^{\ast }$ are the correction functions accounting for the RSL effect, and z_∗ is the RSL height. It is clear to see that the calculation of the momentum and sensible heat fluxes requires C_M and C_H (or u_∗ and θ_∗), and there are three key points to consider in obtaining them:

1. z_0m and z_0h are two key parameters in the bulk transfer equations. Their definitions and influences will be discussed in Sect. 4.1. Note that both z_0m and z_0h are taken into account by the Li scheme. In other words, the Li scheme distinguishes the two principal surface parameters effectively as they generate from different mechanisms.

2. The determination of ζ is the most crucial problem in the Li scheme. In fact, this new scheme consists of two parts. The first part is proposed for atmospheric stable stratification conditions (Li et al., 2014), and the second part then extends the scheme to unstable conditions (Li et al., 2015). For stable conditions, the calculation procedure for a given group of Ri_B, z_0m, and z_0h is as follows: (1) find the region according to z_0m and z_0h, (2) find the section according to the region and Ri_B with Eq. (5) and given coefficients, and (3) calculate ζ using Eq. (6) and the given coefficients.

   $$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{Ri}}_{\mathrm{Bcp}}=\sum C_{mn}{\left(\log L_{\mathrm{0}M}\right)}^m{\left(L_{\mathrm{0}H}-L_{\mathrm{0}M}\right)}^n,\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }={\mathit{Ri}}_{\mathrm{B}}\sum C_{ijk}{\mathit{Ri}}_{\mathrm{B}}^i{L}_{\mathrm{0}M}^j{\left(L_{\mathrm{0}H}-L_{\mathrm{0}M}\right)}^k,\end{array}$$

   where C_mn and C_ijk are the coefficients listed in tables in Li et al. (2014). $L_{\mathrm{0}M}=\ln \frac{z}{z_{\mathrm{0}m}}$, $L_{\mathrm{0}H}=\ln \frac{z}{z_{\mathrm{0}h}}$. m, n=0, 1, 2, and $m+n\le \mathrm{3}$; i, j, k=0, 1, 2, 3, and $i+j+k\le \mathrm{4}$. Similarly, for unstable conditions, eight regions are divided according to the method from Li et al. (2015). For each of the regions, ζ is carried out as follows:

   $$\begin{array}{}\text{(7)}& \mathit{\zeta }={\mathit{Ri}}_{\mathrm{B}}{\displaystyle \frac{L_{\mathrm{0}M}^{\mathrm{2}}}{L_{\mathrm{0}H}}}\sum C_{ijk}{\left({\displaystyle \frac{-{\mathit{Ri}}_{\mathrm{B}}}{\mathrm{1}-{\mathit{Ri}}_{\mathrm{B}}}}\right)}^i{L}_{\mathrm{0}M}^{-j}{\mathrm{L}}_{\mathrm{0}H}^{-k},\end{array}$$

   where C_ijk is listed in Y. Li et al. (2016), and i=0, 1; j, k=0, 1, 2, 3; $i+j+k\le \mathrm{4}$.

3. The universal function is also a key factor in flux calculation. The form of universal function here is adopted from Cheng and Brutsaert (2005) under stable conditions (Eq. 8a, b) and it is adopted from Paulson (1970) under unstable conditions (Eq. 9a, b):

   $$\begin{array}{}\text{(8a)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{M}}\left(\mathit{\zeta }\right)=-a\ln \left[\mathit{\zeta }+{\left(\mathrm{1}+{\mathit{\zeta }}^b\right)}^{\frac{\mathrm{1}}{b}}\right],\mathit{\zeta }>\mathrm{0}\text{ (stable)},\text{(8b)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{H}}\left(\mathit{\zeta }\right)=-c\ln \left[\mathit{\zeta }+{\left(\mathrm{1}+{\mathit{\zeta }}^d\right)}^{\frac{\mathrm{1}}{d}}\right],\mathit{\zeta }>\mathrm{0}\text{ (stable)},\end{array}$$
   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{M}}\left(\mathit{\zeta }\right)=\mathrm{2}\ln {\displaystyle \frac{\mathrm{1}+x}{\mathrm{2}}}+\ln {\displaystyle \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}}-\mathrm{2}\arctan \left(x\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\\ \text{(9a)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }<\mathrm{0}\text{ (unstable)},\text{(9b)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{H}}\left(\mathit{\zeta }\right)=\mathrm{2}\ln {\displaystyle \frac{\mathrm{1}+y}{\mathrm{2}}},\mathit{\zeta }<\mathrm{0}\text{ (unstable)},\end{array}$$

   where a=6.1, b=2.5, c=5.3, d=1.1, $x={\left(\mathrm{1}-\mathrm{16}\mathit{\zeta }\right)}^{\mathrm{1}/\mathrm{4}}$, $y={\left(\mathrm{1}-\mathrm{16}\mathit{\zeta }\right)}^{\mathrm{1}/\mathrm{2}}$.

In addition, the RSL effect is taken into account in the Li scheme. The definition and influence of the RSL will also be discussed in Sect. 4.1. De Ridder (2010) proposed the expression of ${\mathit{\psi }}_{\mathrm{M}}^{\ast }$ and ${\mathit{\psi }}_{\mathrm{H}}^{\ast }$:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{M}}^{\ast }\left(\mathit{\zeta },{\displaystyle \frac{z}{z_{\ast }}}\right)={\mathit{\varphi }}_{\mathrm{M}}\left[\left(\mathrm{1}+{\displaystyle \frac{\mathit{\upsilon }}{{\mathit{\mu }}_{\mathrm{M}}z/z_{\ast }}}\right)\mathit{\zeta }\right]\\ \text{(10a)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{{\mathit{\mu }}_{\mathrm{M}}z/z_{\ast }}}\right)e^{-{\mathit{\mu }}_{\mathrm{M}}z/z_{\ast }},{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{H}}^{\ast }\left(\mathit{\zeta },{\displaystyle \frac{z}{z_{\ast }}}\right)={\mathit{\varphi }}_{\mathrm{H}}\left[\left(\mathrm{1}+{\displaystyle \frac{\mathit{\upsilon }}{{\mathit{\mu }}_{\mathrm{H}}z/z_{\ast }}}\right)\mathit{\zeta }\right]\\ \text{(10b)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{{\mathit{\mu }}_{\mathrm{H}}z/z_{\ast }}}\right)e^{-{\mathit{\mu }}_{\mathrm{H}}z/z_{\ast }},\end{array}$$

where υ=0.5, μ_M=2.59, μ_H=0.95, $z_{\ast }=\mathrm{16.7}{z}_{\mathrm{0}m}$, λ=1.5; ϕ_M and ϕ_H are universal functions before integration. Here, set ${\mathit{\chi }}_{\mathrm{M}}=\mathrm{1}+\frac{\mathit{\upsilon }}{{\mathit{\mu }}_{\mathrm{M}}z/z_{\ast }}$, ${\mathit{\chi }}_{\mathrm{H}}=\mathrm{1}+\frac{\mathit{\upsilon }}{{\mathit{\mu }}_{\mathrm{H}}z/z_{\ast }}$:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{M}}\left({\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)=\mathrm{1}+a{\displaystyle \frac{{\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }+{\left({\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)}^b{\left[\mathrm{1}+{\left({\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)}^b\right]}^{\frac{\mathrm{1}-b}{b}}}{{\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }+{\left[\mathrm{1}+{\left({\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)}^b\right]}^{\frac{\mathrm{1}}{b}}}},\\ \text{(11a)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }>\mathrm{0}\text{ (stable)},{\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{H}}\left({\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)=\mathrm{1}+c{\displaystyle \frac{{\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }+{\left({\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)}^d{\left[\mathrm{1}+{\left({\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)}^d\right]}^{\frac{\mathrm{1}-d}{d}}}{{\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }+{\left[\mathrm{1}+{\left({\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)}^d\right]}^{\frac{\mathrm{1}}{d}}}},\\ \text{(11b)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }>\mathrm{0}\text{ (stable)},\end{array}$$

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{M}}\left({\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)={\left(\mathrm{1}-\mathrm{16}{\mathit{\chi }}_{\mathrm{M}}\mathit{\zeta }\right)}^{-\mathrm{1}/\mathrm{4}},\mathit{\zeta }<\mathrm{0}\text{ (unstable)},\text{(12b)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{H}}\left({\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)={\left(\mathrm{1}-\mathrm{16}{\mathit{\chi }}_{\mathrm{H}}\mathit{\zeta }\right)}^{-\mathrm{1}/\mathrm{2}},\mathit{\zeta }<\mathrm{0}\text{ (unstable)}.\end{array}$$

2.3 The MM5 scheme

The MM5 scheme is a classic one which is widely applied in modeling investigations (Hu et al., 2010; Wang et al., 2015a, b; Tymvios et al., 2017). This scheme does not distinguish z_0h from z_0m; thus the roughness length here is expressed as z₀. For unstable conditions, the function forms are given by Eq. (16a, b) following Paulson (1970), and for stable conditions, the atmospheric stratification conditions are subdivided into three cases according to Zhang and Anthes (1982) and the function forms are given by Eqs. (13), (14), and (15).

1. Strongly stable condition (Ri_B≥0.2):

   $$\begin{array}{}\text{(13)}& {\mathit{\psi }}_{\mathrm{M}}={\mathit{\psi }}_{\mathrm{H}}=-\mathrm{10}\ln {\displaystyle \frac{z}{z_{\mathrm{0}}}}.\end{array}$$

2. Weakly stable condition ($\mathrm{0}<{\mathit{Ri}}_{\mathrm{B}}<\mathrm{0.2}$):

   $$\begin{array}{}\text{(14)}& {\mathit{\psi }}_{\mathrm{M}}={\mathit{\psi }}_{\mathrm{H}}=-\mathrm{5}\left({\displaystyle \frac{{\mathit{Ri}}_{\mathrm{B}}}{\mathrm{1.1}-\mathrm{5}{\mathit{Ri}}_{\mathrm{B}}}}\right)\ln {\displaystyle \frac{z}{z_{\mathrm{0}}}}.\end{array}$$

3. Neutral condition (Ri_B=0):

   $$\begin{array}{}\text{(15)}& {\mathit{\psi }}_{\mathrm{M}}={\mathit{\psi }}_{\mathrm{H}}=\mathrm{0}.\end{array}$$

4. Unstable condition (Ri_B<0):

   $$\begin{array}{}\text{(16a)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{M}}=\mathrm{2}\ln {\displaystyle \frac{\mathrm{1}+x}{\mathrm{2}}}+\ln {\displaystyle \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}}-\mathrm{2}\arctan \left(x\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\text{(16b)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{H}}=\mathrm{2}\ln {\displaystyle \frac{\mathrm{1}+y}{\mathrm{2}}},\end{array}$$

   where $x={\left(\mathrm{1}-\mathrm{16}\mathit{\zeta }\right)}^{\mathrm{1}/\mathrm{4}}$, $y={\left(\mathrm{1}-\mathrm{16}\mathit{\zeta }\right)}^{\mathrm{1}/\mathrm{2}}$.

This scheme calculates turbulent fluxes of the momentum and sensible heat with u_∗ and θ_∗. In order to avoid the huge difference between the two computations, u_∗ is arithmetically averaged with its previous value by Eq. (17), and a lower limit of $u_{\ast }=\mathrm{0.1}$ m s⁻¹ is imposed to prevent the heat flux from being zero under very stable conditions. According to the profile functions of wind and temperature near the ground, θ_∗ is then deduced by Eq. (18).

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }{u}_{\ast }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(u_{\ast }+{\displaystyle \frac{ku}{\ln \frac{z}{z_{\mathrm{0}m}}-{\mathit{\psi }}_{\mathrm{M}}}}\right),\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\ast }{\displaystyle \frac{k\left(\mathit{\theta }-{\mathit{\theta }}_g\right)}{R[\ln \frac{z}{z_{\mathrm{0}h}}-{\mathit{\psi }}_{\mathrm{H}}]}}.\end{array}$$

The calculation procedure of the Li scheme is as follows: (1) determine Ri_B, z_0m, and z_0h according to the observation data; (2) calculate ζ with Ri_B, z_0m, and z_0h; and (3) calculate the momentum and sensible heat fluxes under different conditions. The MM5 scheme is summarized as follows: (1) determine the universal functions according to the values of Ri_B and z₀, (2) calculate the u_∗ and θ_∗ with the meteorological variables and flux data, and (3) derive the turbulent fluxes. Compared with other non-iterative schemes including MM5, the Li scheme can be applied to the full range of roughness status $\mathrm{10}\le \frac{z}{z_{\mathrm{0}m}}\le {\mathrm{10}}^{\mathrm{5}}$ and $-\mathrm{0.5}\le \ln \frac{z_{\mathrm{0}m}}{z_{\mathrm{0}h}}\le \mathrm{30}$ under whole conditions $-\mathrm{5}\le {\mathit{Ri}}_{\mathrm{B}}\le \mathrm{2.5}$. In addition, there are three obvious differences between the Li and MM5 schemes: (1) Li distinguishes z_0h from z_0m but MM5 does not, (2) the two schemes apply different universal functions under stable conditions, and (3) Li considers the RSL effect while MM5 ignores it.

Figure 1 Location (a) and geographical environment (b) at Gucheng station. The map is from Bing Maps.

3 Observational data and methods

The observational fluxes used in this study were measured at Gucheng station from 1 December 2016 to 9 January 2017. Gucheng station (115.40^∘ E, 39.08^∘ N) is located in Gucheng County, Baoding, Hebei Province, and it is about 110 km southwest of Beijing (Fig. 1a). This station is on a farmland site where rice is grown in summer and wheat in winter. The surroundings are mainly farmland and scattered villages (Fig. 1b). At Gucheng station, the momentum and sensible heat fluxes near the surface were measured by the eddy correlation flux measurement system. The system is mainly composed of a sonic anemometer (CSAT3) and a gas analyzer (LI-7500). They are set up at a height of 4 m above the ground surface. The measured fluxes are used to evaluate the two schemes as well as estimate the roughness lengths. The measured meteorological variables including wind speed and direction, temperature, humidity, pressure, and radiation are utilized to calculate the momentum and sensible heat fluxes both in the Li and MM5 schemes. Note the observed meteorological data were from Gucheng station and national basic automatic weather stations in Jing-Jin-Ji in eastern China, respectively. Hourly surface PM_2.5 mass concentration in Baoding and Beijing from the China National Environmental Monitoring Centre (http://www.cnemc.cn/, last access: 14 October 2018) was also used in this paper.

Figure 2 Wind rose map at Gucheng station from 1 December 2016 to 9 January 2017.

3.1 Data processing

To obtain accurate flux data, quality control has been performed for the observational data, including (1) the elimination of the outliers and the data on rainy days, (2) double rotation and WPL correction (Webb et al., 1980), and (3) omission of the dataset when the wind speed is less than 0.5 m s⁻¹. In addition, the wind field, especially the wind direction, has a great impact on the value of z_0m, so it is necessary to understand the situation at Gucheng station. Figure 2 shows the distribution frequency of wind speed and wind direction at Gucheng during the observation (1 December 2016–9 January 2017). The wind speed is stable during this period and the maximum is no more than 5 m s⁻¹, and most of them are about 1–2 m s⁻¹. The wind direction is relatively uniform except for the southeast wind (135^∘).

Figure 3 The surface emissivity ε_s dependence of RMSE between observed near-neutral heat fluxes and parameterized heat fluxes (red for Li and blue for MM5) at Gucheng station.

3.2 Determination of surface skin temperature

The surface skin temperature at Gucheng station is calculated from the radiation data by the following formula:

$$\begin{array}{}\text{(19)}& R_{\mathrm{lw}}^{\uparrow }=\left(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{s}}\right)R_{\mathrm{lw}}^{\downarrow }+{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_g^{\mathrm{4}},\end{array}$$

where $R_{\mathrm{lw}}^{\uparrow }$ and $R_{\mathrm{lw}}^{\downarrow }$ are the surface upward longwave radiation and longwave radiation incident on the surface, respectively; σ is the Stephen Boltzmann constant, $\mathit{\sigma }=\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$. T_g is the surface skin temperature and ε_s is the surface emissivity, which is the prerequisite of T_g calculation. Many research studies estimated the value of ε_s and found it is always 0.9–1 (Stewart et al., 1994; Verhoef et al., 1997). According to the semi-empirical method in Yang et al. (2008), ε_s is estimated when the root mean square error (RMSE) is minimal. In this paper, the Li and MM5 schemes were used to estimate the ε_s value (as shown in Fig. 3). It is clear that the ε_s value corresponding to the minimum RMSE is not very sensitive to the choice of the two schemes. When ε_s is 1, the RMSE has a minimum value. Thus, this experiment takes 1 as the optimal value of ε_s.

3.3 Determination of roughness length z_0m (z_0h)

Using the observed momentum and sensible heat fluxes and the meteorological variables including wind speed, temperature, humidity, and pressure after quality control at Gucheng station, z_0m and z_0h were derived from Eq. (20a, b) following Yang et al. (2003) and Sicart et al. (2014).

$$\begin{array}{}\text{(20a)}& {\displaystyle }& {\displaystyle \frac{u_{\ast }}{u}}={\displaystyle \frac{k}{\ln \frac{z}{z_{\mathrm{0}m}}-{\mathit{\psi }}_{\mathrm{M}}}},\text{(20b)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\theta }}_{\ast }}{\left(\mathit{\theta }-{\mathit{\theta }}_g\right)}}={\displaystyle \frac{k}{R[\ln \frac{z}{z_{\mathrm{0}h}}-{\mathit{\psi }}_{\mathrm{H}}]}}.\end{array}$$

During the observation period, the crops stopped growing and the height did not exceed 0.1 m, so the zero-plane displacement height was ignored and the reference height z was taken as 4 m. The observation time was too short (about 1 month) to consider the effect of seasonal variations on the roughness length. Thus, z_0m and z_0h were assumed to be two fixed values. Based on the variables and formulae mentioned above, the two roughness lengths at Gucheng are derived: z_0m=0.0419 m and z_0h=0.0042 m.

4 Results and discussion

The definitions and influences of the RSL on the calculation of the turbulent flux are discussed in detail in this section. The Li and MM5 schemes are tested offline and evaluated during the haze pollution from 13 to 23 December 2016.

Table 1 Typical values of z_0m corresponding to various land-cover types.

4.1 The influences of the RSL and roughness length on the calculation of the turbulent flux

The RSL is usually defined as the region where the flow is influenced by the individual roughness elements as reflected by the spatial inhomogeneity of the mean flow (Florens et al., 2013). In the RSL, turbulence is strongly affected by individual roughness elements, and standard MOST is no longer valid (Simpson et al., 1998). Therefore, it is necessary to consider the RSL effect in the calculation of the turbulent flux, especially for rough terrain such as forest or large cities. z_0m is defined as the height at which the extrapolated wind speed following the similarity theory vanishes. It is mainly determined by land-cover type and canopy height after excluding large obstructions. In models, z_0m is always based on the look-up table which is related to land-cover types. In this study, z_0m is simply classified based on the research of Stull (1988) and listed in Table 1. It can be seen in Table 1 that the rougher underlying surface corresponds to the larger value of z_0m. z_0h is the height at which the extrapolated air temperature is identical to the surface skin temperature. Some early researchers assumed that z_0m was equal to z_0h (Louis, 1979; Louis et al., 1982). However, the assumption is not applicable in reality because z_0m and z_0h have different physical meanings. Different treatments of z_0m and z_0h may introduce considerable changes in the surface flux calculation (Launiainen, 1995; Kot and Song, 1998; Anurose and Subrahamanyam, 2013). Many studies removed the assumption that z_0m was equal to z_0h and made the schemes more applicable in the situation that z_0m was not equal to z_0h or the ratio of z_0m to z_0h was much larger (Wouters et al., 2012; Li et al., 2014, 2015). Some field experiments even indicated the ratio z_0m∕z_0h has a diurnal variation (Sun, 1999; Yang, 2003, 2008). In this study, we make the common assumption that the ratio z_0m∕z_0h is a constant.

Figure 4 The relationships between C_M(C_H) and Ri_B under different z_0m values and treatments of the RSL. Solid lines: considering the RSL effect; dotted lines: without the RSL effect.

Considering the lowest level in mesoscale models is usually about 10 m, z=10 m is set as the reference height in this study. The range of Ri_B is set according to Louis et al. (1982) in the following discussion. Firstly, the study discusses the effects of different land-cover types (different z_0m values) and the RSL on flux calculation. Set z_0m=z_0h, corresponding to four cases: z_0m=1, 0.5, 0.05, and 0.001 m. These cases correspond to large cities, forests, agricultural fields, and wide water surface, respectively. Figure 4 shows the relationship between C_M(C_H) and Ri_B under different z_0m values and treatments of the RSL. It can be seen that both the RSL and z_0m have impacts on C_M and C_H. Ignoring the RSL effect can result in larger C_M and C_H compared with the results of the original scheme considering the RSL effect. The difference induced by the RSL effect is only evident under the rough surface. For example, the difference under z_0m=1 is obviously greater than other z_0m settings, and when z_0m is reduced to 0.05 or less, the RSL has little effect. Furthermore, the RSL contributes more to sensible heat transfer than to momentum transfer under the same setting of z_0m. The effects of different land-cover types on C_M and C_H are much more significant compared with the RSL. The rougher surface (corresponding to the larger z_0m value) brings the larger C_M (C_H) under the same stability. In addition, there is a corresponding relationship between C_M(C_H) and stability. The value of C_M(C_H) drops with stability. Once Ri_B exceeds the critical value (generally 0.2–0.25), the transfer coefficients decline sharply but are still above 0.

Figure 5 The relationships between C_M(C_H) and Ri_B under different ratios of z_0m to z_0h and treatments of the RSL. Solid lines: considering the RSL effect; dotted lines: without the RSL effect.

Secondly, the effects of difference between z_0m and z_0h as well as the RSL on flux calculation are discussed. The relationship between z_0m and z_0h can be expressed as $k{B}^{-\mathrm{1}}=\ln \frac{z_{\mathrm{0}m}}{z_{\mathrm{0}h}}$. Over the sea, z_0m is comparable to z_0h; over the uniform vegetation surface (e.g., grassland, farmland, woodland), kB⁻¹ is about 2 ($z_{\mathrm{0}m}/z_{\mathrm{0}h}\approx \mathrm{10}$) (Garratt and Hicks, 1973; Garratt, 1978; Garratt and Francey, 1978), which coincides with our results in Gucheng (z_0m=0.0419 m, z_0h=0.0042 m); over the surface with bluff roughness elements, the kB⁻¹ value may be very large. For example, in some large cities, kB⁻¹ is even up to 30 ($z_{\mathrm{0}m}/z_{\mathrm{0}h}\approx {\mathrm{10}}^{\mathrm{13}}$) (Sugawara and Narita, 2009). Therefore, the ratio z_0m∕z_0h varies over a wide range. Figure 5 shows the relationship between C_M(C_H) and Ri_B under different treatments of z_0m∕z_0h. Set z_0m=1 as a large city case, z_0h=1, 0.01, 10⁻⁴, and 10⁻⁶ m, and the large differences derived from the different ratios are displayed in Fig. 5. The differences induced by the RSL effect are more obvious than those in Fig. 4. The different treatments of ratio z_0m∕z_0h have great impacts on turbulent flux transfer, particularly for sensible heat transfer. It seems evident that when z_0h is not equal to z_0m ($z_{\mathrm{0}m}/z_{\mathrm{0}h}=\mathrm{100}$–10⁶), the calculated C_H is much smaller compared to the treatment when z_0h is equal to z_0m ($z_{\mathrm{0}m}/z_{\mathrm{0}h}=$1). In addition, C_M(C_H) decreases with stability, and it decreases much slower when z_0h is not equal to z_0m.

Figure 6 Comparison of calculated and observed fluxes at Gucheng station from 1 December 2016 to 9 January 2017. (a) Momentum fluxes (MM5: z₀=0.0419); (b) sensible heat fluxes (MM5: z₀=0.0419); (c) sensible heat fluxes (MM5: z₀=0.0042). Red dots: the Li scheme; green plus signs: the MM5 scheme.

4.2 Comparison of momentum and sensible heat fluxes calculated by the two schemes

Using the obtained roughness lengths and the observations, the momentum and sensible heat fluxes were calculated by the Li and MM5 schemes. Firstly, z_0m and z_0h were set as 0.0419 and 0.0042, respectively, in the Li scheme, and z₀ was equal to z_0m in the MM5 scheme, to calculate the momentum and sensible heat fluxes, and the results are shown in Fig. 6a and b. It can be seen that compared with MM5, Li performs better with a higher regression coefficient and determination coefficient. For the momentum fluxes, the regression coefficient by Li is 0.6795 and that by MM5 is 0.5598, indicating that the error of Li is 12 % lower than that of MM5. For sensible heat fluxes, the regression coefficient by Li is 0.7967 and that by MM5 is 1.7994. The latter is much larger than 1; that is, the MM5 scheme obviously overestimates the sensible heat due to it not distinguishing z_0h from z_0m. Then, z₀ is made equal to 0.0042 in the MM5 scheme to recalculate the sensible heat flux, and the result is shown in Fig. 6c. It can be seen that the result has a great improvement after modifying the z₀ value, and the regression coefficient by MM5 is 0.7363, indicating that the error was reduced by 54 % after considering the z_0h effect. The result indicates that z_0h plays a critical role in both the SL scheme and the sensible heat flux (Chen and Zhang, 2009; Chen et al., 2011). However, the error of MM5 is still 6 % larger than that of Li. This illustrates that in addition to the effect of roughness length, the algorithm of the Li scheme itself is more reasonable than that of the MM5 scheme.

Figure 7 Variations of hourly turbulent fluxes and observed PM_2.5 at Gucheng station in daytime. (a) Momentum fluxes τ (blue line: observations; red line: the Li scheme; green line: the MM5 scheme) and PM_2.5 concentration (black line); (b) sensible heat fluxes H (the same as τ) and PM_2.5 concentration (black line). Yellow box: stage 1; blue box: stage 2; purple box: stage 3.

4.3 The specific performance of the two schemes in severe haze pollution

There were two obvious pollution processes during this observation period, and one occurred during 13 to 23 December 2016. Figure 7 shows the variations of hourly observed PM_2.5 concentration as well as the momentum and sensible heat fluxes calculated by the Li and MM5 schemes at Gucheng station in this process. For research purposes, only the daytime (from 08:00 to 20:00) was taken into account. Note that in MM5, z₀ was 0.0419 when calculating momentum fluxes and it was 0.0042 when calculating sensible heat fluxes. As shown in Fig. 7, the calculated results of momentum and sensible heat fluxes by the two schemes are generally consistent with the trend of the observations. Specifically, for the momentum fluxes (Fig. 7a), the results of two schemes have little difference when the values of observed momentum fluxes are large or at the peak. When the observed momentum fluxes are small, Li results are close to or less than the observations, while MM5 results are always higher than observations because of the limit of $u_{\ast }=\mathrm{0.1}$ in this scheme. For the sensible heat flux (Fig. 7b), MM5 results are always lower, while Li results are closer to observations, especially when the observed values are small. Furthermore, according to the evolution of PM_2.5 concentration, this haze event was then divided into three stages: the clear stage (stage 1: 13–14), the transition stage (stage 2: 16–18), and the maintenance stage (stage 3: 21–22). As shown in Fig. 7, in the clear stage (stage 1), the atmospheric stratification is unstable, PM_2.5 concentration is low, and there is a strong flux transport in the SL; the corresponding observations of the momentum and sensible heat fluxes are relatively high and they vary greatly. In the transition stage (stage 2), the atmosphere is changing from being unstable to stable, corresponding to haze formation; the momentum and sensible heat fluxes gradually decrease; and the daily variation also decreases. In the maintenance stage (stage 3), the atmospheric stratification is very stable, flux transport in the SL is weak, and both the momentum and sensible heat fluxes are at a low level. It can be seen that the Li results are generally closer to the observations compared with MM5 results in all three stages.

Figure 8 Probability distribution functions (PDFs) of the differences between calculated fluxes (momentum fluxes: left; sensible heat fluxes: right) using two schemes (the Li scheme: red bars; the MM5 scheme: green bars) and for observations in different stages (a, b: whole process; c, d: stage 1; e, f: stage 2; g, h: stage 3).

Figure 8 shows the probability distribution functions (PDFs) of the difference between calculated fluxes (using the Li and MM5 schemes) and observations in different stages at Gucheng station. In the whole pollution process, for the momentum fluxes (Fig. 8a), the PDF from Li tends to cluster in a narrower range centered by 0, and the probability within ±0.005 N m⁻² is 46.82 %, while this value from MM5 falls to 23.02 %. For the sensible heat fluxes (Fig. 8b), the PDF from Li is also more concentrated around 0 than that from MM5. The probabilities of bias from Li and MM5 within ±2.5 W m⁻² are 32.54 % and 13.49 %, respectively. In stage 1, for the momentum fluxes (Fig. 8c), the probability of bias from Li within ±0.005 N m⁻² is 38.09 %. The bias from MM5 is mainly larger than 0, and the probability within ±0.005 N m⁻² is 14.29 %. For the sensible heat fluxes (Fig. 8d), the probability of bias from Li within ±2.5 W m⁻² is 38.09 %, the same as momentum fluxes. The bias from MM5 is mainly less than 0, and the probability within ±2.5 W m⁻² is 9.52 %. In stage 2, the differences between the two schemes are more obvious. The PDFs from Li are the most concentrated around 0 in all cases, while those from MM5 are similar to stage 1. Specifically, for the momentum fluxes (Fig. 8e), the probabilities of bias from Li and MM5 within ±0.005 N m⁻² are 56.25 % and 25.00 %. For the sensible heat fluxes (Fig. 8f), the values within ±2.5 W m⁻² are 40.62 % and 6.25 %. In stage 3, the difference between the two schemes is small. For the momentum fluxes (Fig. 8g), the probabilities of bias from Li and MM5 within ±0.005 N m⁻² are 22.73 % and 27.27 %. For the sensible heat fluxes (Fig. 8h), the values from Li and MM5 within ±2.5 W m⁻² are both 36.36 %.

Table 2 Statistics between the Li and MM5 schemes of the calculated turbulent flux at Gucheng station.

τ: momentum flux; H: sensible heat flux; MB: mean bias; NMB: normalized mean bias; NME: normalized mean error; RMSE: root mean square error. The units of MB and RMSE are µg m⁻³.

Figure 9 As in Fig. 7 but for Beijing station.

Mean bias (MB), normalized mean bias (NMB), normalized mean error (NME), and root mean square error (RMSE) were calculated to test the results of the two schemes. Table 2 shows that the Li scheme generally estimates better than the MM5 scheme. In the whole haze process, the Li scheme underestimates the momentum fluxes by 3.63 % relative to the observations, while the MM5 scheme overestimates the momentum fluxes by 34.03 %. The Li and MM5 schemes underestimate the sensible heat fluxes by 15.69 % and 50.22 %, respectively. In the three stages, the Li scheme performs much better than the MM5 scheme in stage 1 and stage 2; especially in stage 2 when atmospheric stratification transforms from unstable to stable conditions, the difference between the Li and MM5 schemes is particularly significant. That is, the Li and MM5 schemes overestimate the momentum fluxes by 7.68% and 45.56 %, respectively, and they underestimate the sensible heat fluxes by 33.84 % and 76.88 %. The error of Li is much less than that of MM5. In view of the important role of atmospheric stratification in the generation and accumulation of PM_2.5 in stage 2, the Li scheme is expected to show better performance in the online simulation of PM_2.5 than MM5.

Figure 10 The mean momentum and sensible heat fluxes calculated using two schemes (a, b: the Li scheme; c, d: the MM5 scheme) and their differences (Li minus MM5. e: momentum fluxes; f: sensible heat fluxes) in Jing-Jin-Ji during the haze episode (13 to 23 December 2016).

Based on the good behavior of the Li scheme in Gucheng, the same experiment was performed at Beijing station to discuss the effect of different land-cover types on flux calculation. For Beijing station, the assumption z_0m=1 m, $z_{\mathrm{0}m}/z_{\mathrm{0}h}={\mathrm{10}}^{\mathrm{6}}$ was made to represent the surface conditions of the megacity due to a lack in situ measurements of surface turbulent flux. As shown in Fig. 9, the evolution of PM_2.5 concentration at Beijing station was also divided into three stages (stage 1: 13–15; stage 2: 17–19; stage 3: 20–21), like Gucheng shown in Fig. 7. Compared with Gucheng, the momentum transfer at Beijing station is obviously larger due to a great increase of the urban aerodynamic roughness length (z_0m). Meanwhile, the difference between Li and MM5 is greater at Beijing station. The sensible heat transfer of the Li scheme shows a great difference between clear days and pollution days; that is, the sensible heat transfer changes acutely in stage 1, while it changes smoothly in stage 2 and stage 3. However, the result of the MM5 scheme is significantly different from the Li result due to MM5 ignoring the z_0m effect, and the small number of z_0h keeps the sensible heat fluxes at a low level in all three stages.

To quantify the difference between the two schemes, a relative difference is defined as a percentage:

$$\begin{array}{}\text{(21)}& \mathrm{\Delta }V=\left|{\displaystyle \frac{V_{\mathrm{Li}}-V_{\mathrm{MM}\mathrm{5}}}{V_{\mathrm{MM}\mathrm{5}}}}\right|\times \mathrm{100}\mathit{\%},\end{array}$$

where V_Li and V_MM5 are the momentum (or sensible heat) fluxes calculated by the Li and MM5 schemes, respectively. We obtained the relative differences at the two stations in the three stages through the statistics. It is clear that the largest relative difference at Gucheng station is in stage 2 and at Beijing station it is in stage 1. The differences in Beijing are always larger than those in Gucheng for each of the three stages. Specifically, the relative differences of the momentum flux in stage 1, stage 2, and stage 3 increase by 73 %, 34 %, and 27 %, respectively, and the results of the sensible heat flux are 289 %, 52 %, and 68 %, respectively.

We further estimated the surface fluxes in the whole Jing-Jin-Ji region using the two schemes. Figure 10 shows the mean momentum and sensible heat fluxes calculated by the Li and MM5 schemes and their differences in Jing-Jin-Ji during the pollution episode. The assumption (z_0m=0.1 m, $z_{\mathrm{0}m}/z_{\mathrm{0}h}={\mathrm{10}}^{\mathrm{3}}$) was used to represent the average conditions of the underlying surface of the Jing-Jin-Ji region. As shown in Fig. 10, the momentum fluxes calculated by Li are less than those by MM5 in most stations; the sensible heat fluxes calculated by Li are usually larger than those by MM5. The result is consistent with the experiment at Gucheng station, which further indicates the importance of considering both z_0m and z_0h.

5 Conclusions

Using the observed momentum and sensible heat fluxes, together with conventional meteorological data including pressure, temperature, humidity, and wind speed from 1 December 2016 to 9 January 2017, including a severe pollution episode from 13 to 23 December 2016, the differences between the Li and MM5 schemes and the specific performances of the two were discussed and evaluated in this paper. The evolution process of atmospheric stratification from unstable to stable conditions corresponding to PM_2.5 accumulation was mainly discussed. The contributions of roughness lengths (z_0m and z_0h) as well as other factors in the SL schemes to the flux calculation for the momentum and sensible heat were also discussed in detail. The results are summarized as follows:

1. z_0m and z_0h have important effects on turbulent flux calculation in the SL schemes. Different values of z_0m and z_0h could induce great changes in the flux calculation, indicating that it is very necessary and important to distinguish z_0h from z_0m. Ignoring the difference between the two in the MM5 scheme led to large error in the calculation of the sensible heat flux, and this error in Gucheng was 54 %. Besides the roughness length, the algorithms in schemes are also important factors. In addition, ignoring the effect of the RSL in schemes may also result in certain bias of momentum and sensible heat fluxes in megacity regions which represent the rough underlying surface.

2. The effect of z_0m∕z_0h on turbulent fluxes is closely related to land-cover types (z_0m). A rough land-cover type (large z_0m) should be accompanied by a large value of z_0m∕z_0h. The differences between the two schemes for the momentum and sensible heat fluxes in Beijing were much larger than those in Gucheng. This suggests that the MM5 scheme probably induces greater error in megacities with a rough surface (e.g., Beijing) than in suburban areas with a smooth surface (e.g., Gucheng) due to the irrational algorithm of the MM5 scheme itself and the fact that the difference between z_0m and z_0h is ignored.

3. The Li scheme generally performed better than the MM5 scheme in the calculation of both the momentum flux and the sensible heat flux at Gucheng station. The Li scheme described atmospheric stratification, which is closely related to haze pollution, better compared with the MM5 scheme. This advantage was the most prominent in the transition stage from unstable to stable atmospheric stratification, corresponding to the PM_2.5 accumulation. In this stage, the momentum flux calculated by Li was overestimated by 7.68 % and this overestimation by MM5 was up to 45.56 %; the sensible heat flux by Li was underestimated by 33.84 % while this underestimation by MM5 was as much as 76.88 %. In most of the Jing-Jin-Ji region, the momentum fluxes calculated by Li were less than those by MM5 and the sensible heat fluxes by Li were larger than those by MM5, which is consistent with Gucheng.

The offline study of the two SL schemes in this paper showed the superiority of the Li scheme for surface flux calculation corresponding to the PM_2.5 evolution during haze episodes in Jing-Jin-Ji in eastern China. The study results offer a prerequisite and a possible way to improve PBL diffusion simulation and then PM_2.5 prediction, which will be achieved in the follow-up work of integrating the Li scheme into atmosphere chemical models.