An improved non-iterative surface layer flux scheme for atmospheric stable stratification conditions

Parameterization of turbulent fluxes under stably stratified conditions has always been a challenge. Current surface fluxes calculation schemes either need iterations or suffer low accuracy. In this paper, a non-iterative scheme is proposed to approach the classic iterative computation results using multiple regressions. It can be applied to the full range of roughness status 10 ≤ z/z0 ≤ 105 and−0.5 ≤ log(z0/z0h)≤ 30 under stable conditions 0 < RiB ≤ 2.5. The maximum (average) relative errors for the turbulent transfer coefficients for momentum and sensible heat are 12 % (1 %) and 9 % (1 %), respectively.


Introduction
In weather or climate models, the earth's surface is the boundary that needs to be resolved physically (Chen and Dudhia, 2001).The condition of atmosphere aloft (e.g., wind, temperature and humidity) is highly dependent on the momentum, sensible heat and latent heat fluxes at surface.Currently, the exchanges of momentum and heat between the earth's surface and the atmosphere are usually calculated with various schemes based on Monin-Obukhov similarity theory (hereinafter MOST; Monin and Obukhov, 1954) in models.These schemes (e.g., Paulson, 1970;Businger et al., 1971;Dyer, 1974;Holtslag and De Bruin, 1988;Beljaars and Holtslag, 1991;Janjić, 1994;Launiainen, 1995;Högström, 1996) are similar to each other, but the differences among them exist due to different observational data and/or mathematical solutions that were used in retrieving the schemes.One commonly used scheme is Businger-Dyer (BD) equation (Businger, 1966;Dyer, 1967).However, the BD equation suppresses fluxes under stable conditions too quickly and is not applicable when the Richardson number exceeds a critical value (Louis, 1979).Holtslag and De Bruin (1988) and Beljaars and Holtslag (1991) proposed alternative schemes that can be used under very stable conditions.With data collected in the field program CASES-99 (Cooperative Atmosphere-Surface Exchange Study-99) (Poulos et al., 2002), Cheng and Brutsaert (2005, CB05 hereinafter) further provided a new scheme, and it is confirmed to perform better by later research (Guo and Zhang, 2007;Jiménez et al., 2012).Based on the measurements made during experiment SHEBA in Arctic and Halley 2003 experiment in Antarctica, Grachev et al. (2007) and Sanz Rodrigo and Anderson (2013) proposed different similarity functions, respectively.Through systematic mathematical analysis, Sharan and Kumar (2011) proved that similarity functions of CB05 and Grachev et al. (2007) were applicable in the whole stable stratification region.However, all of these studies are based on MOST, and application of MOST in very stable conditions is in doubt since it assumes that turbulence is continuous and stationary, while in very stable conditions turbulence is weak, sporadic and patchy (Sharan and Kumar,  2011).Grachev et al. (2013) indicate that the applicability of local MOST in stable conditions is limited by the inequalities, when both gradient and flux Richardson numbers are below their "critical values", about 0.20-0.25.Further, MOST predicts that mean gradients of turbulence become independent of z (height) in very stable conditions.Wyngaard and Coté (1972) first referred to this limit as "z-less stratification".BD equations follow this prediction, but CB05 and Grachev et al. (2007) do not.To avoid these holdbacks and self-correlation of MOST, Sorbjan (2010) and Sorbjan and Grachev (2010) discussed an alternative local scaling for the stable boundary layer (referred to as gradient-based scaling) when different universal functions were plotted versus the  Another critical issue regarding the flux calculation with MOST is the numerical iteration.Under unstable conditions, the iteration normally converges within five steps (Fairall et al., 1996).By taking advantage of a bulk Richardson number parameterization for an improved first guess (Grachev and Fairall, 1997), the iteration can be reduced to three steps (Fairall et al., 2003).In the Weather Research Forecasting (WRF) model (Skamarock et al., 2008) MM5 similarity surface module, the flux variables from the previous time step are used to calculate the fluxes at current time step, and such an approach can yield reasonable results (Jiménez et  al., 2012).On the other hand, under stable conditions, the flux calculation takes many more steps to converge and hence is time-consuming.To avoid the iteration process, a series of non-iterative schemes are proposed (e.g., Louis, 1979;Garratt, 1992;Launiainen, 1995;Song, 1998;De Bruin et al., 2000;Yang et al., 2001;Li et al., 2010), but they all fail to cover the full range of −0.5 ≤ kB −1 ≤ 30, 10 ≤ z/z 0 ≤ 10 5 and −5.0 ≤ Ri B ≤ 2.5, which is pointed out by Wouters et al. (2012, WRL12 hereinafter).Here kB −1 = ln(z 0 /z 0h ).z is the reference height; and z 0 and z 0h are the aerodynamic and thermal roughness lengths, respectively.Ri B is the bulk Richardson number.Following WRL12, the condition that Ri B > 2.5 is not considered in this study, because it represents extremely stable stratification with very weak wind and little flux exchange.To calculate fluxes under all conditions, and also to include the roughness sublayer effect, WRL12 proposed an updated scheme based on the iterated results of CB05 under stable conditions.However, for a given Ri B , WRL12 uses only one equation to cover the whole large range of z/z 0 and kB −1 , which results in biases at some z/z 0 and kB −1 conditions.Therefore, to avoid the iteration process and keep the accuracy at the same time, the objective of this paper is to propose a group of equations that divide the calculation into eight regions according to z 0 and z 0h values.To compare with WRL12, and with the fact that CB05 equations are currently widely accepted, the new equations are also based on the iterated results of CB05 equations.WRL12.Section 3 introduces the new equations, and Sect. 4 intercompares these schemes.Summary and conclusions are presented in Sect. 5.

Revisiting CB05 and WRL12
The momentum flux τ and sensible heat flux H are defined as Here u * is the friction velocity, θ * is the temperature scale, ρ the air density, and c p the specific heat capacity at constant pressure.Based on MOST, the friction velocity u * and Here u and θ are the wind speed and potential temperature at the reference height z. k is the von Kármán constant.z * is the roughness sublayer height.θ 0 is the potential temperature at the height of z 0h .ψ m and ψ h are the integrated stability functions for momentum and heat, respectively.ψ * m and ψ * h are the correction functions accounting for roughness sublayer effect.L is the Obukhov length defined as ψ * m and ψ * h are given by De Ridder (2010): µ m = 2.59, µ h = 0.95, and φ m,h are the stability functions for momentum and heat.Following Sarkar and De Ridder (2010) and WRL12, z * /z 0 = 16.7 is adopted in this study.
CB05 gives the form of φ m,h and ψ m,h :   3), ( 4), ( 5), and (6), φ m,h and ψ m,h of CB05, fluxes can be calculated through iterations: with a first guess of ζ , u * and θ * can be calculated from Eqs. ( 3) and (4); then ζ again can be derived from Eq. ( 5).This procedure iterates until the results converge.The relationships of ζ ∼ Ri B , ζ ∼ ln(z/z 0 ), and ζ ∼ ln(z 0 /z 0h ) from CB05 are shown in Figs. 1, 2 and 3, respectively.Conditions with Ri B = 0.05, 0.2, 0.5, z/z 0 = 10, 1000, 10 5 and kB −1 = −0.5, 15, 30 are plotted.However, due to the limitation of computational time in numerical weather and climate models, the calculation results after five steps are always taken to approximate the fluxes (e.g., MYJ and MYNN surface module in WRF model; Janjić, 1996;Nakanishi and Niino, 2006).It is found that with the first guess of ln(z/z 0h ) and five steps of iteration, the results are still far away from the precise value.Figure 4 , and here n indicates the iteration step).Figure 5 shows the steps needed to converge into 5 % relative error with CB05 equations for various Ri B with z/z 0 = 10, 1000, 10 5 and kB −1 = −0.5, 15, 30.It shows that when Ri B = 0.74, z/z 0 = 10 and kB −1 = 30, more than 80 steps of iteration are needed to reduce the calculation error within 5 %.The iteration takes more steps to converge when there is a larger aerodynamic roughness length z 0 and a smaller thermal roughness length z 0h , which is common over an urban surface (Sugawara and Narita, 2009).When z/z 0 = 10 and kB −1 = 30, the largest error can reach 75 % after five-step iteration (Fig. 4) and 82 steps are needed for the results to converge (Fig. 5).However, when z/z 0 becomes large, for example z/z 0 = 10 5 (i.e., a representative value for a smooth sea surface), five steps are enough for the results to be within 5 % error under all kB −1 and Ri B conditions (Fig. 5).

Derivation of the new scheme
It can be seen from Figs. 1, 2 and 3 that ζ varies with Ri B , log(z/z 0 ) and kB −1 with remarkable nonlinearity.Specially, when kB −1 is large, ζ ∼ z 0 relationship can hardly be approximated by a cubic equation at some Ri B values (Fig. 2).Correspondingly, when z 0 is large, ζ ∼ z 0h also needs a high power series equation to approximate (at least cubic fit is not enough, Fig. 3).In order to reduce the complexity, weakly and strongly stable conditions were treated separately in previous studies (e.g., Launiainen, 1995;Li et al., 2010;WRL12).Analogously, multiple regions are considered for z 0 and z 0h for the regression of ζ = f Ri B , L 0M , kB −1 in this paper.In this way, the complexity of the equations can be reduced and at the same time their accuracy can be maintained.Although the total number of equations is increased due to the division of z 0 and z 0h , the calculation efficiency is still enhanced since the logical judgment of the region according to z 0 and z 0h values in program codes takes much less time than iterations.The critical issue here is how to divide the z 0 and z 0h regions in a reasonable way to obtain the smallest number of regions but the highest accuracy.For this purpose, the z 0 and z 0h are first divided into 13 and 14 sections according to the values of z/z 0 and z 0 /z 0h , respectively.For z/z 0 , the sections are 10 ∼ 20, 20 ∼ 40, 40 ∼ 80, . . ., 10 240 ∼ 20 480, 20 480 ∼ 40 960 and 40 960 ∼ 10 5 ; for z 0 /z 0h , the sections are 0.607 ∼ 1, 1 ∼ 10, 10 ∼ 100, 100 ∼ 10 3 , 10 3 ∼ 10 4 , . . ., 10 11 ∼ 10 12 and 10 12 ∼ 1.07 × 10 13 .z/z 0 ∈ 10 ∼ 20 and z 0 /z 0h ∈ 10 12 ∼ 1.07×10 13 is the region that needs the highest power series equation to approximate.This region is firstly chosen to find a maximum critical value of ζ c1 that can make the regression Then some of the z 0 and z 0h regions can be merged with each other for the section ζ ∈ 0 ∼ 0.33 and a total of 8 z 0 − z 0h regions are left in the z 0 − z 0h plane.In other words, the regression error of Eq. ( 23) can be kept within 5 % in any of the eight regions when ζ ∈ 0 ∼ 0.33 (Table 1).Thus, for these eight regions, it can be found that with the sections divided by the specified critical values ζ cp (where p is 1, 2, 3, . . ., it indicates the section and its maximum value depends on the z 0 − z 0h region), the regression error with Eq. ( 23) can be kept within 5 % for ζ ≤ 0.5 and 10 % or ζ > 0.5.For a given pair of z 0 and z 0h , the division by ζ cp can be transformed to Ri Bcp : Here m, n = 0, 1, 2, and m + n ≤ 3; p is 1, 2, 3, . . ., which indicates the section and its maximum value depends on the z 0 − z 0h region.For region 1 and 7, the maximum p is 6, while for other regions it varies between 3 and 5.The coefficients for Eq. ( 24) are shown in Table 2.The Ri Bcp then cuts the 0-2.5 Ri B range into several sections: Sect. 1 is from 0 to Ri Bc1 , Sect. 2 from Ri Bc1 to Ri Bc2 , and so on.The coefficients for Eq. ( 23) in each section are given in Tables 3-10.The procedure to obtain these coefficients is summarized below: 1  Method: when ζ ∈ 0 ∼ ζ c1 , regression with Eq. ( 23) is kept within 5 % error.
Method: variations of combinations of the 13 sections of z/z 0 and 14 sections of z 0 /z 0h are tested to minimize the numbers of regions, and regression with Eq. ( 23) and ζ ∈ 0 ∼ 0.33 is kept within 5 % error.
The calculation procedure for a given group of z 0 , z 0h and Ri B is the following: (1) find the region according to z 0 and z 0h with Table 1; (2) find the section according to the region and Ri B with Eq. ( 24) and coefficients in Table 2; and (3) in Tables 3-10 find the coefficients for the particular region and section and use Eq. ( 23) to calculate ζ .Figure 7 presents the relative error of ζ with new equations compared with iterated results of CB05 for various Ri B with z/z 0 = 10, 1000, 10 5 and kB −1 = −0.5, 15, 30.With the new equations, the relative error is controlled to be within 10 % for the whole range.Especially, when ζ ≤ 0.5, the relative error is within 5 % since it happens more often in the real conditions (Fig. 8). .
Here "Error (ζ )" indicates ζ or C M,H at a particular ζ , z 0 and z 0h .Although Eq. ( 29) presents the form of continuous integral, it is actually calculated discretely with interval 0.035 for log( z z ) and 0.1 for log( z 0 z 0h ).The results indicate that the maximum ζ exceeds 50 % when using CB05 with five-step iteration or WRL12, while the averaged ζ for the two methods both exceeds 15 %.On the contrary, the maximum ζ of the new scheme is always smaller than 5 % (when ζ ≤ 0.5) and 10 % (when ζ > 0.5), and the average ζ is always smaller than 2 % in the whole range.The maximum C M from CB05 with five-step iteration (WRL12) exceeds 50 % (40 %), and average C M exceeds 30 % (8 %).The maximum C H from CB05 with fivestep iteration (WRL12) exceeds 50 % (24 %), and average C H exceeds 18 % (6 %).Comparatively, the new scheme controls the maximum C M ( C H ) to be within 12 % (9 %) and the average C M ( C H ) within 1 % (1 %).Table 11 summarizes the characteristics of the four methods.

Summary and conclusions
Although CB05 provides a way to calculate surface fluxes under stable conditions, its practical usage is confined due to the involved iteration process.It has been shown that iteration with five steps will result in large calculation errors, especially when z/z 0 is small and kB −1 is large, which is common over an urban surface.WRL12 proposed a way to avoid the iteration, but it introduces large error in the calculation procedure so that its calculation accuracy needs to be improved.Through dividing the z 0 − z 0h plane into eight regions, the new scheme develops a group of equations with higher accuracy.The calculation error of ζ = f Ri B , L 0M , kB −1 is always controlled to be within 5 % (when ζ ≤ 0.5) and 10 % (when ζ > 0.5).The calculation procedure is also simple; for a small Ri B (i.e., Ri B < Ri Bc1 ), only one time computation of Eqs. ( 23) and (24) will suffice.The maximum computation step is six times for Eq. ( 24) and one time for Eq. ( 23) when it is in region 1 or 7 and at the same time Ri B is large (i.e., Ri B > Ri Bc6 ).Note that the Eq. ( 24) has only a maximum of eight elements and a minimum of four elements, so the calculation is still efficient.
The new equations involve a large number of parameters, which increase the complexity of coding.However, the effort of coding the new scheme is minimal as compared to its potential gain, which includes the accuracy of the new scheme and the avoidance of iterations.Besides, a compromise can be made between accuracy and complexity.For models that are not interested in high kB −1 values, region 1 and 2 (i.e., 10 ≤ z/z 0 ≤ 10 5 and −0.607 ≤ z 0 /z 0h ≤ 100) have provided reasonable coverage (see Garratt, 1992;Launiainen, 1995), and the other six regions can be ignored.For example, in WRF model MM5 surface module, z 0h = z 0 is assumed during the calculation of frictional velocity (Jiménez et al., 2012).While for models that include urban surface effects, it is better to keep all the regions.Further, CB05 probably is not the final solution for the surface flux calculation under stable stratification.The method used to derive non-iterative equations presented here can be used in future studies to transfer the new iterative algorithm to non-iterative equations.Overall, the new equations cover the full range of −0.5 ≤ kB −1 ≤ 30, 10 ≤ z/z 0 ≤ 10 5 and stable conditions (i.e., 0 < Ri B ≤ 2.5), and maintain high accuracy and efficiency.It is expected that its usage in climate and weather forecasting models can lead to better performance in surface flux calculation under stable conditions, especially over urban surfaces.

Fig. 1 .
Fig. 1.The relationship between Ri B and ζ from the precise results of CB05.

Fig. 2 .
Fig. 2. The relationship between log(z/z 0 ) and ζ from the precise results of CB05.Black line indicates the cubic fit of curve kB −1 = 30 and Ri B = 0.5 with least square method.

Fig. 3 .
Fig. 3.The relationship between log(z/z 0h ) (i.e., kB −1 ) and ζ from the precise results of CB05.Black line indicates the cubic fit of curve z/z 0 = 10 and Ri B = 0.5 with least square method.

Fig. 4 .
Fig. 4. Relative error after five steps of iteration with CB05 equations under certain z 0 and z 0h conditions.

Fig. 5 .
Fig. 5. Steps needed to converge into 5 % relative error with CB05 equations under certain z 0 and z 0h conditions.The inset shows the whole perspective.

4
Comparison of the results from CB05 with five-step iteration, WRL12 and the new scheme The maximum and average relative error of ζ , C M and C H calculated from CB05 with five-step iteration, WRL12 and the new scheme are shown in Figs. 8, 9 and 10 for various ζ with z/z 0 = 10, 1000, 10 5 and kB −1 = −0.5, 15, 30.C M and C H are the transfer coefficients for momentum and sensible heat respectively, and

Table 1 .
The eight regions divided by z/z 0 and z 0 /z 0h values.
where ζ (cal) is the calculation result, and ζ (precise) is the precise result from the ultimate iteration of CB05 (when

Table 4 .
Similar to Table 3, but for region 2.

Table 5 .
Similar to Table3, but for region 3.

Table 6 .
Similar to Table 3, but for region 4.

Table 7 .
Similar to Table 3, but for region 5.

Table 8 .
Similar to Table3, but for region 6.

515-529, 2014 526 Y. Li et al.: An improved non-iterative surface layer flux schemeTable 9 .
The relative error for C M and C H is calculated from Similar to Table3, but for region 7.

Table 10 .
Similar to Table 3, but for region 8.

Table 11 .
Summarization of the characteristics of the four methods.Calculation time is the time each method needs for computing ζ from Ri B , z 0 and z 0h in the range 0 < Ri B ≤ 2.5, 10 ≤ z/z 0 ≤ 10 5 and −0.5 ≤ log(z 0 /z 0h ) ≤ 30 with the interval of 0.01 for Ri B , 0.035 for log(z/z 0 ) and 0.1 for log(z 0 /z 0h ).The calculation is performed on a desktop computer with an Intel Core i5 processor, and note that the calculation time can vary with different computer.